Algebraic Multigrid Solver Research Articles

A Fast Algebraic Multigrid Solver and Accurate Discretization for Highly Anisotropic Heat Flux I: Open Field Lines

A Fast Algebraic Multigrid Solver and Accurate Discretization for Highly Anisotropic Heat Flux I: Open Field Lines

DC3DPAFEM: An efficient and accurate 3-D direct current resistivity anisotropic forward modeling software for complex geological settings

Nowadays, there is a growing trend that direct current (DC) field surveys are shifting towards challenging areas characterized by mountainous topography and electrical anisotropy. Given these complex geological settings, there is an urgent need for 3-D DC forward modeling software capable of effectively addressing large-scale problems and delivering accurate modeling results to interpret field data. However, most open-source software packages face certain limitations, such as the high numerical cost to handle complex surface topography, the lack of consideration for anisotropic conductivity, the absence of mesh refinement techniques to guarantee accuracy in forward modeling, and the lack of parallel computing techniques to solve large-scale problems. In this study, we develop an efficient and highly accurate 3-D DC anisotropic forward modeling software, namely DC3DPAFEM, using the adaptive finite element algorithm based on the unstructured tetrahedral mesh. Firstly, we construct a strong compatible boundary value problem (BVP) for 3-D anisotropic DC problems by adopting a specialized secondary potential approach to handle the surface topography efficiently. Then, we develop a goal-oriented adaptive mesh refinement (AMR) technique to ensure accurate forward modeling results, even with a coarse initial mesh. To ensure time and memory efficiency, we employ a robust conjugate gradient (CG) algorithm preconditioned by the algebraic multigrid (AMG) solver to solve the large-scale linear system of equations resulting from complex geological structures. We aim to investigate the performance of the AMG scheme in anisotropic DC cases. Furthermore, we incorporate the domain decomposition technique into the iterative solution scheme for further efficiency gains. This technique significantly improves computing efficiency for large-scale problems in parallel clusters. Finally, we conduct comprehensive performance tests for DC3DPAFEM using a two-layer anisotropic model and a 3-D complex model with undulating terrain. The results of both examples validate the accuracy of DC3DPAFEM, as they closely align with the analytical solutions and the solutions obtained from the existing 3-D DC forward modeling code. Compared to traditional direct solver MUMPS and ILU-preconditioned iterative solvers, DC3DPAFEM exhibits highly scalable performance for large-scale problems, offering significant advantages in terms of memory and time consumption. Overall, DC3DPAFEM demonstrates substantial advances in efficiency, accuracy, and practicality through a series of numerical examples. This open-source code provides an efficient and available tool for developing a 3-D DC inversion method that can deal with large-scale problems involving intricate topography and anisotropic media.

A multiscale preconditioner for crack evolution in porous microstructures: Accelerating phase‐field methods

AbstractPhase‐field methods are attractive for simulating the mechanical failure of geometrically complex porous microstructures described by 2D/3D x‐ray CT images in subsurface (e.g., CO storage) and manufacturing (e.g., Li‐ion battery) applications. They capture the nucleation, growth, and branching of fractures without prior knowledge of the propagation path or having to remesh the domain. Their drawback lies in the high computational cost for the typical domain sizes encountered in practice. We present a multiscale preconditioner that significantly accelerates the convergence of Krylov solvers in computing solutions of linear(ized) systems arising from the sequential discretization of the momentum and crack‐evolution equations in phase‐field methods. The preconditioner is an algebraic reformulation of a recent pore‐level multiscale method (PLMM) by the authors and consists of a global preconditioner and a local smoother . Together, and attenuate low‐ and high‐frequency errors simultaneously. The proposed , used in the momentum equation only, is a simplification of a recent variant proposed by the authors that is much cheaper and easier to deploy in existing solvers. The smoother , used in both the momentum and crack‐evolution equations, is built such that it is compatible with and more robust and efficient than black‐box smoothers like ILU(). We test and systematically for static‐ and evolving‐crack problems on complex 2D/3D porous microstructures, and show that they outperform existing algebraic multigrid solvers. We also probe different strategies for updating as cracks evolve and show the associated cost can be minimized if is updated adaptively and infrequently. Both and are scalable on parallel machines and can be implemented non‐intrusively in existing codes.

