Parallel Direct Solver Research Articles

NUMA-aware parallel sparse LU factorization for SPICE-based circuit simulators on ARM multi-core processors

In circuit simulators that resemble the Simulation Program with Integrated Circuit Emphasis (SPICE), one of the most crucial steps is the solution of numerous sparse linear equations generated by frequency domain analysis or time domain analysis. The sparse direct solvers based on lower-upper (LU) factorization are extremely time-consuming, so their performance has become a significant bottleneck. Despite the existence of some parallel sparse direct solvers for circuit simulation problems, they remain challenging to adapt in terms of performance and scalability in the face of rapidly evolving parallel computers with multiple NUMA hardware based on ARM architecture. In this paper, we introduce a parallel sparse direct solver named HLU, which re-examines the performance of the parallel algorithm from the viewpoint of parallelism in pipeline mode and the computing efficiency of each task. To maximize task-level parallelism and further minimize the thread waiting time, HLU devises a fine-grained scheduling method based on an elimination tree in pipeline mode, which employs depth-first search (DFS-like) to iteratively search for parent tasks and then place dependent tasks in the same task queue. HLU also suggests two NUMA node affinity strategies: thread affinity optimization based on NUMA nodes topology to guarantee computational load balancing and data affinity optimization to enable effective memory placement when threads access data. The rationality and effectiveness of the sparse solver HLU are validated by the SuiteSparse Matrix Collection. In comparison with KLU and NICSLU, the experimental results and analysis show that HLU attains a speedup of up to 9.14× and 1.26x (geometric mean) on a Huawei Kunpeng 920 Server, respectively.

On cost-efficient parallel iterative solvers for 3D frequency-domain seismic multisource viscoelastic anisotropic wave modeling

Solving large sparse linear systems in 3D frequency-domain seismic wave modeling, especially in viscoelastic anisotropic media, poses significant challenges due to the increasing number of discrete moduli and nonzero elements in the linear system matrix. The computational load surpasses that of acoustic or viscoacoustic media, making it even more challenging when dealing with multisource problems. Popular scientific tools for solving a linear system, such as multifrontal massively parallel sparse direct solver (MUMPS), structured matrix package (STRUMPACK), and portable extensible toolkit for scientific computation (PETSc), can be used, but their applicability to our specific problem has not been comprehensively evaluated. Our study aims to tackle the challenges in solving large, sparse, complex-valued symmetric linear systems with multiple right-hand side vectors for 3D frequency-domain seismic wave modeling. We have adopted the preconditioned conjugate gradient iterative algorithms as the foundation for our research, introducing two highly cost-effective parallel iterative solvers: the parallel symmetric successive overrelaxation conjugate gradient (P-SSORCG) and the parallel incomplete Cholesky conjugate gradient (P-ICCG). These novel solvers are subjected to a comprehensive comparative analysis against well-established scientific tools, such as MUMPS, STRUMPACK, and PETSc, in the context of 3D frequency-domain seismic wave modeling. We show their promising performances in a practical 3D SEG/EAGE overthrust model and demonstrate that the grouped P-SSORCG offers an efficient alternative to parallel direct solvers, particularly in situations wherein computational resources are limited.

A Novel Parallel Direct Solver Using Adaptive Cross Approximation for Analysis of Electromagnetic Radiation Problems With Complex Structures

A Novel Parallel Direct Solver Using Adaptive Cross Approximation for Analysis of Electromagnetic Radiation Problems With Complex Structures

A parallel preconditioning Schur complement approach for large scale industrial problems

The purpose of this paper is to introduce a new strategy to improve the convergence and efficiency of the class of domain decomposition known as Schur complement techniques related to interface variables for the simulation of mechanical, electrical and thermal problems in presence of cross points. More precisely, we are interested not only in domain decomposition with two equal parts having the same physical properties but rather in more general splitting components. It is known that in the first case, the optimal convergence with good pre-conditioner is obtained in two iterations and the problem is still challenging in the second case. The primary goal then is to fill part of the gap in such problem domain decomposition techniques and to contribute to solve extremely difficult industrial problems of large scale by using parallel sparse direct solver of the multi-core environment of the whole system and handling each part of the system independently of the change of the mesh or the shifting of the mathematical method of resolution and subsequently, we treat the interface as boundary conditions. The numerical experiments of our algorithm are performed on few models arising from discretization of partial differential equations using Finite Element Method (FEM).

