Matrix-free multigrid block-preconditioners for higher order Discontinuous Galerkin discretisations

Efficient and suitably preconditioned iterative solvers for elliptic partial differential equations (PDEs) of the convection-diffusion type are used in all fields of science and engineering. To achieve optimal performance, solvers have to exhibit high arithmetic intensity and need to exploit every form of parallelism available in modern manycore CPUs. The computationally most expensive components of the solver are the repeated applications of the linear operator and the preconditioner. For discretisations based on higher-order Discontinuous Galerkin methods, sum-factorisation results in a dramatic reduction of the computational complexity of the operator application while, at the same time, the matrix-free implementation can run at a significant fraction of the theoretical peak floating point performance. Multigrid methods for high order methods often rely on block-smoothers to reduce high-frequency error components within one grid cell. Traditionally, this requires the assembly and expensive dense matrix solve in each grid cell, which counteracts any improvements achieved in the fast matrix-free operator application. To overcome this issue, we present a new matrix-free implementation of block-smoothers. Inverting the block matrices iteratively avoids storage and factorisation of the matrix and makes it is possible to harness the full power of the CPU. We implemented a hybrid multigrid algorithm with matrix-free block-smoothers in the high order DG space combined with a low order coarse grid correction using algebraic multigrid where only low order components are explicitly assembled. The effectiveness of this approach is demonstrated by solving a set of representative elliptic PDEs of increasing complexity, including a convection dominated problem and the stationary SPE10 benchmark.Comment: 28 pages, 10 figures, 10 tables; accepted for publication in Journal of Computational Physic

Multigrid preconditioners for the hybridised discontinuous Galerkin discretisation of the shallow water equations

17 USC 105 interim-entered record; under review.The article of record as published may be found at https://doi.org/10.1016/j.jcp.2020.109948Numerical climate- and weather-prediction models require the fast solution of the equations of fluid dynamics. Discontinuous Galerkin (DG) discretisations have several advantageous properties. They can be used for arbitrary domains and support a structured data layout, which is particularly important on modern chip architectures. For smooth solutions, higher order approximations can be particularly efficient since errors decrease exponentially in the polynomial degree. Due to the wide separation of timescales in atmospheric dynamics, semi-implicit time integrators are highly efficient, since the implicit treatment of fast waves avoids tight constraints on the time step size, and can therefore improve overall efficiency. However, if implicit-explicit (IMEX) integrators are used, a large linear system of equations has to be solved in every time step. A particular problem for DG discretisations of velocity-pressure systems is that the normal Schur-complement reduction to an elliptic system for the pressure is not possible since the numerical fluxes introduce artificial diffusion terms. For the shallow water equations, which form an important model system, hybridised DG methods have been shown to overcome this issue. However, no attention has been paid to the efficient solution of the resulting linear system of equations. In this paper we address this issue and show that the elliptic system for the flux unknowns can be solved efficiently by using a non-nested multigrid algorithm. The method is implemented in the Firedrake library and we demonstrate the excellent performance of the algorithm both for an idealised stationary flow problem in a flat domain and for nonstationary setups in spherical geometry from the well-known testsuite in Williamson et al. (1992) [23]. In the latter case the performance of our bespoke multigrid preconditioner (although itself not highly optimised) is comparable to that of a highly optimised direct solver.EPSRCEPSRCEP/L015684/1UK-Fluids network (EPSRC grant EP/N032861/1

Hybridised multigrid preconditioners for a compatible finite element dynamical core

Compatible finite element discretisations for the atmospheric equations of motion have recently attracted considerable interest. Semi-implicit timestepping methods require the repeated solution of a large saddle-point system of linear equations. Preconditioning this system is challenging since the velocity mass matrix is non-diagonal, leading to a dense Schur complement. Hybridisable discretisations overcome this issue: weakly enforcing continuity of the velocity field with Lagrange multipliers leads to a sparse system of equations, which has a similar structure to the pressure Schur complement in traditional approaches. We describe how the hybridised sparse system can be preconditioned with a non-nested two-level preconditioner. To solve the coarse system, we use the multigrid pressure solver that is employed in the approximate Schur complement method previously proposed by the some of the authors. Our approach significantly reduces the number of solver iterations. The method shows excellent performance and scales to large numbers of cores in the Met Office next-generation climate- and weather prediction model LFRic.Comment: 24 pages, 13 figures, 5 tables; accepted for publication in Quarterly Journal of the Royal Meteorological Societ

