Multigrid Preconditioner Research Articles

Multigrid Preconditioning for Discontinuous Galerkin Discretizations of an Elliptic Optimal Control Problem with a Convection-Dominated State Equation

Multigrid Preconditioning for Discontinuous Galerkin Discretizations of an Elliptic Optimal Control Problem with a Convection-Dominated State Equation

Simulation of hydro-deformation coupling problem in unsaturated porous media using exponential SWCC and hybrid improved iteration method with multigrid and multistep preconditioner

Simulation of hydro-deformation coupling problem in unsaturated porous media using exponential SWCC and hybrid improved iteration method with multigrid and multistep preconditioner

Fast 3D solvers for interactive computational mechanics

While interactive simulations have been mostly limited to Computer Graphics applications, new generations of Graphics Processing Units (GPUs) allow the realization of industrial-grade interactive 3D physics simulations. By combining an immersed boundary method with efficient GPU-based MINRES and CG solvers using a GPU-based geometric multigrid preconditioner, we demonstrate a fast industrial 3D computational mechanics solver. The various implementation aspects - specifically how they differ from similar concepts used in the Computer Graphics community - are discussed in detail. The proposed concept opens up new classes of industrial simulation applications allowing a democratization beyond today’s expert users, from designer centric simulation to operational and service decisions based on 3D simulations. To support this, we provide various benchmark cases including a real-world study of a simulation-based service decision for a damaged gear-box mount.

Multigrid Preconditioning for Regularized Least-Squares Problems

Multigrid Preconditioning for Regularized Least-Squares Problems

Gauge-equivariant pooling layers for preconditioners in lattice QCD

We demonstrate that gauge-equivariant pooling and unpooling layers can perform as well as traditional restriction and prolongation layers in multigrid preconditioner models for lattice QCD. These layers introduce a gauge degree of freedom on the coarse grid, allowing for the use of explicitly gauge-equivariant layers on the coarse grid. We investigate the construction of coarse-grid gauge fields and study their efficiency in the preconditioner model. We show that a combined multigrid neural network using a Galerkin construction for the coarse-grid gauge field eliminates critical slowing down. Published by the American Physical Society 2024

A symmetric multigrid-preconditioned Krylov subspace solver for Stokes equations

Numerical solution of discrete PDEs corresponding to saddle point problems is highly relevant to physical systems such as Stokes flow. However, scaling up numerical solvers for such systems is often met with challenges in efficiency and convergence. Multigrid is an approach with excellent applicability to elliptic problems such as the Stokes equations, and can be a solution to such challenges of scalability and efficiency. The degree of success of such methods, however, is highly contingent on the design of key components of the multigrid scheme, including the hierarchy of discretizations, and the relaxation scheme used. Additionally, in many practical cases, it may be more effective to use a multigrid scheme as a preconditioner to an iterative Krylov subspace solver, as opposed to striving for maximum efficacy of the relaxation scheme in all foreseeable settings. In this paper, we propose an efficient symmetric multigrid preconditioner for the Stokes Equations on a staggered finite difference discretization. Our contribution is focused on crafting a preconditioner that (a) is symmetric indefinite, matching the property of the Stokes system itself, (b) is appropriate for preconditioning the SQMR iterative scheme [1], and (c) has the requisite symmetry properties to be used in this context. In addition, our design is efficient in terms of computational cost and facilitates scaling to large domains.

Alya toward exascale: algorithmic scalability using PSCToolkit

In this paper, we describe an upgrade of the Alya code with up-to-date parallel linear solvers capable of achieving reliability, efficiency and scalability in the computation of the pressure field at each time step of the numerical procedure for solving a Large Eddy Simulation formulation of the incompressible Navier–Stokes equations. We developed a software module in the Alya’s kernel to interface the libraries included in the current version of PSCToolkit, a framework for the iterative solution of sparse linear systems, on parallel distributed-memory computers, by Krylov methods coupled to Algebraic MultiGrid preconditioners. The Toolkit has undergone various extensions within the EoCoE-II project with the primary goal of facing the exascale challenge. Results on a realistic benchmark for airflow simulations in wind farm applications show that the PSCToolkit solvers significantly outperform the original versions of the Conjugate Gradient method available in the Alya’s kernel in terms of scalability and parallel efficiency and represent a very promising software layer to move the Alya code toward exascale.

