Geometric Multigrid Preconditioner Research Articles

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.

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.

Optimal Geometric Multigrid Preconditioners for HDG-P0 Schemes for the reaction-diffusion equation and the Generalized Stokes equations

We present the lowest-order hybridizable discontinuous Galerkin schemes with numerical integration (quadrature), denoted as HDG-P0 for the reaction-diffusion equation and the generalized Stokes equations on conforming simplicial meshes in two- and three-dimensions. Here by lowest order, we mean that the (hybrid) finite element space for the global HDG facet degrees of freedom (DOFs) is the space of piecewise constants on the mesh skeleton. A discontinuous piecewise linear space is used for the approximation of the local primal unknowns. We give the optimal a priori error analysis of the proposed HDG-P0 schemes, which hasn’t appeared in the literature yet for HDG discretizations as far as numerical integration is concerned. Moreover, we propose optimal geometric multigrid preconditioners for the statically condensed HDG-P0 linear systems on conforming simplicial meshes. In both cases, we first establish the equivalence of the statically condensed HDG system with a (slightly modified) nonconforming Crouzeix–Raviart (CR) discretization, where the global (piecewise-constant) HDG finite element space on the mesh skeleton has a natural one-to-one correspondence to the nonconforming CR (piecewise-linear) finite element space that live on the whole mesh. This equivalence then allows us to use the well-established nonconforming geometry multigrid theory to precondition the condensed HDG system. Numerical results in two- and three-dimensions are presented to verify our theoretical findings.

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.

A multigrid preconditioner for spatially adaptive high-order meshless method on fluid–solid interaction problems

We present a monolithic geometric multigrid preconditioner for solving fluid–solid interaction problems in Stokes limit. The problems are discretized by a spatially adaptive high-order meshless method, the generalized moving least squares (GMLS) with adaptive h-refinement. For solving fluid–solid interaction problems, we need to deal with a tightly coupled system consisting of the flow field and solid bodies, resulting in a linear system of equations with a block structure. In Stokes limit, solid kinematics can be dominated by the singularities governing the lubrication effects. Resolving those singularities with adaptive h-refinement can lead to an ill-conditioned linear system of equations. The key ingredients of the multigrid preconditioner include the interpolation and restriction operators and the smoothers. For constructing the interpolation and restriction operators, we utilize the geometric information of hierarchical sets of GMLS nodes generated in adaptive h-refinement. We build decoupled smoothers through physics-based splitting and then combine them via a multiplicative overlapping Schwarz approach. Through numerical examples with the inclusion of different numbers and shapes of solid bodies, we demonstrate the performance and assess the scalability of the designed preconditioner. As the total degrees of freedom and the number of solid bodies increase, the proposed monolithic geometric multigrid preconditioner can ensure convergence and good scalability when using the Krylov iterative method for solving the linear systems of equations generated from the spatially adaptive GMLS discretization. More specifically, for a fixed number of solid bodies, as the discretization resolution is incrementally refined, the number of iterations of the linear solver can be maintained at the same level, indicating nearly linear scalability of our preconditioner with respect to the total degrees of freedom. When the number of solid bodies Ns increases, the number of iterations is nearly proportional to Ns, implying the sublinear optimality with respect to the number of solid bodies.

A multigrid preconditioner for tensor product spline smoothing

Penalized spline smoothing is a well-established, nonparametric regression method that is efficient for one and two covariates. Its extension to more than two covariates is straightforward but suffers from exponentially increasing memory demands and computational complexity, which brings the method to its numerical limit. Penalized spline smoothing with multiple covariates requires solving a large-scale, regularized least-squares problem where the occurring matrices do not fit into storage of common computer systems. To overcome this restriction, we introduce a matrix-free implementation of the conjugate gradient method. We further present a matrix-free implementation of a simple diagonal as well as more advanced geometric multigrid preconditioner to significantly speed up convergence of the conjugate gradient method. All algorithms require a negligible amount of memory and therefore allow for penalized spline smoothing with multiple covariates. Moreover, for arbitrary but fixed covariate dimension, we show grid independent convergence of the multigrid preconditioner which is fundamental to achieve algorithmic scalability.

A high-throughput hybrid task and data parallel Poisson solver for large-scale simulations of incompressible turbulent flows on distributed GPUs

