Parallel Sparse Direct Solver Research Articles

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 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).

Virtual battery electric vehicle development via 1D tools

Actual testing and prototyping costs make up a significant portion of engineering budgets. Virtual demonstration mainly relies on fast and accurate models with robust performance prediction capability as a cost-effective solution. In this manuscript, GT-Suite, a one-dimensional simulation tool is preferred for developing an isothermal battery model. The developed battery model is implemented to a Light Commercial Vehicle model to run Worldwide Harmonized Light Vehicles Test Cycle. The critical outputs such as state of charge, energy depletion, heat rejection, and generated power are reported. Histograms such as pack current, charge, and discharge current with respect to its nominal capacity (C-rate) are also created to examine the operation capability of battery. In the results, it is seen that the state of charge diminishes from 90% to 82%, as an expected behavior. It is also found that the Parallel Sparse Direct and Multi-Recursive Iterative Linear Solvers methodology reduces the simulation duration by approximately 100 times, in comparison with Singular Value Decomposition Electrical Inversion Scheme. The runtime of the battery pack modeled with the cellular approach combined with the SVD method is more than 90 h. However, the runtime drops to 0.75 h when the PARDISO technique is applied. The method developed within this study can be used for rapid and accurate development of batteries. The verification can be completed in virtual environment and a vital reduction in engineering/prototype/test costs can be guaranteed. The innovation the developed methodology propose is ensuring a fast and reliable performance assessment via taking the holistic effect of integrated models into consideration. Hence it is possible for original equipment manufacturer to completely or gradually eliminate the actual tests and hence prototype costs, test system, and engineer allocations.

Newly Released Capabilities in the Distributed-Memory SuperLU Sparse Direct Solver

We present the new features available in the recent release of SuperLU_DIST , Version 8.1.1. SuperLU_DIST is a distributed-memory parallel sparse direct solver. The new features include (1) a 3D communication-avoiding algorithm framework that trades off inter-process communication for selective memory duplication, (2) multi-GPU support for both NVIDIA GPUs and AMD GPUs, and (3) mixed-precision routines that perform single-precision LU factorization and double-precision iterative refinement. Apart from the algorithm improvements, we also modernized the software build system to use CMake and Spack package installation tools to simplify the installation procedure. Throughout the article, we describe in detail the pertinent performance-sensitive parameters associated with each new algorithmic feature, show how they are exposed to the users, and give general guidance of how to set these parameters. We illustrate that the solver’s performance both in time and memory can be greatly improved after systematic tuning of the parameters, depending on the input sparse matrix and underlying hardware.

A sequential direct hybrid algorithm to compute stationary distribution of continuous-time Markov chain

In modeling a wide variety of problems of the real world using Markov processes, we usually get large sparse Markov chains. In the case of non-Markovian process, the behavior of a stochastic phenomenon or system can be approximated by continuous-time Markov chains (CTMC) using Phase-type distributions. However, this substantially increases the number of states, which requires a large amount of computation time and operational memory resources. Even the most recent computer hardware is not enough to provide a solution in almost real time, unless a specific algorithm is designed.In this paper, a new hybrid algorithm for the computation of a stationary distribution of a large ergodic homogeneous Markov chains with a finite state space and continuous time is suggested. This algorithm would solve the aforementioned challenges. Depending on model size and sparsity, it significantly reduces a computational time when compared with other recent methods, such as the SparseLU solver from Eigen library (version 3.9.9) and the MUMPS solver (MUltifrontal Massively Parallel sparse direct Solver, version 5.4.0). It is important to note that the algorithm works well independently on the structure of the generator matrix of a continuous-time Markov chain. For the study of the developed hybrid algorithm, the generator matrix is constructed randomly using a special algorithm, based on combining two methods: the slightly modified Grassmann, Taskar and Heyman (GTH) algorithm and the Gauss elimination method. The new idea is to remove the majority of states in a near-optimal order by the GTH method and complete computations for the remaining states by the Gauss method. The robust approach to identify a switching position from one algorithm to another one has been proposed and tested.

A Segregated Approach for Modeling the Electrochemistry in the 3-D Microstructure of Li-Ion Batteries and Its Acceleration Using Block Preconditioners