One-shot procedures for efficient minimum compliance topology optimization

In this paper, a one-shot approach for minimum compliance topology optimization is investigated. In the convex case of variable thickness sheet optimization, an optimality criteria scheme with a single step of an iterative state solver converges to the same solutions as an accurate procedure. This remarkable behavior is explained by the crucial role of the geometric multigrid preconditioner that generates accurate design sensitivities on a coarse scale. When applied to SIMP-based topology optimization, a one-shot procedure yields designs with the same primary load-transferring features as an accurate procedure. The difference is in the absence and presence of thin features that have a relatively minor impact on the objective. Several practical remedies are suggested that can close the small performance gaps between the one-shot procedure and an accurate one. The solution scheme is directly applicable to large-scale problems executed on parallel computers, allowing for significant computational savings.

A comparison of Algebraic Multigrid Bidomain solvers on hybrid CPU–GPU architectures

The numerical simulation of cardiac electrophysiology is a highly challenging problem in scientific computing. The Bidomain system is the most complete mathematical model of cardiac bioelectrical activity. It consists of an elliptic and a parabolic partial differential equation (PDE), of reaction–diffusion type, describing the spread of electrical excitation in the cardiac tissue. The two PDEs are coupled with a stiff system of ordinary differential equations (ODEs), representing ionic currents through the cardiac membrane. Developing efficient and scalable preconditioners for the linear systems arising from the discretization of such computationally challenging model is crucial in order to reduce the computational costs required by the numerical simulations of cardiac electrophysiology. In this work, focusing on the Bidomain system as a model problem, we have benchmarked two popular implementations of the Algebraic Multigrid (AMG) preconditioner embedded in the PETSc library and we have studied the performance on the calibration of specific parameters. We have conducted our analysis on modern HPC architectures, performing scalability tests on multi-core and multi-GPUs settings. The results have shown that, for our problem, although scalability is verified on CPUs, GPUs are the optimal choice, since they yield the best performance in terms of solution time.

Reconstructing relaxed configurations in elastic bodies: Mathematical formulations and numerical methods for cardiac modeling

Modeling the behavior of biological tissues and organs often necessitates the knowledge of their shape in the absence of external loads. However, when their geometry is acquired in-vivo through imaging techniques, bodies are typically subject to mechanical deformation due to the presence of external forces, and the load-free configuration needs to be reconstructed. This paper addresses this crucial and frequently overlooked topic, known as the inverse elasticity problem (IEP), by delving into both theoretical and numerical aspects, with a particular focus on cardiac mechanics. In this work, we extend Shield’s seminal work to determine the structure of the IEP with arbitrary material inhomogeneities and in the presence of both body and active forces. These aspects are fundamental in computational cardiology, and we show that they may break the variational structure of the inverse problem. In addition, we show that the inverse problem might have no solution even in the presence of constant Neumann boundary conditions and a polyconvex strain energy functional. We then present the results of extensive numerical tests to validate our theoretical framework, and to characterize the computational challenges associated with a direct numerical approximation of the IEP. Specifically, we show that this framework outperforms existing approaches both in terms of robustness and optimality, such as Sellier’s iterative procedure, even when the latter is improved with acceleration techniques. A notable discovery is that multigrid preconditioners are, in contrast to standard elasticity, not efficient, where a one-level additive Schwarz and generalized Dryja–Smith–Widlund provide a much more reliable alternative. Finally, we successfully address the IEP for a full-heart geometry, demonstrating that the IEP formulation can compute the stress-free configuration in real-life scenarios where Sellier’s algorithm proves inadequate.

Relaxed Jacobi method as multigrid smoother and preconditioner

Relaxed Jacobi method as multigrid smoother and preconditioner

Multigrid for two-sided fractional differential equations discretized by finite volume elements on graded meshes

It is known that the solution of a conservative steady-state two-sided fractional diffusion problem can exhibit singularities near the boundaries. As a consequence of this, and due to the conservative nature of the problem, we adopt a finite volume elements discretization approach over a generic non-uniform mesh. We focus on grids mapped by a smooth function which consists in a combination of a graded mesh near the singularity and a uniform mesh where the solution is smooth. Such a choice gives rise to Toeplitz-like discretization matrices and thus allows a low computational cost of the matrix–vector product and detailed spectral analysis. The obtained spectral information is used to develop an ad-hoc parameter-free multigrid preconditioner for GMRES, which is numerically shown to yield good convergence results in presence of graded meshes mapped by power functions that accumulate points near the singularity. The approximation order of the considered graded meshes is numerically compared with the one of a certain composite mesh given in literature that still leads to Toeplitz-like linear systems and is then still well-suited for our multigrid method. Several numerical tests confirm that power-graded meshes result in lower approximation errors than composite ones and that our solver has a wide range of applicability.