Linear-Time Solution of 3D Magnetostatics With Aggregation-Based Algebraic Multigrid Solvers

Linear-Time Solution of 3D Magnetostatics With Aggregation-Based Algebraic Multigrid Solvers

Light Weight Coarse Grid Aggregation for Smoothed Aggregation Algebraic Multigrid Solver

Light Weight Coarse Grid Aggregation for Smoothed Aggregation Algebraic Multigrid Solver

Bringing physics into the coarse‐grid selection: Approximate diffusion distance/effective resistance measures for network analysis and algebraic multigrid for graph Laplacians and systems of elliptic partial differential equations

AbstractIn a recent paper, the author examined a correlation affinity measure for selecting the coarse degrees of freedom (CDOFs) or coarse nodes (C nodes) in systems of elliptic partial differential equations (PDEs). This measure was applied to a set of relaxed vectors, which exposed the near‐nullspace components of the PDE operator. Selecting the CDOFs using this affinity measure and constructing the interpolation operators using a least‐squares procedure, an algebraic multigrid (AMG) method was developed. However, there are several noted issues with this AMG solver. First, to capture strong anisotropies, a large number of test vectors may be needed; and second, the solver's performance can be sensitive to the initial set of random test vectors. Both issues reflect the sensitive statistical nature of the measure. In this article, we derive several other statistical measures that ameliorate these issues and lead to better AMG performance. These measures are related to a Markov process, which the PDE itself may model. Specifically, the measures are based on the diffusion distance/effective resistance for such process, and hence, these measures incorporate physics into the CDOF selection. Moreover, because the diffusion distance/effective resistance can be used to analyze graph networks, these measures also provide a very economical scheme for analyzing large‐scale networks. In this article, the derivations of these measures are given, and numerical experiments for analyzing networks and for AMG performance on weighted‐graph Laplacians and systems of elliptic boundary‐value problems are presented.

Fast large-scale gravity modeling using a high-order compact cascadic multigrid strategy

A new high-order compact (HOC) cascadic multigrid (MG) strategy is developed for modeling large-scale complex gravity problems on nonuniform rectilinear grids. First, the governing 3D Poisson’s gravitational potential equation is discretized by a novel HOC difference scheme, which allows for Robin boundary conditions to further reduce the computational domain. Then, the resulting asymmetric linear system of equations is solved by a generalized extrapolation cascadic MG (gEXCMG) method augmented with a novel MG prolongation operator. The quartic interpolation and completed Richardson extrapolation assist this MG prolongation operator to produce a good initial guess for the biconjugate gradient stabilized smoother. Two synthetic and SEG/EAGE salt models are examined to verify the presented approach. Results indicate that our algorithm can efficiently and accurately produce the gravity signals for complex geologic models. Moreover, the comparisons with the state-of-the-art algebraic MG solvers demonstrate that our gEXCMG solver offers substantially better efficiency with significantly less memory consumption. Therefore, our newly developed method has the potential for serving as a forward engine in large-scale gravity inversions.

Three-precision algebraic multigrid on GPUs

Recent research has demonstrated that using low precision inside some levels of an algebraic multigrid (AMG) solver can improve performance without negatively impacting the AMG quality. In this paper, we build upon previous research and implement an AMG that can use double, single, and half precision for the distinct multigrid levels. The implementation is platform-portable across GPU architectures from AMD, Intel, and NVIDIA. In an experimental analysis, we demonstrate that the use of half precision can be a viable option in multigrid. We evaluate the performance of different AMG configurations and demonstrate that mixed precision AMG can provide runtime savings compared to a double precision AMG.