Effective elastic properties of three-dimensional multiple crack problems with the isogeometric boundary element parallel fast direct solver

We propose an isogeometric boundary element (IGABEM) parallel fast direct solver based on OpenMP to solve three-dimensional multiple crack problems. In this paper, the orderly and randomly distributed cracks are studied to observe the influence of the multiple crack interaction on the effective elastic properties. This paper designs the parallel implementation of the matrix low rank approximation and the matrix assembly for the hierarchical off-diagonal low rank (HODLR) matrix. In the matrix low rank approximation, we propose the m-n tree to improve the accelerated hybrid algorithm to balance the CPU performance and memory usage. The numerical results show that the effective elastic properties can respond to the strong interaction to a certain extent, but the maximum stress intensity factors (SIFs) of multiple cracks still need to be calculated to diagnose the possible crack initiation points. On the other hand, the numerical results also show that the IGABEM parallel fast direct solver can greatly improve the computational efficiency on the premise of ensuring the computational accuracy.

Is 3D frequency-domain FWI of full-azimuth/long-offset OBN data feasible? The Gorgon data FWI case study

Frequency-domain full-waveform inversion (FWI) is potentially amenable to efficient processing of full-azimuth long-offset stationary-recording seabed acquisition carried out with a sparse layout of ocean-bottom nodes (OBNs) and broadband sources because the inversion can be performed with a few discrete frequencies. However, computing the solution of the forward (boundary-value) problem efficiently in the frequency domain with linear algebra solvers remains a challenge for large computational domains involving tens to hundreds of millions of parameters. We illustrate the feasibility of 3D frequency-domain FWI with a subset of the 2015/2016 Gorgon OBN data set in the North West Shelf, Australia. We solve the forward problem with the massively parallel multifrontal direct solver MUMPS, which includes four key features to reach high computational efficiency: an efficient parallelism combining message-passing interface and multithreading, block low-rank compression, mixed-precision arithmetic, and efficient processing of sparse sources. The Gorgon subdata set involves 650 OBNs that are processed as reciprocal sources and 400,000 sources. Monoparameter FWI for vertical wavespeed is performed in the viscoacoustic vertically transverse isotropic approximation with a classical frequency continuation approach proceeding from a starting frequency of 1.7 Hz to a final frequency of 13 Hz. The target covers an area ranging from 260 km2 (frequency ≥ 8.5 Hz) to 705 km2 (frequency ≤ 8.5 Hz) for a maximum depth of 8 km. Compared to the starting model, FWI dramatically improves the reconstruction of the bounding faults of the Gorgon horst at reservoir depths as well as several intrahorst faults and several horizons of the Mungaroo Formation down to a depth of 7 km. Seismic modeling reveals a good kinematic agreement between recorded and simulated data, but amplitude mismatches between the recorded and simulated reflection from the reservoir suggest elastic effects. Therefore, future works involve multiparameter reconstruction for density and attenuation before considering elastic FWI from hydrophone and geophone data.

Integration of a Gaussian quadrature grid discretization approach with a generalized stiffness reduction method and a parallelized direct solver for 3-D frequency-domain seismic wave modelling in viscoelastic anisotropic media

