Algebraic Multigrid Preconditioner Research Articles

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.

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.

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.

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.

A matrix–free high–order solver for the numerical solution of cardiac electrophysiology

We propose a matrix–free solver for the numerical solution of the cardiac electrophysiology model consisting of the monodomain nonlinear reaction–diffusion equation coupled with a system of ordinary differential equations for the ionic species. Our numerical approximation is based on the high–order Spectral Element Method (SEM) to achieve accurate numerical discretization while employing a much smaller number of Degrees of Freedom than first–order Finite Elements. We combine vectorization with sum–factorization, thus allowing for a very efficient use of high–order polynomials in a high performance computing framework. We validate the effectiveness of our matrix–free solver in a variety of applications and perform different electrophysiological simulations ranging from a simple slab of cardiac tissue to a realistic four–chamber heart geometry. We compare SEM to SEM with Numerical Integration (SEM–NI), showing that they provide comparable results in terms of accuracy and efficiency. In both cases, increasing the local polynomial degree p leads to better numerical results and smaller computational times than reducing the mesh size h. We also implement a matrix–free Geometric Multigrid preconditioner that results in a comparable number of linear solver iterations with respect to a state–of–the–art matrix–based Algebraic Multigrid preconditioner. As a matter of fact, the matrix–free solver proposed here yields up to 45× speed–up with respect to a conventional matrix–based solver.

Integral boundary conditions in phase field models

Modeling the chemical, electric and thermal transport as well as phase transitions and the accompanying mesoscale microstructure evolution within a material in an electronic device setting involves the solution of partial differential equations often with integral boundary conditions. Employing the familiar Poisson equation describing the electric potential evolution in a material exhibiting insulator to metal transitions, we exploit a special property of such an integral boundary condition, and we properly formulate the variational problem and establish its well-posedness. We then compare our method with the commonly-used Lagrange multiplier method that can also handle such boundary conditions. Numerical experiments demonstrate that our new method achieves optimal convergence rate in contrast to the conventional Lagrange multiplier method. Furthermore, the linear system derived from our method is symmetric positive definite, and can be efficiently solved by Conjugate Gradient method with algebraic multigrid preconditioning.

Scalable computational algorithms for geospatial COVID-19 spread using high performance computing.

A nonlinear partial differential equation (PDE) based compartmental model of COVID-19 provides a continuous trace of infection over space and time. Finer resolutions in the spatial discretization, the inclusion of additional model compartments and model stratifications based on clinically relevant categories contribute to an increase in the number of unknowns to the order of millions. We adopt a parallel scalable solver that permits faster solutions for these high fidelity models. The solver combines domain decomposition and algebraic multigrid preconditioners at multiple levels to achieve the desired strong and weak scalabilities. As a numerical illustration of this general methodology, a five-compartment susceptible-exposed-infected-recovered-deceased (SEIRD) model of COVID-19 is used to demonstrate the scalability and effectiveness of the proposed solver for a large geographical domain (Southern Ontario). It is possible to predict the infections for a period of three months for a system size of 186 million (using 3200 processes) within 12 hours saving months of computational effort needed for the conventional solvers.

An efficient semi-implicit solver for solid electrolyte interphase growth in Li-ion batteries

Growth of the Solid Electrolyte Interphase is one of the major degradation processes in Li-ion batteries and yields an increase in cell resistance, as well as a reduction in capacity. In order to facilitate computer-aided design of novel cell concepts and omit costly prototyping, it is crucial to incorporate this mechanism into electrochemical simulations on the microscale and to develop efficient numerical solvers. Electrochemical battery models without degradation are commonly discretized with a fully implicit discretization resulting in a monolithic scheme. Solving the resulting systems with a Newton method allows for quadratic convergence, where the major part of the computational effort is spent on solving the underlying linear systems. Hence, finding a suitable preconditioner for the linear solver is essential for overall solver performance. While an efficient algebraic multigrid preconditioner results in a good performance for the model without degradation, we encounter a noticeable loss in performance of the solver when including the coupled degradation model. We attribute the observed performance issues to the implicitly defined interface currents introduced by the degradation model. By taking into account the slow dynamics of the degradation, we propose an alternative semi-implicit solution approach separating the degradation dynamics and eliminating the implicit interface currents, in order to promote the efficient utilization of the algebraic multigrid preconditioner. The monolithic and semi-implicit solvers are evaluated by performing a series of full cell simulations featuring calendar and cyclic aging scenarios and a complex geometry including the binder domain in each electrode. Based on the results, the semi-implicit solver significantly improves the performance, while maintaining the convergence behaviour.