Algebraic Multigrid Using a Stencil–CSR Hybrid Format on GPUs

We propose a new sparse matrix format which captures the matrix structure typical for discretized partial differential equations with piecewise-constant coefficients. The format uses a stencil representation for some blocks of matrix rows, typically corresponding to regions with constant coefficients, whereas other rows are encoded in the general compressed sparse row format. The stencil representation saves memory and is suitable for SIMD-like parallelism as available on GPUs. Further, this format is well suited for the implementation of algebraic multigrid methods, and we present a proof-of-concept GPU-accelerated aggregation-based algebraic multigrid solver based on this format. This solver is compared on a few model problems (2-dimension and 3-dimension Poisson-like) with the compressed sparse row–based solver AmgX from NVIDIA and with the CUDA version of the BoomerAMG solver. For the considered tested problems with one million unknowns or more, the presented solver outperforms AmgX and BoomerAMG in terms of both run time and memory usage, and the performance gap increases with the system size.

Controlled-source electromagnetic modeling using a versatile secondary electric-field formulation and efficient multigrid-based preconditioner

Frequency-domain controlled-source electromagnetic (CSEM) surveying is a widely used geophysical electromagnetic exploration method. To avoid the singularity near the transmitting sources and improve the accuracy of three-dimensional (3D) CSEM forward modeling, we adopt the vector finite element method and secondary electric field formulation combined with hexahedral grids to solve the frequency-domain Maxwell equations. Primary field computation applies the fast Hankel transform using digital linear filters to quickly calculate the electromagnetic field generated by different electromagnetic dipoles or loops in the layered medium. Therefore, our forward modeling algorithm is highly versatile and compatible with various electromagnetic sources and scenarios of CSEM application (e.g., airborne, offshore, or marine). To improve the convergence rate, we use an algebraic multigrid solver based on the established Hiptmair–Xu decomposition. Using direct solvers requires a vast amount of memory, while conventional iterative solvers usually fail or converge only very slowly. Our new iterative scheme does not require the divergence corrections commonly applied in finite difference or finite volume electromagnetic modeling to accelerate the convergence. We implement complete parallelism in our finite element programming to efficiently solve the CSEM forward problems on massively parallel clusters. We use different numerical methods to obtain high-precision comparison results for one-dimensional and typical 3D models to demonstrate our algorithm's effectiveness and correctness. Finally, a public model of 3 million cells is used to simulate the responses of marine CSEM surveying. The proposed new solver shows robust scalability when using nearly 1000 cores on the cluster, demonstrating that our algorithm is particularly suitable for large-scale 3D forward modeling.

Stationary Stokes solver for single-phase flow in porous media: A blastingly fast solution based on Algebraic Multigrid Method using GPU

The paper is focused on high efficiency Stokes solver that is applied to the incompressible flow in porous media. Computational domains are represented by binarized 3D computed tomography voxel models. The solution procedure is constructed around a fast algebraic multigrid solver that utilizes the power of graphics processing units (GPUs). In order to minimize memory footprint and accelerate the solver, a simple MAC-type staggered finite difference discretization is used and a coupled Stokes saddle-point type system is solved directly on a GPU. The MAC discretization is discussed, taking particular care of internal solid boundaries around pores. The method includes topological domain analysis on a GPU, which removes isolated volumes with no flow in the domain (based on the connected component labeling), depending on the boundary conditions. We consider various types of boundary conditions and efficient parallel strategies for the GPUs, including fast matrix assembly and residual regularization. Our method is extensively benchmarked both against analytical solutions and applied problems from digital rock physics in terms of computational wall time, precision, approximation order, and convergence. We demonstrate that it takes up to 5–23 s on a modern GPU card to obtain a solution with 1⋅10−6 error residuals for 3D geometries with 300–4503 voxels and a porosity range of 5–37%.

An efficient extrapolation multigrid method based on a HOC scheme on nonuniform rectilinear grids for solving 3D anisotropic convection–diffusion problems