SUMMARY We integrate three advanced numerical techniques—Gaussian quadrature grid (GQG) discretization, a new generalized stiffness reduction method and the latest version of an efficient parallelized direct solver to achieve accurate 3-D frequency-domain seismic wave modelling in viscoelastic anisotropic media. A GQG is employed to sample and interpolate both model parameters and wavefield quantities as well as to fit with arbitrary free-surface topography and subsurface interfaces of a geological model. A new version of the generalized stiffness reduction method is utilized to effectively remove the artificial boundary edge effects for which the common perfectly matched layer method fails. The most recent version of a multifrontal massively parallel direct solver is applied to tackle the notoriously expensive computation of frequency-domain 3-D wave modelling. We validate the 3-D modelling by comparing with the exact solutions for homogeneous viscoelastic isotropic, vertically transversely isotropic and orthorhombic media. All the results show very close matches between the numerical and analytical solutions. Then, we investigate the computational efficiency of the parallelized direct solver, compare its performance using different ordering schemes, in-core and out-of-core factorization modes and the block low-rank approximation in the factorization for different grid sizes. Our modelling results show that the ordering scheme of the so-called ‘Metis’ is the best for reducing computer memory and run time, and the parallelized direct solver is remarkably faster than iterative solvers for similar workloads but at the expense of higher memory requirements. The out-of-core factorization mode can effectively reduce the memory cost without a compromising on run time. The block low-rank approximation is able to significantly reduce the run time in both the factorization and solving process (up to 56 per cent in total), but will increase the memory cost when using the out-of-core factorization mode. Efficient application of this parallel direct solver should use ‘Metis’ as the ordering scheme and select the out-of-core factorization mode without the block low-rank approximation as the best scheme to save the memory cost, or the in-core factorization mode with the block low-rank approximation for the fastest computation. Finally, we demonstrate the excellent applicability of the 3-D wave modelling scheme for a practical and complex heterogeneous geological model.

Fine-grained heterogeneous parallel direct solver for finite element problems

The improvement of computational effectiveness is a vital issue in the field of large-scale finite element analysis. The performance is fundamentally determined by the efficiency of solving sparse linear system equations using the implicit finite element method. This paper presents a direct linear solver based on heterogeneous hybrid parallel computing on CPUs and GPUs. This can efficiently utilize computing resources of multiple devices to achieve performance improvement. Initially, we partition the elimination tree into several subtrees to accomplish the task decomposition. Based on this, we build a dynamic programming mathematical model to balance the computational load of the various devices. Then, we develop a numerical decomposition strategy by combining node parallelism and tree parallelism for the CPUs. In addition, efficient numerical decomposition is achieved on the GPU through batch processing and maximizing the overlap between computations and data transfers. Numerical experiments show that, compared with MKL PARDISO, the performance of numerical factorization can be improved by up to 10 times by using CPU and dual-path GPU hybrid calculations, and the computation time of simulation can be reduced by one-third for the multicondition analysis of Body In White and by 20% for the large-scale nonlinear finite element deformation analysis.

Iterative solver with folded preconditioner for finite element simulation of magnetotelluric fields

We propose an iterative solver with a folded preconditioner for large sparse systems of linear algebraic equations (SLAEs), which are obtained as a result of a vector finite element approximation of electromagnetic fields in magnetotelluric problems. The solver is based on the computational scheme of the Conjugate A-Orthogonal Conjugate Residual method (COCR), diagonal preconditioner, and residual smoothing. To build a finite element approximation, the primary-secondary field approach is used. To construct a preconditioner, we use a matrix, the elements of which are the weights of the representation of gradients of nodal basis functions as a linear combination of edge basis functions. For Dublin Test Model 1 in a wide frequency range from 0.0001 Hz to 10 Hz, we compare the accuracy and computational efficiency of our code with the results of other authors, as well as compare the computational efficiency of the proposed diagonal folded preconditioned (DFP) COCR solver with other solvers. It is shown that the apparent resistivities obtained on the coarsest mesh with our code differ by no more than a few percent from the solution on a finer mesh and from some solutions obtained by other authors using the finite element, finite difference, and integral equation methods. The computational time required for our code is 8 times less than the minimum time presented by other authors (with a corresponding recalculation of the performance of the computers used). It is also shown that the DFP COCR solver practically does not increase the number of iterations at low frequencies and its computational time is about 20 times less than the time required for the COCR solver with a standard diagonal preconditioner. A comparison of the DFP COCR solver with the Parallel Direct Sparse solver (PARDISO) shows that even at low and medium frequencies, the computational time of the DFP COCR solver is about 20 times less than the computational time of the PARDISO, and at high frequencies this difference reaches already 50 times.