Battery performance is strongly correlated with electrode microstructure. Electrode materials for lithium-ion batteries have complex microstructure geometries that require millions of degrees of freedom to solve the electrochemical system at the microstructure scale. A fast-iterative solver with an appropriate preconditioner is then required to simulate large representative volume in a reasonable time. In this work, a finite element electrochemical model is developed to resolve the concentration and potential within the electrode active materials and the electrolyte domains at the microstructure scale, with an emphasis on numerical stability and scaling performances. The block Gauss-Seidel (BGS) numerical method is implemented because the system of equations within the electrodes is coupled only through the nonlinear Butler–Volmer equation, which governs the electrochemical reaction at the interface between the domains. The best solution strategy found in this work consists of splitting the system into two blocks—one for the concentration and one for the potential field—and then performing block generalized minimal residual preconditioned with algebraic multigrid, using the FEniCS and the Portable, Extensible Toolkit for Scientific Computation libraries. Significant improvements in terms of time to solution (six times faster) and memory usage (halving) are achieved compared with the MUltifrontal Massively Parallel sparse direct Solver. Additionally, BGS experiences decent strong parallel scaling within the electrode domains. Last, the system of equations is modified to specifically address numerical instability induced by electrolyte depletion, which is particularly valuable for simulating fast-charge scenarios relevant for automotive application.

Parallel solution of saddle point systems with nested iterative solvers based on the Golub‐Kahan Bidiagonalization

SummaryThe Golub‐Kahan bidiagonalization is widely used in the singular value decomposition of rectangular matrices and has been generalized to an iterative solver for symmetric indefinite linear systems with a two‐by‐two block structure. In this work, we present a scalability study of this generalized solver as implemented in a recent release of the parallel numerical library PETSc (Portable, Extensible Toolkit for Scientific Computation). We present an improved solver performance for the two‐dimensional (2D) Stokes equations as compared to previous work. Furthermore, we investigate the performance of different parallel inner solvers in the outer Golub‐Kahan iteration for a three‐dimensional Stokes problem. The study includes parallel sparse direct solvers and multigrid methods. When increasing the number of cores for a fixed total problem size, the solver exhibits good speedups of up to 50% at the 1024 core count. For the tests in which the total problem size grows while the workload in each core stays constant, the parallel performance of the solver scales almost linearly with the increase in the core counts. In particular, the computation time increases only by about 15% when the number of cores increases from 80 to 1024 for a 2D test case.

Auto-Gauging of Vector Potential by Parallel Sparse Direct Solvers—Numerical Observations

When using magnetic vector potential (MVP)-based formulations for magnetostatic or eddy-current problems, either gauge conditions specifying the divergence of the MVP or tree gauging by eliminating redundant degrees of freedom of the MVP is usually imposed to ensure uniqueness of solutions. Explicit gauging of the MVP is not always necessary since classical iterative solvers can automatically and implicitly fix the gauge as long as the right-hand side vectors are consistent. Besides iterative solvers, implicit gauging is also observed when using state-of-the-art parallel sparse direct solvers (PSDSs), thanks to the built-in functions of handling null-spaces of either real symmetric positive semi-definite matrix systems or those complex symmetric systems from eddy-current problems. Both static and eddy-current examples are solved by PSDS to demonstrate results of local physical quantities or global quantities such as magnetic energy or joule losses. High-order edge/nodal elements are also considered in our numerical examples and it is observed that PSDS can also easily and correctly handle the delicate discrete null spaces.

Improved Equilibrated Error Estimates for Open Boundary Magnetostatic Problems Based on Dual A and H Formulations

Calculating the bounds of global energy is an important issue in computational electromagnetism, which can provide guaranteed results when extracting inductance parameters. In this paper, an improved equilibrated type a posteriori error estimate for open boundary magnetostatic problems is proposed. We derive our error estimator based on vector dual formulations, which can be efficiently solved using parallel sparse direct solvers. The new estimator can provide a sharp and guaranteed estimate of the finite-element spatial discretization error. Moreover, the computational cost is cheaper than using existing equilibrated error estimators. Numerical experiments are carried out to showcase the performance of our error estimator, including the modified TEAM workshop problem 13 and the benchmark IEEJ problem.

Modeling kV X-ray-Induced Coloration in Radiochromic Films