A Comparison of Linear Solvers for Resolving Flow in Three‐Dimensional Discrete Fracture Networks

AbstractWe compare various methods for resolving steady flow within three‐dimensional discrete fracture networks, including direct methods, Krylov subspace methods with and without preconditioning, and multi‐grid methods. We compared the performance of the methods based on compute times and scaling of the solution as a function of the number of grid nodes and log‐variance of the hydraulic aperture. The methods are applied to three test cases: (a) variable density of networks with a truncated power‐law distribution of fracture lengths, (b) a fixed network composed of monodisperse fracture sizes but varied permeability/aperture heterogeneity, (c) and a network based on field site in Nevada, US. We chose these cases to allow us to study the impact of the mesh size and flow properties, as well as to demonstrate our conclusions on a large‐scale, realistic problem (more than 40 million mesh nodes). A direct solution using Cholesky factorization outperformed other methods for every example but was closely followed in performance by some algebraic multigrid (AMG) preconditioned Krylov subspace methods. Among the Krylov methods, conjugate gradients (CG) with an AMG preconditioner performs the best. Generally, Cholesky factorization is recommended, but CG with an AMG preconditioner may be suitable for very large problems beyond 40 million nodes where the entire linear system cannot reside in memory.

A Monolithic Divergence-Conforming HDG Scheme for a Linear Fluid-Structure Interaction Model

We present a novel monolithic divergence-conforming HDG scheme for a linear fluid-structure interaction problem with a thick structure. A pressure-robust optimal energy-norm estimate is obtained for the semidiscrete scheme. When combined with a Crank--Nicolson time discretization, our fully discrete scheme is energy stable and produces an exactly divergence-free fluid velocity approximation. The resulting linear system, which is symmetric and indefinite, is solved using a preconditioned MinRes method with a robust block algebraic multigrid preconditioner.

A monolithic algebraic multigrid framework for multiphysics applications with examples from resistive MHD

We consider monolithic algebraic multigrid (AMG) algorithms for the solution of block linear systems arising from multiphysics simulations. While the multigrid idea is applied directly to the entire linear system, AMG operators are constructed by leveraging the matrix block structure. In particular, each block corresponds to a set of physical unknowns and physical equations. Multigrid components are constructed by first applying existing AMG procedures to matrix sub-blocks. The resulting AMG sub-components are then composed together to define a monolithic AMG preconditioner. Given the problem-dependent nature of multiphysics systems, different blocking choices may work best in different situations, and so software flexibility is essential. We apply different blocking strategies to systems arising from resistive magnetohydrodynamics in order to demonstrate the associated trade-offs.

Partial Reuse AMG Setup Cost Amortization Strategy for the Solution of Non-Steady State Problems

The partial reuse algebraic multigrid (AMG) setup cost amortization strategy is presented for the solution of non-steady state problems. The transfer operators are reused from the previous time steps, and the system matrices and the smoother operators are rebuilt on each of the AMG hierarchy levels. It is shown on the example of modelling a two-fluid dam break scenario that the strategy may decrease the AMG preconditioner setup cost by 40% to 200%. The total compute time is decreased by up to 20%, but the specific outcome depends on the fraction of time that the setup step initially takes.

Comparison between algebraic and matrix‐free geometric multigrid for a Stokes problem on adaptive meshes with variable viscosity