Effective 3D-transient electromagnetic inversion using finite-element method with a parallel direct solver

The correct interpretation of time-domain or transient electromagnetic (TEM) data relies on effective inversion techniques. Herein, we have developed an effective 3D-TEM inversion method based on the vector finite-element method. The domain was discretized with unstructured tetrahedra to simulate complex geologic structures. We develop a second-order backward-Euler scheme, rigorously with nonuniform time-step sizes, to approximate the time derivative of the electric field. The waveform was directly simulated and considered in the forward modeling and inversion. In a 3D inversion, the most time-consuming element of the workflow is solving a large system of linear equations. To alleviate this problem, we adopt a parallel direct solver for solving a large system of equations involved in forward modeling and sensitivity calculations. We use the Gauss-Newton approach to solve the inverse problem and formulate the least-squares problem to obtain model updates. We avoided the explicit sensitivity calculation by solving the matrix-vector multiplication between the sensitivity matrix (or its transpose) and vector. The developed inversion algorithm is validated using several realistic synthetic models. Finally, we apply the method to practical TEM data with complex topography for coal mine goaf detection. We thus obtain a detailed 3D resistivity model that revealed the layering structure and some localized anomalies related to the coal mine goaf.

Modeling of unsteady flows of multiphase viscous fluid in a pore space

The authors have developed and implemented a numerical algorithm to model unsteady flows of a viscous multiphase isothermal fluid by finite difference method using the projection method for the numerical solution of the Navier-Stokes equation. The projection method implies splitting the initial system of equations by physical processes, in which convective transport and the effect of the pressure gradient are separately taken into account. As a result, at each step, it is necessary to solve the Poisson equation to find the pressure field. The solution of SLAE is performed by a parallel direct solver based on LU decomposition. An explicit scheme is used to solve the Cahn-Hilliard equation to update the phase field, the parameter of which is taken into account when adding surface forces to the Navier-Stokes equation. Computational experiments showing qualitative and quantitative agreement with experimental and numerical data from the literature are presented.

Robust and scalable adaptive BDDC preconditioners for virtual element discretizations of elliptic partial differential equations in mixed form

The Virtual Element Method (VEM) is a recent numerical technology for the solution of partial differential equations on computational grids constituted by polygonal or polyhedral elements of very general shape. The aim of this work is to develop effective linear solvers for a general order VEM approximation designed to approximate three-dimensional scalar elliptic equations in mixed form. The proposed Balancing Domain Decomposition by Constraints (BDDC) preconditioner allows to use conjugate gradient iterations, albeit the algebraic linear systems arising from the discretization of the problem are indefinite, ill-conditioned, and of saddle point nature. The condition number of the resulting positive definite preconditioned system is adaptively controlled by means of deluxe scaling operators and suitable local generalized eigenvalue problems for the selection of optimal primal constraints. Numerical results confirm the theoretical estimates and the reliability of the adaptive procedure, with the experimental condition numbers always very close to the prescribed adaptive tolerance parameter. The scalability and quasi-optimality of the preconditioner are demonstrated, and the performances of the proposed solver are compared with state-of-the-art parallel direct solvers and block preconditioning techniques in a distributed memory setting.

Three-dimensional regularized inversion of magnetotelluric data with a minimum support gradient constraint