We develop an efficient multigrid method combined with a high-order compact (HOC) finite difference scheme on nonuniform rectilinear grids for solving 3D diagonal anisotropic convection–diffusion problems with boundary/interior layers. Firstly, we derive a fourth-order compact finite difference scheme to discretize the 3D anisotropic convection–diffusion equation on a rectilinear grid. Then, the resulting large-scale asymmetric linear system of equations is solved by a generalized extrapolation cascadic multigrid (gEXCMG) method based on two novel multigrid (MG) prolongation operators. The highlight of this paper is the application of the quintic Lagrange interpolation and the completed Richardson extrapolation in the design of the new MG prolongation operator on nonuniform rectilinear grids, which can produce a good initial guess (sixth-order approximation to the finite difference solution) for the SSOR-preconditioned biconjugate gradient stabilized (BiCGStab) smoother. In the end, numerical experiments show that the gEXCMG method combined with the HOC scheme can achieve fourth-order accuracy for 3D anisotropic convection–diffusion problems with few smoothing steps on the finest grid. Moreover, the proposed gEXCMG method can offer substantially better efficiency than the state-of-the-art algebraic MG solver, namely, aggregation-based algebraic multigrid (AGMG) method, for large linear systems arising from the discretization of second order elliptic PDEs.

Thermodynamic characterization of various funnels subjected to conjugate natural convection

A thermal system becomes energy competent only when the entropy generation becomes as less as possible. Therefore, the thermodynamic characterization of different funnels used as IRS devices is investigated. Also, the outcome of the heat transfer study is also presented. The main aim is to select the best funnel for maximum heat transfer rate and minimum entropy production. Different equations, namely Navier-Stokes equation, turbulence equation, energy equations, and entropy production equation are solved using finite volume method based algebraic multigrid solver of ANSYS-Fluent. Surface temperature, Rayleigh number, and different shapes of funnels are used as parameters to conduct the study. The results reveal that the cylindrical louvered funnel outperforms the rest of the considered funnel in terms of heat transfer rate, increasing with surface temperature and Rayleigh number. Irreversibility caused by fluid friction is found to be negligible as compared to heat transfer. A performance evaluation criterion (ratio of irreversibility to heat transfer) measures the performance of different funnels. Results indicated that it decreases with Ra, indicating that the rise in heat transfer becomes much more than the rise in irreversibility at high Ra. The contours of velocity, temperature, and entropy generation are also provided for better understanding. Also, the time-temperature graph shows that the louvered cylindrical funnel cools fastest while the louvered conical funnel is the slowest.

A Unified Arbitrary Lagrangian\u2013Eulerian Model for Fluid\u2013Structure Interaction Problems Involving Flows in Flexible Channels

In this work a finite element-based model for analyzing incompressible flows in flexible channels is presented. The model treats the fluid–solid interaction problem in a monolithic way, where the governing equations for both sub-domains are solved on a single moving grid taking advantage of an arbitrary Lagrangian/Eulerian framework (ALE). The unified implementation of the governing equations for both sub-domains is developed, where these are distinguished only in terms of the mesh-moving strategy and the constitutive equation coefficients. The unified formulation is derived considering a Newtonian incompressible fluid and a hypoelastic solid. Hypoelastic constitutive law is based on the strain rate and thus naturally facilitates employing velocity as a kinematic variable in the solid. Unifying the form of the governing equations and defining a semi-Lagrangian interface mesh-motion algorithm, one obtains the coupled problem formulated in terms of a unique kinematic variable. Resulting monolithic system is characterized by reduced variable heterogeneity resembling that of a single-media problem. The model used in conjunction with algebraic multigrid linear solver exhibits attractive convergence rates. The model is tested using a 2D and a 3D example.

Porting an aggregation-based algebraic multigrid method to GPUs