The solution of the pressure Poisson equation arising in the numerical solution of incompressible Navier–Stokes equations (INSE) is by far the most expensive part of the computational procedure, and often the major restricting factor for parallel implementations. Improvements in iterative linear solvers, e.g. deploying Krylov-based techniques and multigrid preconditioners, have been successfully applied for solving the INSE on CPU-based parallel computers. These numerical schemes, however, do not necessarily perform well on GPUs, mainly due to differences in the hardware architecture. Our previous work using many P100 GPUs of a flagship supercomputer showed that porting a highly optimized MPI-parallel CPU-based INSE solver to GPUs, accelerates significantly the underlying numerical algorithms, while the overall acceleration remains limited (Zolfaghari et al. [3]). The performance loss was mainly due to the Poisson solver, particularly the V-cycle geometric multigrid preconditioner. We also observed that the pure compute time for the GPU kernels remained nearly constant as grid size was increased. Motivated by these observations, we present herein an algebraically simpler, yet more advanced parallel implementation for the solution of the Poisson problem on large numbers of distributed GPUs. Data parallelism is achieved by using the classical Jacobi method with successive over-relaxation and an optimized iterative driver routine. Task parallelism is enhanced via minimizing GPU-GPU data exchanges as iterations proceed to reduce the communication overhead. The hybrid parallelism results in nearly 300 times less time-to-solution and thus computational cost (measured in node-hours) for the Poisson problem, compared to our best-case scenario CPU-based parallel implementation which uses a preconditioned BiCGstab method. The Poisson solver is then embedded in a flow solver with explicit third-order Runge-Kutta scheme for time-integration, which has been previously ported to GPUs. The flow solver is validated and computationally benchmarked for the transition and decay of the Taylor-Green Vortex at Re=1600 and the flow around a solid sphere at ReD=3700. Good strong scaling is demonstrated for both benchmarks. Further, nearly 70% lower electrical energy consumption than the CPU implementation is reported for Taylor-Green vortex case. We finally deploy the solver for DNS of systolic flow in a bileaflet mechanical heart valve, and present new insight into the complex laminar-turbulent transition process in this prosthesis.

Phase-field modeling of multivariant martensitic transformation at finite-strain: Computational aspects and large-scale finite-element simulations

Large-scale 3D martensitic microstructure evolution problems are studied using a finite-element discretization of a finite-strain phase-field model. The model admits an arbitrary crystallography of transformation and arbitrary elastic anisotropy of the phases, and incorporates Hencky-type elasticity, a penalty-regularized double-obstacle potential, and viscous dissipation. The finite-element discretization of the model is performed in Firedrake and relies on the PETSc solver library. The large systems of linear equations arising are efficiently solved using GMRES and a geometric multigrid preconditioner with a carefully chosen relaxation. The modeling capabilities are illustrated through a 3D simulation of the microstructure evolution in a pseudoelastic CuAlNi single crystal during nano-indentation, with all six orthorhombic martensite variants taken into account. Robustness and a good parallel scaling performance have been demonstrated, with the problem size reaching 150 million degrees of freedom.

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.

Algebraic multigrid preconditioning of the Hessian in optimization constrained by a partial differential equation

SummaryWe construct an algebraic multigrid (AMG) based preconditioner for the reduced Hessian of a linear‐quadratic optimization problem constrained by an elliptic partial differential equation. While the preconditioner generalizes a geometric multigrid preconditioner introduced in earlier works, its construction relies entirely on a standard AMG infrastructure built for solving the forward elliptic equation, thus allowing for it to be implemented using a variety of AMG methods and standard packages. Our analysis establishes a clear connection between the quality of the preconditioner and the AMG method used. The proposed strategy has a broad and robust applicability to problems with unstructured grids, complex geometry, and varying coefficients. The method is implemented using the Hypre package and several numerical examples are presented.

A matrix‐free approach for finite‐strain hyperelastic problems using geometric multigrid

SummaryThis work investigates matrix‐free algorithms for problems in quasi‐static finite‐strain hyperelasticity. Iterative solvers with matrix‐free operator evaluation have emerged as an attractive alternative to sparse matrices in the fluid dynamics and wave propagation communities because they significantly reduce the memory traffic, the limiting factor in classical finite element solvers. Specifically, we study different matrix‐free realizations of the finite element tangent operator and determine whether generalized methods of incorporating complex constitutive behavior might be feasible. In order to improve the convergence behavior of iterative solvers, we also propose a method by which to construct level tangent operators and employ them to define a geometric multigrid preconditioner. The performance of the matrix‐free operator and the geometric multigrid preconditioner is compared to the matrix‐based implementation with an algebraic multigrid (AMG) preconditioner on a single node for a representative numerical example of a heterogeneous hyperelastic material in two and three dimensions. We find that matrix‐free methods for finite‐strain solid mechanics are very promising, outperforming linear matrix‐based schemes by two to five times, and that it is possible to develop numerically efficient implementations that are independent of the hyperelastic constitutive law.

Multigrid solvers for immersed finite element methods and immersed isogeometric analysis