Commercially available radiochromic films are primarily designed for clinical X-ray dosimetry. These films change color upon exposures to radiation as a result of solid-state polymerization (SSP). Built on a previous model developed for SSP upon exposures to ultraviolet (UV) radiation, a new model was developed in the present work for X-ray-induced coloration in Gafchromic EBT3 films. Monte Carlo simulations using the Monte Carlo N-Particle (MCNP) code were employed to model the transport and interaction of photons and the generated secondary electrons within the film active layer. The films were exposed to continuous-energy photon beams. The dose DE in the external radiation detector (i.e., ionization chamber) was determined and the realistic dose DA in the film active layer was then obtained using the calibration coefficient R (=DA/DE). The finite element method (FEM) was used to solve the classical steady-state Helmholtz equation using the multifrontal massively parallel sparse direct solver (MUMPS). An extensive grid independence test was carried out and the numerical stability of the present model was ensured. The reflected light intensity from the film surface was used to theoretically obtain the net reflective optical density of the film exposed to X-ray. Good agreement was obtained between the experimental and theoretical results of the net reflective optical density of the film. For X-ray doses >~600 cGy, due to the already formed densely cross-linked structure in the active layer of the EBT3 film so further bond formation was less likely, the reflected light intensity from the film surface increased at a relatively lower rate when compared to those for dose values <~600 cGy.

Selected inversion as key to a stable Langevin evolution across the QCD phase boundary

We present new results of full QCD at nonzero chemical potential. In PRD 92, 094516 (2015) the complex Langevin method was shown to break down when the inverse coupling decreases and enters the transition region from the deconfined to the confined phase. We found that the stochastic technique used to estimate the drift term can be very unstable for indefinite matrices. This may be avoided by using the full inverse of the Dirac operator, which is, however, too costly for four-dimensional lattices. The major breakthrough in this work was achieved by realizing that the inverse elements necessary for the drift term can be computed efficiently using the selected inversion technique provided by the parallel sparse direct solver package PARDISO. In our new study we show that no breakdown of the complex Langevin method is encountered and that simulations can be performed across the phase boundary.

2.5D marine CSEM modeling in the frequency-domain based on an improved interpolation scheme at receiver positions

In the marine controlled-source electromagnetic (CSEM) survey, the receivers are usually placed at the seafloor. The resistivity contrast between the seawater and seafloor sediments is large, which can cause difficulties in numerical modeling of CSEM fields at receiver locations. In this paper, we present an improved interpolating method for calculating electric and magnetic fields at the seafloor with a resistivity contrast. This method is applied to the 2. 5 dimensional (2. 5D) frequency-domain CSEM modeling with towed transmitters and receivers located at the seafloor. Considering the discontinuity of the normal electric fields, we use the normal current electric density for interpolation. We simulate the 2. 5D marine CSEM responses by the staggered finite-difference (SFD) method with Fourier transform to the strike direction. The final SFD equations are solved by the direct solver MUMPS (MUltifrontal Massively Parallel Sparse direct Solver). To avoid the source singularities, the secondary-field approach is used and the primary fields excited by the electric dipole source can be calculated quasi-analytically for the one-dimensional (1D) layered background model. We focus on interpolating of electric and magnetic fields in the wavenumber domain to the receiver locations at the seafloor interface between the conductive seawater and resistive seafloor formation. The secondary electric and magnetic fields are used for interpolation instead of the total fields for high numerical accuracy. After performing the inverse Fourier transform to the wavenumbers, the electric and magnetic fields in the space domain are obtained. To check the accuracy of our 2. 5D marine CSEM SFD modeling algorithm with the improved receiver interpolating technique, we compare our results with both the 1D analytical results and the adaptive finite element results. The SFD numerical results are approved to be accurate. We also compare the numerical accuracy between our improved interpolation scheme and others, i.e., the conventional linear interpolation and the rigorous interpolation. The proposed interpolation only utilizes the nodes below/above the seafloor interface, and is proved to be much more accurate than the other two interpolating methods used.

Modeling Coloration of a Radiochromic Film with Molecular Dynamics-Coupled Finite Element Method

Radiochromic films change color upon exposures to radiation doses as a result of solid-state polymerization (SSP). Commercially available radiochromic films are primarily designed for, and have become widely used in, clinical X-ray dosimetry. However, many intriguing properties of radiochromic films are not yet fully understood. The present work aimed at developing a theoretical model at both atomic and macroscopic scales to provide a platform for future works to understand these intriguing properties. Despite the fact that radiochromic films were primarily designed for clinical X-ray dosimetry, dose-response curves for the Gafchromic EBT3 film obtained for ultraviolet (UV) radiation were employed to develop our model in order to avoid complications of ionization, non-uniform energy deposition, as well as dispersed doses caused by secondary electrons set in motion by the indirectly ionizing X-ray photons, which might introduce added uncertainties to the model and overshadow the basic SSP processes. The active layer in the EBT3 film consisted of diacetylene (DA) pentacosa-10,12-diynoate monomers, which were modelled using molecular dynamics (MD). The degrees of SSP in the atomic scale upon different UV exposures were obtained to determine the absorption coefficients of the active layer, which were then input into the finite element method (FEM). The classical steady-state Helmholtz equation was engaged to model the reflection from the active layer using the FEM technique. The multifrontal massively parallel sparse direct solver (MUMPS) was employed to solve the present numerical problem. Very good agreement between experimentally and theoretically obtained coloration in terms of net reflective optical density was achieved for different UV exposures. In particular, for UV exposures larger than ~40 J/cm2, the reflected light intensity decreased at a lower rate when compared to other UV exposure values, which was explained by the densely cross-linked structure under near-complete polymerization, and thus the lower efficiency for further bond formation between DA monomer strands.