Regularization is an effective technique to obtain stable solutions for 3D magnetotelluric (MT) inversion. Smooth constraints are commonly used in traditional inversion algorithms. It is unsuitable if anomalies occur in block patterns. In some cases, it is necessary to obtain sharp boundaries. Some studies have proposed stabilizing functionals to obtain clear and sharp boundaries, such as minimum support (MS) and minimum gradient support (MGS) functionals; however, the optimal value of the focusing parameter should be selected to produce reliable inversion results. In this study, we expanded a new stabilizing functional that we called minimum support gradient (MSG) stabilizing functional as a stabilizer in 3D MT inversion. We developed a Gauss-Newton (GN) algorithm based on the edge-based finite-element method (FEM) to invert the MT data. In the inversion process, the forward modeling and the adjoint problem are solved by the Math Kernel Library (MKL) PARDISO parallel direct solver. we also parallelized the inversion program over frequencies by using Message Passing Interface (MPI) to further speed up the inversion process. Two synthetic models are tested to compare the inversion results for different stabilizing functionals. Compared to traditional smooth inversion, the MSG inversion can clearly distinguish the boundaries of the anomalies. The inverted resistive targets are closer to the true test model. Additionally, the MSG inversion results are stable and robust for a wide range of focusing parameter values, indicating that our method is effective and easy to implement. Finally, the proposed algorithm is applied to the long MT data from the USArray component of EarthScope project in the United States.

3-D Finite Volume Modeling for LWD Azimuthal Propagation Resistivity Tool With Multiple Annular Antenna Recesses Using Coupled Potentials on Cylindrical Grids

The sensors of logging-while-drilling (LWD) azimuth propagation resistivity (APR) tool are usually installed in the annular antenna recesses filled with ferrite or other magnetic materials on the grooved metal mandrel. In this article, we will improve the 3-D finite volume method (3-D FVM) on cylindrical grids to solve the Helmholtz equations about coupled potentials for precise simulation of the tool responses with multiple antenna recesses. The irregular and novel Yee’s staggered grids are designed to conform to the grooved mandrel surface. The electrical continuation technique is adopted to approximate the perfectly electrical conductor (PEC) boundary condition on the mandrel surface for efficient enhancement of the discretization precision of the potentials around the mandrel. The equivalent conductivity and permeability in heterogeneous elements are determined by the volume average. After that, the parallel direct sparse solver (PARDISO) is used to solve the discretized equation about coupled potentials for effectively ensuring the accuracy and stability of numerical solutions. Finally, we validate the 3-D FVM against results by the hybrid algorithm in homogeneous formation. Through numerical results, we further investigate the influence of both the mandrel and antenna recesses on the LWD tool responses, and its response characteristics in several different environments.

3-D Inversion of Z-Axis Tipper Electromagnetic Data Using Finite-Element Method With Unstructured Tetrahedral Grids

<inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$Z$ </tex-math></inline-formula> -axis tipper electromagnetic (ZTEM) technique is an airborne electromagnetic method that detects anomalies in the deep earth that are induced by natural sources. Conventional ZTEM forward modeling is generally conducted using structured grids that have limited accuracy and cannot be used to invert complex underground structures and topography. However, unstructured grids can accurately model complex underground structures with fewer cells. Therefore, to effectively model and recover the complex underground structures, we developed a 3-D framework with unstructured tetrahedral grids for ZTEM forward modeling and inversions. The forward problem was formulated using the finite-element method with unstructured tetrahedral grids. To solve the inverse problem, we used a limited-memory quasi-Newton algorithm (L-BFGS) with a parallel direct solver for optimization in order to avoid the explicit calculation of the Hessian matrix and save the memory and computational time. To validate our forward and inversion algorithms, we conducted numerical experiments on three synthetic models and inverted a survey dataset acquired for a tunnel project in Tibet, China. The experimental results demonstrate the effectiveness of our unstructured finite-element and L-BFGS method for modeling and inverting ZTEM data.

3-D Magnetotelluric Inversion and Application Using the Edge-Based Finite Element With Hexahedral Mesh