Ill-conditioning of the system matrix is a well-known complication in immersed finite element methods and trimmed isogeometric analysis. Elements with small intersections with the physical domain yield problematic eigenvalues in the system matrix, which generally degrades efficiency and robustness of iterative solvers. In this contribution we investigate the spectral properties of immersed finite element systems treated by Schwarz-type methods, to establish the suitability of these as smoothers in a multigrid method. Based on this investigation we develop a geometric multigrid preconditioner for immersed finite element methods, which provides mesh-independent and cut-element-independent convergence rates. This preconditioning technique is applicable to higher-order discretizations, and enables solving large-scale immersed systems at a computational cost that scales linearly with the number of degrees of freedom. The performance of the preconditioner is demonstrated for conventional Lagrange basis functions and for isogeometric discretizations with both uniform B-splines and locally refined approximations based on truncated hierarchical B-splines.

Camellia: A Rapid Development Framework for Finite Element Solvers

Abstract The discontinuous Petrov–Galerkin (DPG) methodology of Demkowicz and Gopalakrishnan guarantees the optimality of the finite element solution in a user-controllable energy norm, and provides several features supporting adaptive schemes. The approach provides stability automatically; there is no need for carefully derived numerical fluxes (as in DG schemes) or for mesh-dependent stabilization terms (as in stabilized methods). In this paper, we focus on features of Camellia that facilitate implementation of new DPG formulations; chief among these is a rich set of features in support of symbolic manipulation, which allow, e.g., bilinear formulations in the code to appear much as they would on paper. Many of these features are general in the sense that they can also be used in the implementation of other finite element formulations. In fact, because DPG’s requirements are essentially a superset of those of other finite element methods, Camellia provides built-in support for most common methods. We believe, however, that the combination of an essentially “hands-free” finite element methodology as found in DPG with the rapid development features of Camellia are particularly winsome, so we focus on use cases in this class. In addition to the symbolic manipulation features mentioned above, Camellia offers support for one-irregular adaptive meshes in 1D, 2D, 3D, and space-time. It provides a geometric multigrid preconditioner particularly suited for DPG problems, and supports distributed parallel execution using MPI. For its load balancing and distributed data structures, Camellia relies on packages from the Trilinos project, which simplifies interfacing with other computational science packages. Camellia also allows loading of standard mesh formats through an interface with the MOAB package. Camellia includes support for static condensation to eliminate element-interior degrees of freedom locally, usually resulting in substantial reduction of the cost of the global problem. We include a discussion of the variational formulations built into Camellia, with references to those formulations in the literature, as well as an MPI performance study.

Space-time discretizations using constrained first-order system least squares (CFOSLS)

This paper studies finite element discretizations for three types of time-dependent PDEs, namely heat equation, scalar conservation law and wave equation, which we reformulate as first order systems in a least-squares setting, subject to a space-time conservation constraint (coming from the original PDE). Available piecewise polynomial finite element spaces in (n+1)-dimensions for functional spaces from the (n+1)-dimensional de Rham sequence for n=2,3 are used for the implementation of the method. Computational results illustrating the error behavior, iteration counts and performance of block-diagonal and monolithic geometric multigrid preconditioners are presented for the discrete CFOSLS system. The results are obtained from a parallel implementation of the methods for which we report reasonable scalability.

Matrix‐Free Locally Adaptive Finite Element Solution of Density‐Functional Theory With Nonorthogonal Orbitals and Multigrid Preconditioning

In this paper, we propose a new numerical method to find the ground state energy of a given physical system within the Kohn–Sham density functional theory. The h‐adaptive finite element method is adopted for spatial discretization and implemented with matrix‐free operator evaluation. The ground state energy is found by performing unconstrained minimization with non‐orthogonal orbitals using the limited memory Broyden–Fletcher–Goldfarb–Shanno (BFGS) method. A geometric multigrid preconditioner is applied to improve the convergence. The clear advantage of the proposed approach is demonstrated on selected examples by comparing the performance to other methods such as preconditioned steepest descent minimization. The proposed method provides a solid framework toward (N) complexity for the locally adaptive real‐space solution of density functional theory with finite elements.

Preconditioners for the Nonconforming Voxel Edge Element Method

The nonconforming voxel finite element method builds a structured grid of nested cube elements. It has been used both with scalar (nodal) and vector (edge) elements. For large problems, iterative methods such as the conjugate gradient method are needed to solve the global matrix equation. Preconditioning the iterative method can greatly decrease the computation time. A geometric multigrid preconditioner is investigated, for the scattering of plane waves by conducting objects, using either a Gauss-Seidel smoother or a smoother that incorporates a projection onto the gradient subspace.

ParMooN—A modernized program package based on mapped finite elements