Dimensionless modeling for convective drying of tuberous crop (Solanum tuberosum) by considering shrinkage

The shrinkage of biological materials (agri-food), is a common phenomenon developed during a drying operation that directly affects the product quality. In this work a dimensionless mathematical model including the reduction at the interface of a high moisture content food (tuber) in 1D geometry was developed. The mass and heat conservation equations coupled to the Leibniz-Reynolds Transport Theorem (LRTT) at the interface of a high moisture content food were solved. The system of nonlinear unsteady state equations was solved by using a multi-frontal massively parallel sparse direct Solver. In order to validate this model, tuberous crop slices were dried at three temperatures: 40, 50 and 60 °C, relative humidity of 25% and a constant airflow of 1.5 m/s. Drying kinetics and the evolution of temperature within the product were logged. The model is able to simulate the moisture loss and predict the highest thickness change. The simulations show a good agreement with the experimental drying kinetics, as well as for the evolution of temperature. Simulations show 1-D shrinkage is not linear and the model improves its prediction when the drying temperature is higher (intensive mass flux).

On Using Cholesky-Based Factorizations and Regularization for Solving Rank-Deficient Sparse Linear Least-Squares Problems

By examining the performance of modern parallel sparse direct solvers and exploiting our knowledge of the algorithms behind them, we perform numerical experiments to study how they can be used to efficiently solve rank-deficient sparse linear least-squares problems arising from practical applications. The Cholesky factorization of the normal equations breaks down when the least-squares problem is rank-deficient, while applying a symmetric indefinite solver to the augmented system can give an unacceptable level of fill in the factors. To try to resolve these difficulties, we consider a regularization procedure that modifies the diagonal of the unregularized matrix. This leads to matrices that are easier to factorize. We consider both the regularized normal equations and the regularized augmented system. We employ the computed factors of the regularized systems as preconditioners with an iterative solver to obtain the solution of the original (unregularized) problem. Furthermore, we look at using limited-me...

PSPIKE+: A family of parallel hybrid sparse linear system solvers

We present PSPIKE+, as a family of parallel hybrid linear system solvers that are more robust than other available preconditioned iterative methods, and more scalable than parallel sparse direct solvers. A critical step in using PSPIKE+ is the extraction of a preconditioner that is in the form of a generalized central band which consists of a set of overlapping diagonal blocks that encapsulate as many of the heaviest elements of the coefficient matrix. To address this issue, we propose the G-PAVER reordering scheme that produces an effective preconditioner in the form of overlapping diagonal blocks while simultaneously accommodating load balancing on distributed-memory architectures and controlling overlap sizes. In each outer Krylov subspace iteration for solving the linear system, we need to solve linear systems involving the preconditioner. This is accomplished using a parallel tearing-based method, which may be regarded as an algebraic domain decomposition scheme. The tearing-based method requires the solution of a much smaller linear system (we call the balance system) whose size is equal to the sum of the overlaps. In this paper, we form the balance system explicitly and solve it directly. Our parallel implementation of PSPIKE+ based on G-PAVER, PARDISO 5.0.0, and solving the explicit balance system results in impressive robustness and good parallel scalability on computing platforms with a limited number of multi-core nodes.

Laplace-domain waveform modeling and inversion for the 3D acoustic–elastic coupled media