A Primal-Dual Box-Constrained QP Pressure Poisson Solver With Topology-Aware Geometry-Inspired Aggregation AMG.

We propose a new barrier-based box-constrained convex QP solver based on a primal-dual interior point method to efficiently solve large-scale pressure Poisson problems with non-negative pressure constraints, which commonly arise in liquid animation. The performance of prior active-set-based approaches is limited by the need to repeatedly update the active set. Our solver eliminates this issue by entirely avoiding the use of an active set, which in turn makes the inner problems of our Newton iteration process fully unconstrained. For efficiency, exploiting the solution uniqueness of convex QPs and the fact that the pressure constraints are simple box constraints, we aggressively update solution candidates without performing any step selection procedure (such as line search) and instead directly clamp candidates back to the bounds wherever constraint violations occur. Additionally, to accelerate the inner linear solves, we present a topology-aware geometry-inspired aggregation algebraic multigrid preconditioner and describe in detail several key performance optimizations that we incorporate. We demonstrate the efficacy of our solver in various practical scenarios and show that it often surpasses various alternatives in terms of speed and memory usage.

A hybrid shifted Laplacian multigrid and domain decomposition preconditioner for the elastic Helmholtz equations

In this work we extend the shifted Laplacian approach to the elastic Helmholtz equation. The shifted Laplacian multigrid method is a common preconditioning approach for the discretized acoustic Helmholtz equation. In some cases, like geophysical seismic imaging, one needs to consider the elastic Helmholtz equation, which is harder to solve: it is three times larger and contains a nullity-rich grad-div term. These properties make the solution of the equation more difficult for multigrid solvers. The key idea in this work is combining the shifted Laplacian with approaches for linear elasticity. We provide local Fourier analysis and numerical evidence that the convergence rate of our method is independent of the Poisson's ratio. Moreover, to better handle the problem size, we complement our multigrid method with the domain decomposition approach, which works in synergy with the local nature of the shifted Laplacian, so we enjoy the advantages of both methods without sacrificing performance. We demonstrate the efficiency of our solver on 2D and 3D problems in heterogeneous media.

COSMOGLOBE: Towards end-to-end CMB cosmological parameter estimation without likelihood approximations

We implement support for a cosmological parameter estimation algorithm in Commander and quantify its computational efficiency and cost. For a semi-realistic simulation similar to Planck LFI 70 GHz, we find that the computational cost of producing one single sample is about 20 CPU-hours and that the typical Markov chain correlation length is ∼100 samples. The net effective cost per independent sample is ∼2000 CPU-hours, in comparison with all low-level processing costs of 812 CPU-hours for Planck LFI and WMAP in COSMOGLOBE Data Release 1. Thus, although technically possible to run already in its current state, future work should aim to reduce the effective cost per independent sample by one order of magnitude to avoid excessive runtimes, for instance through multi-grid preconditioners and/or derivative-based Markov chain sampling schemes. This work demonstrates the computational feasibility of true Bayesian cosmological parameter estimation with end-to-end error propagation for high-precision CMB experiments without likelihood approximations, but it also highlights the need for additional optimizations before it is ready for full production-level analysis.

Targeting a faster time‐to‐solution of mortar‐based contact problems

AbstractWhile mortar methods are undoubtedly among the most prominent and accurate approaches in computational contact mechanics, they come with considerable cost. The required numerical effort originates not only from the evaluation of the mortar integrals but also from the solution of the system of linear equations within Newton‐type solvers. While each of these cost factors has been analyzed and addressed individually in the past, this contribution outlines approaches and first results to comprehensively reduce the time to solution by combining dynamic load balancing techniques with algebraic multigrid preconditioners. We present a proof of concept for a three‐dimensional unilateral contact simulation with an evolving contact zone.

A comparative study of scalable multilevel preconditioners for cardiac mechanics