AbstractProblems arising in Earth's mantle convection involve finding the solution to Stokes systems with large viscosity contrasts. These systems contain localized features which, even with adaptive mesh refinement, result in linear systems that can be on the order of 109 or more unknowns. One common approach for preconditioning to the velocity block of these systems is to apply an Algebraic Multigrid (AMG) V‐cycle (as is done in the ASPECT software, for example), however, we find that AMG is lacking robustness with respect to problem size and number of parallel processes. Additionally, we see an increase in iteration counts with refinement when using AMG. In contrast, the Geometric Multigrid (GMG) method, by using information about the geometry of the problem, should offer a more robust option.Here we present a matrix‐free GMG V‐cycle which works on adaptively refined, distributed meshes, and we will compare it against the current AMG preconditioner (Trilinos ML) used in the ASPECT1 software. We will demonstrate the robustness of GMG with respect to problem size and show scaling up to 114,688 cores and 217 billion unknowns. All computations are run using the open‐source, finite element library deal.II.2

AMG Preconditioners for Linear Solvers towards Extreme Scale

Linear solvers for large and sparse systems are a key element of scientific applications, and their efficient implementation is necessary to harness the computational power of current computers. Algebraic MultiGrid (AMG) preconditioners are a popular ingredient of such linear solvers; this is the motivation for the present work where we examine some recent developments in a package of AMG preconditioners to improve efficiency, scalability, and robustness on extreme-scale problems. The main novelty is the design and implementation of a parallel coarsening algorithm based on aggregation of unknowns employing weighted graph matching techniques; this is a completely automated procedure, requiring no information from the user, and applicable to general symmetric positive definite (s.p.d.) matrices. The new coarsening algorithm improves in terms of numerical scalability at low operator complexity over decoupled aggregation algorithms available in previous releases of the package. The preconditioners package is built on the parallel software framework \texttt{PSBLAS}, which has also been updated to progress towards exascale. We present weak scalability results on one of the most powerful supercomputers in Europe, for linear systems with sizes up to $O(10^{10})$ unknowns.

On the use of multigrid preconditioners for topology optimization

Topology optimization for large-scale problems continues to be a computational challenge. Several works exist in the literature to address this topic, and all make use of iterative solvers to handle the linear system arising from the finite element analysis (FEA). However, the preconditioners used in these works vary, and in many cases are notably suboptimal. A handful of works have already demonstrated the effectiveness of geometric multigrid (GMG) preconditioners in topology optimization. We provide a direct comparison of GMG preconditioners with algebraic multigrid (AMG) preconditioners. We demonstrate that AMG preconditioners offer improved robustness over GMG preconditioners as topologies evolve, albeit with a higher overhead cost. In 2D the gain from increased robustness more than offsets the overhead cost. However, in 3D the overhead becomes prohibitively large. We thus demonstrate the benefits of mixing geometric and algebraic methods to limit overhead cost while improving robustness, particularly in 3D.

Preconditioning strategies for vectorial finite element linear systems arising from phase-field models for fracture mechanics

Phase-field models are frequently adopted to simulate fracture mechanics problems in the context of the finite element method. To depict fracture, this method involves solving a coupled set of Helmholtz-like damage-field equation and augmented linear momentum balance equation. Solutions to these coupled equations are then used as descriptions of crack propagation phenomena within solids. However, this method imposes a constrain of using extremely fine meshing for properly predicting cracks. For practical problems of interest, this very often leads to linear systems with large sizes that have to be repetitively assembled and solved. As such, iterative solution procedures such as the Krylov subspace-based methods for solving these large linear systems within the framework of serial/parallel computing environments become mandatory to obtain results in a feasible time. In this work, the vectorial finite discretization for a hybrid phase-field formulation – a monolithic solving scheme – is presented. The underlying nonlinearity present in the coupled set of equations of the hybrid phase-field model is dealt through Picard iteration that helps to preserve the symmetry of the linearized system to solve. Due to the symmetric positive definite nature of the finite element linear systems obtained for this problem, the conjugate gradient method makes a standard choice of iterative solution algorithm. In this article, to improve convergence rates, consequently time to solution, of the conjugate gradient method applied to crack propagation problems, different preconditioning strategies are analyzed, tuned, and discussed. Brittle fracture benchmarks are used to measure the performance of preconditioners which are then applied to massively parallel simulations with millions of unknowns. A series of numerical experiments show that the algebraic multigrid preconditioner is well suited for solving the phase-field model for fracture, being superior to the Jacobi and the block Jacobi preconditioning in all regards: ease of solving the problem, iterations to converge, time to solution, and parallel scaling on more than a thousand processes.