We present a hybrid GPU-CPU version of the AGMG software, a popular algebraic multigrid (AMG) solver which implements an aggregation-based AMG method. With the new implementation, the solution stage runs on a GPU, except operations on the coarsest grid, which are executed on a CPU. To maximize the speedup, two novel %new features are introduced. On the one hand, ℓ1-Jacobi smoothing is combined with polynomial acceleration (or polynomial smoothing), leading to improved performance compared with standard ℓ1-Jacobi smoothing, while not requiring to compute eigenvalue estimates as standard polynomial smoothing does. On the other hand, besides the K-cycle used in standard AGMG, we introduce the relaxed W-cycle, which tends to combine the advantages of the K-cycle and the standard W-cycle. Numerical results show that the new implementation inherits the robustness of the original AGMG software, while bringing significant speedups on GPUs. A comparison with AmgX, a reference AMG solver from NVIDIA, suggests that the presented hybrid GPU-CPU version of AGMG is more robust and often significantly faster in the solution stage.

Fast Solution of 3-D Eddy-Current Problems in Multiply Connected Domains by a, v-φ and t-φ Formulations With Multigrid-Based Algorithm for Cohomology Generation

The fast solution of three-dimensional eddy current problems is still an open problem, especially when real-size finite element models with millions of degrees of freedom are considered. In order to lower the number of degrees of freedom a magnetic scalar potential can be used in the insulating parts of the model. This may become difficult when the model geometry presents some conductive parts which are multiply connected. In this work a multigrid-based algoritm is proposed that allows for a calculation in linear-time of cohomology, which is needed to introduce the scalar potential without cuts. This algorithm relies on an algebraic multigrid solver for curl-curl field problems, which ensures optimal computational complexity. Numerical results show that the novel algorithm outperforms state-of-the-art methods for cohomology generation based on homological algebra. In addition, based on this algoritm, novel <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$a,v$ </tex-math></inline-formula> - <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$\varphi $ </tex-math></inline-formula> and <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$t$ </tex-math></inline-formula> - <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$\varphi $ </tex-math></inline-formula> formulations to analyze three-dimensional eddy current problems in multiply connected domains are proposed. Both formulations, after discretization by the cell method, lead to a complex symmetric system of linear equations amenable to fast iterative solution by Krylov-subspace solvers. These formulations are able to provide very accurate numerical results, with a minimum amount of degrees of freedoms to represent the eddy current model. In this way the computational performance is improved compared to the classical <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$A,V$ </tex-math></inline-formula> - <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$A$ </tex-math></inline-formula> formulation typically implemented in finite element software for electromagnetic design.

A mass conserving arbitrary Lagrangian–Eulerian formulation for three‐dimensional multiphase fluid flows

AbstractThe arbitrary Lagrangian–Eulerian (ALE) framework presented in [Sahin and Guventurk, An arbitrary Lagrangian–Eulerian framework with exact mass conservation for the numerical simulation of 2D rising bubble problem.Int J Numer Methods Eng. 2017;112:2110–2134] has been extended to simulate three‐dimensional immiscible incompressible multiphase fluid flows. The div‐stable side centered unstructured finite volume formulation is used for the discretization of the incompressible Navier–Stokes equations. A novel kinematic boundary condition that satisfies the local discrete geometric conservation law (DGCL) is implemented to ensure the mass conservation of each species at the machine precision. The pressure field is treated to be discontinuous across the interface with the discontinuous treatment of density and viscosity. Surface tension force is treated as a tangent force and discretized in a semi‐implicit form. Two different approaches for the computation of unit normal vector have been implemented: the least squares biquadratic surface fitting (LSBSF) and the mean weighted by sine and edge length reciprocals (MWSELR). The resulting large system of algebraic equations is solved in a fully coupled manner in order to improve the time step restrictions. As a preconditioner, an approximate matrix factorization similar to that of the projection method is employed and the parallel algebraic multigrid solver BoomerAMG provided by the HYPRE library, which is accessed through the PETSc library, has been utilized for the scaled discrete Laplacian of pressure and the diagonal blocks of mesh deformation equations. The accuracy of the proposed method is initially validated on the static bubble problem, since the surface tension force is highly sensitive to the accurate evaluation of the unit normal vector and the inaccuracies significantly contribute to unphysical velocities, called parasitic currents. The calculations indicate that the parasitic currents can be reduced to machine precision for the MWSELR method. Then, the proposed approach is applied to the single bubble rising in a viscous quiescent liquid for both low and high density ratios. The calculations produce accurate predictions of the bubble shape, center of mass, rise velocity, and so forth. Furthermore, the mass of each species is conserved at machine precision and discontinuous pressure field is obtained in order to avoid errors due to the incompressibility restriction in the vicinity of liquid‐liquid interfaces at large density and viscosity ratios.