Multigrid Preconditioners for the Discontinuous Galerkin Spectral Element Method : Construction and Analysis

Discontinuous Galerkin (DG) methods offer a great potential for simulations of turbulent and wall bounded flows with complex geometries since these high-order schemes offer a great potential in handling eddies. Recently, space-time DG methods have become more popular. These discretizations result in implicit schemes of high order in both spatial and temporal directions. In particular, we consider a specific DG variant, the DG Spectral Element Method (DG-SEM), which is suitable to construct entropy stable solvers for conservation laws. Since the size of the corresponding nonlinear systems is dependent on the order of the method in all dimensions, the problem arises to efficiently solve these huge nonlinear systems with regards to CPU time as well as memory consumption.Currently, there is a lack of good solvers for three-dimensional DG applications, which is one of the major obstacles why these high order methods are not used in e.g. industry. We suggest to use Jacobian-free Newton- Krylov (JFNK) solvers, which are advantageous in memory minimization. In order to improve the convergence speed of these solvers, an efficient preconditioner needs to be constructed for the Krylov sub-solver. However, if the preconditioner requires the storage of the DG system Jacobian, the favorable memory consumption of the JFNK approach is obsolete.We therefore present a multigrid based preconditioner for the Krylov sub-method which retains the low mem- ory consumption, i.e. a Jacobian-free preconditioner. To achieve this, we make use of an auxiliary first order finite volume replacement operator. With this idea, the original DG mesh is refined but can still be implemented algebraically. As smoother, we consider the pseudo time iteration W3 with a symmetric Gauss-Seidel type approx- imation of the Jacobian. Numerical results are presented demonstrating the potential of the new approach.In order to analyze multigrid preconditioners, a common tool is the Local Fourier Analysis (LFA). For a space- time model problem we present this analysis and its benefits for calculating smoothing and two-grid convergence factors, which give more insight into the efficiency of the multigrid method

Composable code generation for high order, compatible finite element methods

It has been widely recognised in the HPC communities across the world, that exploiting modern computer architectures, including exascale machines, to a full extent requires software commu- nities to adapt their algorithms. Computational methods with a high ratio of floating point op- erations to bandwidth are favorable. For solving partial differential equations, which can model many physical problems, high order finite element methods can calculate approximations with a high efficiency when a good solver is employed. Matrix-free algorithms solve the corresponding equations with a high arithmetic intensity. Vectorisation speeds up the operations by calculating one instruction on multiple data elements. Another recent development for solving partial differential are compatible (mimetic) finite ele- ment methods. In particular with application to geophysical flows, compatible discretisations ex- hibit desired numerical properties required for accurate approximations. Among others, this has been recognised by the UK Met office and their new dynamical core for weather and climate fore- casting is built on a compatible discretisation. Hybridisation has been proven to be an efficient solver for the corresponding equation systems, because it removes some inter-elemental coupling and localises expensive operations. This thesis combines the recent advances on vectorised, matrix-free, high order finite element methods in the HPC community on the one hand and hybridised, compatible discretisations in the geophysical community on the other. In previous work, a code generation framework has been developed to support the localised linear algebra required for hybridisation. First, the framework is adapted to support vectorisation and further, extended so that the equations can be solved fully matrix-free. Promising performance results are completing the thesis.Open Acces

Higher-order compatible finite element schemes for the nonlinear rotating shallow water equations on the sphere