In this work, we provide a performance comparison between the Balancing Domain Decomposition by Constraints (BDDC) and the Algebraic Multigrid (AMG) preconditioners for cardiac mechanics on both structured and unstructured finite element meshes. The mechanical behavior of myocardium can be described by the equations of three-dimensional finite elasticity, which are discretized by finite elements in space and yield the solution of a large scale nonlinear algebraic system. This problem is solved by a Newton-Krylov method, where the solution of the Jacobian linear system is accelerated by BDDC/AMG preconditioners. We thoroughly explore the main parameters of the BDDC preconditioner in order to make the comparison fair. We focus on: the performance of different direct solvers for the local and coarse problems of the BDDC algorithm; the impact of the different choices of BDDC primal degrees of freedom; and the influence of the finite element degree. Scalability tests are performed on Linux clusters up to 1024 processors, and we conclude with a performance study on a realistic electromechanical simulation.

A Multi-GPU Aggregation-Based AMG Preconditioner for Iterative Linear Solvers

We present and release in open source format a sparse linear solver which efficiently exploits heterogeneous parallel computers. The solver can be easily integrated into scientific applications that need to solve large and sparse linear systems on modern parallel computers made of hybrid nodes hosting Nvidia Graphics Processing Unit (GPU) accelerators. The work extends previous efforts of some of the authors in the exploitation of a single GPU accelerator and proposes an implementation, based on the hybrid MPI-CUDA software environment, of a Krylov-type linear solver relying on an efficient Algebraic MultiGrid (AMG) preconditioner already available in the <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">BootCMatchG</i> library. Our design for the hybrid implementation has been driven by the best practices for minimizing data communication overhead when multiple GPUs are employed, yet preserving the efficiency of the GPU kernels. Strong and weak scalability results of the new version of the library on well-known benchmark test cases are discussed. Comparisons with the Nvidia AmgX solution show a speedup, in the solve phase, up to 2.0x.

Algebraic Multigrid Preconditioner for Statically Condensed Systems Arising from Lowest-Order Hybrid Discretizations

We address the numerical solution of linear systems arising from the hybrid discretizations of second-order elliptic partial differential equations. Such discretizations hinge on a hybrid set of degrees of freedom (DoFs), respectively, defined in cells and faces, which naturally gives rise to a global hybrid system of linear equations. Assuming that the cell unknowns are only locally coupled, they can be efficiently eliminated from the system, leaving only face unknowns in the resulting Schur complement, which is also called the statically condensed matrix. We propose in this work an algebraic multigrid (AMG) preconditioner specifically targeting condensed systems corresponding to lowest-order discretizations (piecewise constant). Like traditional AMG methods, we retrieve geometric information on the coupling of the DoFs from algebraic data. However, as the condensed matrix only gives information on the faces, we use the uncondensed version to reconstruct the connectivity graph between elements and faces. An aggregation-based coarsening strategy mimicking a geometric coarsening or semicoarsening can then be set up to build coarse levels. Numerical experiments are performed on diffusion problems discretized by the hybrid high-order method at the lowest order. Our approach uses a K-cycle to precondition an outer flexible Krylov method. The results demonstrate similar performances, in most cases, compared to a standard AMG method and a notable improvement on anisotropic problems with Cartesian meshes.

End-to-end GPU acceleration of low-order-refined preconditioning for high-order finite element discretizations

In this article, we present algorithms and implementations for the end-to-end GPU acceleration of matrix-free low-order-refined preconditioning of high-order finite element problems. The methods described here allow for the construction of effective preconditioners for high-order problems with optimal memory usage and computational complexity. The preconditioners are based on the construction of a spectrally equivalent low-order discretization on a refined mesh, which is then amenable to, for example, algebraic multigrid preconditioning. The constants of equivalence are independent of mesh size and polynomial degree. For vector finite element problems in H(curl) and H(div) (e.g., for electromagnetic or radiation diffusion problems), a specially constructed interpolation–histopolation basis is used to ensure fast convergence. Detailed performance studies are carried out to analyze the efficiency of the GPU algorithms. The kernel throughput of each of the main algorithmic components is measured, and the strong and weak parallel scalability of the methods is demonstrated. The different relative weighting and significance of the algorithmic components on GPUs and CPUs is discussed. Results on problems involving adaptively refined nonconforming meshes are shown, and the use of the preconditioners on a large-scale magnetic diffusion problem using all spaces of the finite element de Rham complex is illustrated.