Numerical investigation of flow and heat transfer characteristics of a full-scale infrared suppression device with cylindrical funnels

Numerical investigation of fluid flow and heat transfer characteristics of a full-scale Infrared suppression (IRS) device for ocean liners is carried out in the present study. The main objective of the work is to reduce the temperature of the turbine exhaust by incorporating the IRS device, which can be effectively used in stealth technology for naval ships. In this study, the IRS device having multiple cylindrical funnels are used for conducting the computational experiment. The numerical solution of the Navier-Stokes equation along with the turbulence as well as the energy equation is solved using the algebraic multigrid solver of ANSYS Fluent 15.0. Many pertinent input parameters, like the inlet Reynolds number, fluid temperature, the number of funnels, and their sizes with overlapping or without it, have been varied in a practical range to study the mass suction into the IRS device. It is observed that the mass entrainment increases with an increase in the nozzle fluid temperature as well as its velocity. Visualization of temperature plume and flow field over the IRS device with multiple funnels has been picturized in this analysis. Furthermore, correlations for the prediction of mass entrainment rate as a function of pertinent input parameters have been developed, which can be used for academic and industrial purposes.

Bucket-based multigrid preconditioner for solving pressure Poisson equation using a particle method

The execution time for incompressible particle methods is dominated by solving large sparse linear systems. In this study, a novel simple bucket-based multigrid (BMG) preconditioner is presented to retrieve linear scaling. In the algorithm, the domain is decomposed into cubic boxes, to enable recursive aggregation of the closest neighboring particles. Under a moving particle semi-implicit (MPS) discretization scheme, parametric, verification, and performance studies were conducted for basic problems. From the parametric study, the BMG preconditioner was accompanied with a Krylov subspace accelerated multigrid cycle strategy and incorporated into a generalized conjugate residual method. Regardless of dimension and degree of particle distribution irregularity, the method significantly outperformed a conjugate gradient (CG), an incomplete Cholesky CG, and a conventional plain aggregation-based algebraic multigrid preconditioned solver. For a dam-break problem with 1.3 M and 2.3 M particles in 2D and 3D, the proposed method, as a linear solver for the MPS method, was 27 and 5.4 times faster than the CG method.

Computational homogenization with million-way parallelism using domain decomposition methods

Parallel computational homogenization using the well-knwon $$\hbox {FE}^2$$ approach is described and combined with domain decomposition and algebraic multigrid solvers. It is the purpose of this paper to show that and how the $$\hbox {FE}^2$$ method can take advantage of the largest supercomputers available and those of the upcoming exascale era for virtual material testing of micro-heterogeneous materials such as advanced steel. The $$\hbox {FE}^2$$ method is a computational micro-macro homogenization approach where at each Gauss integration point of the macroscopic finite element problem a microscopic finite element problem, defined on a representative volume element (RVE), is attached. Note that the $$\hbox {FE}^2$$ method is not embarrassingly parallel since the RVE problems are coupled through the macroscopic problem. Numerical results considering different grids on both, the macroscopic and microscopic level as well as weak scaling results for up to a million parallel processes are presented.