We describe a compatible finite element discretisation for the shallow water equations on the rotating sphere, concentrating on integrating consistent upwind stabilisation into the framework. Although the prognostic variables are velocity and layer depth, the discretisation has a diagnostic potential vorticity that satisfies a stable upwinded advection equation through a Taylor-Galerkin scheme; this provides a mechanism for dissipating enstrophy at the gridscale whilst retaining optimal order consistency. We also use upwind discontinuous Galerkin schemes for the transport of layer depth. These transport schemes are incorporated into a semi-implicit formulation that is facilitated by a hybridisation method for solving the resulting mixed Helmholtz equation. We illustrate our discretisation with some standard rotating sphere test problems.Comment: accepted versio

Efficient Solvers for Space-Time Discontinuous Galerkin Spectral Element Methods

In this thesis we study efficient solvers for space-time discontinuous Galerkin spectral element methods (DG-SEM). These discretizations result in fully implicit schemes of variable order in both spatial and temporal directions. The popularity of space-time DG methods has increased in recent years and entropy stable space-time DG-SEM have been constructed for conservation laws, making them interesting for these applications. The size of the nonlinear system resulting from differential equations discretized with space-time DG-SEM is dependent on the order of the method, and the corresponding Jacobian is of block form with dense blocks. Thus, the problem arises to efficiently solve these huge nonlinear systems with regards to CPU time as well as memory consumption. The lack of good solvers for three-dimensional DG applications has been identified as one of the major obstacles before high order methods can be adapted for industrial applications.It has been proven that DG-SEM in time and Lobatto IIIC Runge-Kutta methods are equivalent, in that both methods lead to the same discrete solution. This allows to implement space-time DG-SEM in two ways: Either as a full space-time system or by decoupling the temporal elements and using implicit time-stepping with Lobatto IIIC methods. We compare theoretical properties and discuss practical aspects of the respective implementations.When considering the full space-time system, multigrid can be used as solver. We analyze this solver with the local Fourier analysis, which gives more insight into the efficiency of the space-time multigrid method. The other option is to decouple the temporal elements and use implicit Runge-Kutta time-stepping methods. We suggest to use Jacobian-free Newton-Krylov (JFNK) solvers since they are advantageous memory-wise. An efficient preconditioner for the Krylov sub-solver is needed to improve the convergence speed. However, we want to avoid constructing or storing the Jacobian, otherwise the favorable memory consumption of the JFNK approach would be obsolete. We present a preconditioner based on an auxiliary first order finite volume replacement operator. Based on the replacement operator we construct an agglomeration multigrid preconditioner with efficient smoothers using pseudo time integrators. Then only the Jacobian of the replacement operator needs to be constructed and the DG method is still Jacobian-free. Numerical experiments for hyperbolic test problems as the advection, advection-diffusion and Euler equations in several dimensions demonstrate the potential of the new approach

h-multigrid agglomeration based solution strategies for discontinuous Galerkin discretizations of incompressible flow problems

In this work we exploit agglomeration based $h$-multigrid preconditioners to speed-up the iterative solution of discontinuous Galerkin discretizations of the Stokes and Navier-Stokes equations. As a distinctive feature $h$-coarsened mesh sequences are generated by recursive agglomeration of a fine grid, admitting arbitrarily unstructured grids of complex domains, and agglomeration based discontinuous Galerkin discretizations are employed to deal with agglomerated elements of coarse levels. Both the expense of building coarse grid operators and the performance of the resulting multigrid iteration are investigated. For the sake of efficiency coarse grid operators are inherited through element-by-element $L^2$ projections, avoiding the cost of numerical integration over agglomerated elements. Specific care is devoted to the projection of viscous terms discretized by means of the BR2 dG method. We demonstrate that enforcing the correct amount of stabilization on coarse grids levels is mandatory for achieving uniform convergence with respect to the number of levels. The numerical solution of steady and unsteady, linear and non-linear problems is considered tackling challenging 2D test cases and 3D real life computations on parallel architectures. Significant execution time gains are documented.Comment: 78 pages, 7 figure