Three-dimensional (3-D) inversion technique has become an important and practical approach for magnetotelluric (MT) data interpretation. In this article, we developed a 3-D parallelized MT inversion scheme using the edge-based finite element method and applied the developed method to the newly collected MT data in the Xinjiang Luntai area. The distorted hexahedral element is adopted to incorporate topography into the forward modeling and inversion for complicated scenarios. We use the Gauss–Newton optimization method to minimize the objective functional for MT inversion. The developed algorithm is parallelized using MPI over frequencies and parallel direct solvers when solving the forward and adjoint problems for each frequency. We compare the performance of the least-square QR (LSQR) factorization and preconditioned conjugate gradient (PCG) solvers for the model update within each Gauss–Newton iteration and found that the LSQR solver is more stable. The developed inversion algorithm is validated using several synthetic models. Finally, we applied the inversion algorithm to the subsurface resistivity imaging in the Luntai area. The recovered geoelectric model from full 3-D inversion fits well with the known geological and geophysical information. The recovered model shows a low resistivity layer which may be caused by the salt strata. Besides, the inversion results reveal the movement tectonic in this survey area within a depth of 9 km.

Scalable domain decomposition preconditioner for Navier–Stokes equations coupled with the heat equation

In this article, we study the thermal instability that appears from time to time while printing using a 3D printer. To solve the semi-discretized problem at each time-step, we use a scalable parallel algorithm based on a two-level Optimized Restricted Additive Schwarz (ORAS) domain decomposition preconditioner for GMRES. Parallel scalability tests are conducted with comparison against the parallel direct solver MUMPS and the one-level Schwarz method, which show lack of robustness for larger number of processors. 2D numerical tests illustrate that the number of iterations to reach GMRES convergence depends on the state of the physical simulation during time, and that the second level of preconditioning is needed to achieve robustness.

3D MT Anisotropic Inversion Based on Unstructured Finite-element Method

The conventional 3D magnetotelluric (MT) forward modeling and inversions generally assume an isotropic earth model. However, wrong results can be obtained when using an isotropic model to interpret the data influenced by the anisotropy. To effectively model and recover the earth structures including anisotropy, we develop a 3D MT inversion framework for a triaxial anisotropic model. We use the unstructured finite-element method for our forward modeling. This offers more possibility to simulate more complex underground geology and topography. To solve the inverse modeling problem, we use a limited-memory quasi-Newton algorithm (L-BFGS) with a parallel direct solver for optimization that avoids the explicit calculation of the Hessian matrix and saves the memory and computational time. We validate our algorithm via numerical experiments on both synthetic data and MT survey data from the US Array project.

A Parallel Direct Domain Decomposition Solver Based on Schur Complement for Electromagnetic Finite Element Analysis

In this letter, a parallel direct solver for the analysis of electromagnetic problems with minimum memory cost is presented. The solver makes use of the finite element method to perform the analysis dividing the entire computational domain into several subdomains. Then, the partial lower/upper decomposition of each subdomain is performed obtaining the Schur complement matrix of each subdomain. This process is performed one by one using all the computer power available for each decomposition. Subsequently, the original 3-D problem is reduced to a 2-D problem on the interface skeleton through assembling the obtained Schur complement matrices. After obtaining the interface solution by using a dense direct solver, each well-posed subdomain is solved independently. Thus, the resulting parallel solver does not suffer from convergence issues since is based on direct method strategies instead of iterative methodologies. Through comparisons with commonly used parallel direct solvers, the benefits and the accuracy of the proposed solver are demonstrated. Finally, a waveguide wide-edge slot antenna with 14 million tetrahedrons is analyzed showing the power of the presented parallel direct solver.

Alternating directions parallel hybrid memory iGRM direct solver for non-stationary simulations

The three-dimensional isogeometric analysis (IGA-FEM) is a modern method for simulation. The idea is to utilize B-splines or NURBS basis functions for both computational domain descriptions and the engineering computations. Refined isogeometric analysis (rIGA) employs a mixture of patches of elements with B-spline basis functions, and $C^0$ separators between them. It enables a reduction of the computational cost of direct solvers. Both IGA and rIGA come with challenging sparse matrix structure, that is expensive to generate. In this paper, we show a hybrid parallelization method to reduce the computational cost of the integration phase using hybrid-memory parallel machines. The two-level parallelization includes the partitioning of the computational mesh into sub-domains on the first level (MPI), and loop parallelization on the second level (OpenMP). We show that hybrid parallelization of the integration reduces the contribution of this phase significantly. Thus, alternative algorithms for fast isogeometric integration are not necessary.