Laplace-domain waveform inversion reconstructs long-wavelength subsurface models by using the zero-frequency component of damped seismic signals. Despite the computational advantages of Laplace-domain waveform inversion over conventional frequency-domain waveform inversion, an acoustic assumption and an iterative matrix solver have been used to invert 3D marine datasets to mitigate the intensive computing cost. In this study, we develop a Laplace-domain waveform modeling and inversion algorithm for 3D acoustic–elastic coupled media by using a parallel sparse direct solver library (MUltifrontal Massively Parallel Solver, MUMPS). We precisely simulate a real marine environment by coupling the 3D acoustic and elastic wave equations with the proper boundary condition at the fluid-solid interface. In addition, we can extract the elastic properties of the Earth below the sea bottom from the recorded acoustic pressure datasets. As a matrix solver, the parallel sparse direct solver is used to factorize the non-symmetric impedance matrix in a distributed memory architecture and rapidly solve the wave field for a number of shots by using the lower and upper matrix factors. Using both synthetic datasets and real datasets obtained by a 3D wide azimuth survey, the long-wavelength component of the P-wave and S-wave velocity models is reconstructed and the proposed modeling and inversion algorithm are verified. A cluster of 80CPU cores is used for this study.

Frequency-domain elastic full-waveform multiscale inversion method based on dual-level parallelism

The complexity of an elastic wavefield increases the nonlinearity of inversion. To some extent, multiscale inversion decreases the nonlinearity of inversion and prevents it from falling into local extremes. A multiscale strategy based on the simultaneous use of frequency groups and layer stripping method based on damped wave field improves the stability of inversion. A dual-level parallel algorithm is then used to decrease the computational cost and improve practicability. The seismic wave modeling of a single frequency and inversion in a frequency group are computed in parallel by multiple nodes based on multifrontal massively parallel sparse direct solver and MPI. Numerical tests using an overthrust model show that the proposed inversion algorithm can effectively improve the stability and accuracy of inversion by selecting the appropriate inversion frequency and damping factor in lowfrequency seismic data.

Parallel direct solver for solving systems of linear equations resulting from finite element method on multi-core desktops and workstations

This paper presents the PARFES—parallel sparse direct solver for solving systems of linear equations resulting from the finite element method applied to problems of structural and solid mechanics. PARFES is intended for use on shared-memory multi-core desktops, laptops and workstations. The solver’s high performance relies on the use of a specific algorithm for the procedure of block matrix multiplication using AVX instructions for the vectorization of computations, and the parallelizing approach that ensures a balanced processor load. The block looking-left factorization algorithm allows for the use of disk memory if the dimensionality of problem exceeds the capacity of RAM.

Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation

Computationally efficient 3-D frequency-domain full waveform inversion (FWI) is applied to ocean-bottom cable data from the Valhall oil field in the visco-acoustic vertical transverse isotropic (VTI) approximation. Frequency-domain seismic modelling is performed with a parallel sparse direct solver on a limited number of computer nodes. A multiscale imaging is performed by successive inversions of single frequencies in the 3.5–10 Hz frequency band. The vertical wave speed is updated during FWI while density, quality factor QP and anisotropic Thomsen's parameters δ and ϵ are kept fixed to their initial values. The final FWI model shows the resolution improvement that was achieved compared to the initial model that was built by reflection traveltime tomography. This FWI model shows a glacial channel system at 175 m depth, the footprint of drifting icebergs on the palaeo-seafloor at 500 m depth, a detailed view of a gas cloud at 1 km depth and the base cretaceous reflector at 3.5 km depth. The relevance of the FWI model is assessed by frequency-domain and time-domain seismic modelling and source wavelet estimation. The agreement between the modelled and recorded data in the frequency domain is excellent up to 10 Hz although amplitudes of modelled wavefields propagating across the gas cloud are overestimated. This might highlight the footprint of attenuation, whose absorption effects are underestimated by the homogeneous background QP model (QP = 200). The match between recorded and modelled time-domain seismograms suggests that the inversion was not significantly hampered by cycle skipping. However, late arrivals in the synthetic seismograms, computed without attenuation and with a source wavelet estimated from short-offset early arrivals, arrive 40 ms earlier than the recorded seismograms. This might result from dispersion effects related to attenuation. The repeatability of the source wavelets inferred from data that are weighted by a linear gain with offset is dramatically improved when they are estimated in the FWI model rather than in the smooth initial model. The two source wavelets, estimated in the FWI model from data with and without offset gain, show a 40 ms time-shift, which is consistent with the previous analysis of the time-domain seismograms. The computational efficiency of our frequency-domain approach is assessed against a recent time-domain FWI case study performed in a similar geological environment. This analysis highlights the efficiency of the frequency-domain approach to process a large number of sources and receivers with limited computational resources, thanks to the efficiency of the substitution step performed by the direct solver. This efficiency can be further improved by using a block-low rank version of the multifrontal solver and by exploiting the sparsity of the source vectors during the substitution step. Future work will aim to update attenuation and density at the same time of the vertical wave speed.