ParMooN is a program package for the numerical solution of elliptic and parabolic partial differential equations. It inherits the distinct features of its predecessor MooNMD (John and Matthies, 2004): strict decoupling of geometry and finite element spaces, implementation of mapped finite elements as their definition can be found in textbooks, and a geometric multigrid preconditioner with the option to use different finite element spaces on different levels of the multigrid hierarchy. After having presented some thoughts about in-house research codes, this paper focuses on aspects of the parallelization for a distributed memory environment, which is the main novelty of ParMooN. Numerical studies, performed on compute servers, assess the efficiency of the parallelized geometric multigrid preconditioner in comparison with some parallel solvers that are available in the library PETSc. The results of these studies give a first indication whether the cumbersome implementation of the parallelized geometric multigrid method was worthwhile or not.

Boundary Treatment and Multigrid Preconditioning for Semi-Lagrangian Schemes Applied to Hamilton\u2013Jacobi\u2013Bellman Equations

We analyse two practical aspects that arise in the numerical solution of Hamilton–Jacobi–Bellman equations by a particular class of monotone approximation schemes known as semi-Lagrangian schemes. These schemes make use of a wide stencil to achieve convergence and result in discretization matrices that are less sparse and less local than those coming from standard finite difference schemes. This leads to computational difficulties not encountered there. In particular, we consider the overstepping of the domain boundary and analyse the accuracy and stability of stencil truncation. This truncation imposes a stricter CFL condition for explicit schemes in the vicinity of boundaries than in the interior, such that implicit schemes become attractive. We then study the use of geometric, algebraic and aggregation-based multigrid preconditioners to solve the resulting discretised systems from implicit time stepping schemes efficiently. Finally, we illustrate the performance of these techniques numerically for benchmark test cases from the literature.

High level implementation of geometric multigrid solvers for finite element problems: Applications in atmospheric modelling

The implementation of efficient multigrid preconditioners for elliptic partial differential equations (PDEs) is a challenge due to the complexity of the resulting algorithms and corresponding computer code. For sophisticated (mixed) finite element discretisations on unstructured grids an efficient implementation can be very time consuming and requires the programmer to have in-depth knowledge of the mathematical theory, parallel computing and optimisation techniques on manycore CPUs. In this paper we show how the development of bespoke multigrid preconditioners can be simplified significantly by using a framework which allows the expression of the each component of the algorithm at the correct abstraction level. Our approach (1) allows the expression of the finite element problem in a language which is close to the mathematical formulation of the problem, (2) guarantees the automatic generation and efficient execution of parallel optimised low-level computer code and (3) is flexible enough to support different abstraction levels and give the programmer control over details of the preconditioner. We use the composable abstractions of the Firedrake/PyOP2 package to demonstrate the efficiency of this approach for the solution of strongly anisotropic PDEs in atmospheric modelling. The weak formulation of the PDE is expressed in Unified Form Language (UFL) and the lower PyOP2 abstraction layer allows the manual design of computational kernels for a bespoke geometric multigrid preconditioner. We compare the performance of this preconditioner to a single-level method and hypre's BoomerAMG algorithm. The Firedrake/PyOP2 code is inherently parallel and we present a detailed performance analysis for a single node (24 cores) on the ARCHER supercomputer. Our implementation utilises a significant fraction of the available memory bandwidth and shows very good weak scaling on up to 6,144 compute cores.

Efficient multigrid preconditioners for atmospheric flow simulations at high aspect ratio

SummaryMany problems in fluid modelling require the efficient solution of highly anisotropic elliptic partial differential equations (PDEs) in ‘flat’ domains. For example, in numerical weather and climate prediction, an elliptic PDE for the pressure correction has to be solved at every time step in a thin spherical shell representing the global atmosphere. This elliptic solve can be one of the computationally most demanding components in semi‐implicit semi‐Lagrangian time stepping methods, which are very popular as they allow for larger model time steps and better overall performance. With increasing model resolution, algorithmically efficient and scalable algorithms are essential to run the code under tight operational time constraints. We discuss the theory and practical application of bespoke geometric multigrid preconditioners for equations of this type. The algorithms deal with the strong anisotropy in the vertical direction by using the tensor‐product approach originally analysed by Börm and Hiptmair [Numer. Algorithms, 26/3 (2001), pp. 219–234]. We extend the analysis to three dimensions under slightly weakened assumptions and numerically demonstrate its efficiency for the solution of the elliptic PDE for the global pressure correction in atmospheric forecast models. For this, we compare the performance of different multigrid preconditioners on a tensor‐product grid with a semi‐structured and quasi‐uniform horizontal mesh and a one‐dimensional vertical grid. The code is implemented in the Distributed and Unified Numerics Environment, which provides an easy‐to‐use and scalable environment for algorithms operating on tensor‐product grids. Parallel scalability of our solvers on up to 20 480 cores is demonstrated on the HECToR supercomputer. Copyright © 2015 John Wiley & Sons, Ltd.