Robust multigrid for high-order discontinuous Galerkin methods: A fast Poisson solver suitable for high-aspect ratio Cartesian grids

Highlights

• P-multigrid for discontinuous formulations of Poisson problems on Cartesian grids.

• Robust with respect to the mesh size and polynomial order, at least up to $P=32$.

• Convergence rates between 0.02 and 0.003, improving with increasing P.

• Computational complexity of PN for N unknowns, linear complexity in runtime.

• Excellent performance up to aspect ratios of 1:16.

Abstract

We present a polynomial multigrid method for nodal interior penalty and local discontinuous Galerkin formulations of the Poisson equation on Cartesian grids. For smoothing we propose two classes of overlapping Schwarz methods. The first class comprises element-centered and the second face-centered methods. Within both classes we identify methods that achieve superior convergence rates, prove robust with respect to the mesh spacing and the polynomial order, at least up to $P=32$. Consequent structure exploitation yields a computational complexity of $O(PN)$, where N is the number of unknowns. Further we demonstrate the suitability of the face-centered method for element aspect ratios up to 32.

Introduction

High-order discretization methods are exciting because of their promise to deliver higher accuracy at lower cost than first and second order methods. Much confidence has been put in the discontinuous Galerkin (DG) method because it combines multiple desirable properties of finite element and finite volume methods, including geometric flexibility, variable approximation order, straightforward adaptivity and suitability for conservation laws [6], [19]. Traditionally, DG methods have been used in the numerical solution of hyperbolic and convection-dominated problems. Nevertheless, the need for implicit diffusion schemes and application to other problem classes, such as incompressible flow and elasticity, led to a growing interest in DG methods and related solution techniques for elliptic equations [2], [30].

The most efficient elliptic solvers are based on multigrid (MG) techniques and can be classified into polynomial or p-MG [18], [11], [17], geometric or h-MG [15], [22], [10], [23], [25], [31] or, combining both concepts, hp-MG [35], [1], and algebraic MG [29], [28], [3], [32]. Apart from their different coarsening strategy, polynomial and geometric multigrid are closely related to each other and can be applied with the same smoothing methods. Early work on p-MG goes back to Helenbrook and coworkers [18], [17] who explored various smoothers for DG formulations of the Poisson equation. For isotropic grids they identified block Gauss–Seidel as the best choice, whereas more expensive line smoothing proved necessary on high-aspect-ratio grids. At about the same time, Gopalakrishnan and Kanschat [15] developed h-MG preconditioners for Poisson and convection–diffusion problems, which also use element-based block-Gauss–Seidel methods for smoothing. Kanschat [22], [23] extended this approach to locally refined Cartesian grids in two and three space dimensions. In both cases, polynomial and geometric multigrid, block-Gauss–Seidel smoothing yields acceptable convergence rates for low to moderate polynomial degrees, e.g. $\rho \approx 0.5$ with one pre-smoothing for $P=4$. However, the convergence degrades with increasing P, which renders the approach unfeasible for higher polynomial degrees.

Several researchers proposed algebraic multigrid methods for various DG formulations of elliptic equations. Olson and Schroder [28] presented a preconditioned conjugate gradient (PCG) method based on smoothed aggregation. Using block relaxation combined with energy-minimizing prolongation it attains mesh independent convergence rates corresponding to residual reductions of about 0.3 per smoothing step for $P\le 6$, but degrades with increasing approximation order. Bastian et al. [3] proposed a non-smoothed aggregation approach. For smoothing they use block relaxations which operate on extended aggregates and, hence, can be regarded as overlapping Schwarz methods. This approach yields by far the most efficient DG-MG method reported until now. The method proved robust with respect to the polynomial order, up to at least $P=6$, though the iteration count exhibits a logarithmic dependence on the mesh spacing. For the DG method of Oden, Babuška and Baumann it achieved convergence rates of $\rho \approx 0.04$ with one pre- and post-smoothing, which corresponds to a residual reduction by a factor of 25 in one step. The approach was also shown to work with symmetric and non-symmetric interior penalty methods, although it required nearly twice as many iterations in the latter case. A possible drawback is the rise of cost with increasing polynomial order. The authors did not specify the complexity of their algorithm, however the solver runtimes indicate that the cost per unknown grows as $P^3$ in 2D and $P^5$ in 3D.

To the best of our knowledge, none of the proposed MG methods is robust with respect to both, the polynomial order and the mesh spacing. Computational complexity and robustness against high aspect ratios are further issues that need to be considered to strengthen the competitiveness of DG methods for elliptic equations. As a step into that direction we present a new p-multigrid method for interior penalty and local discontinuous Galerkin discretizations of the Poisson equation on Cartesian grids. Our approach is motivated and strongly influenced by previous work dedicated to the continuous spectral element method [26], [16], [20], [33]. We propose two classes of multiplicative and weighted additive Schwarz methods, which use an adjustable overlap depending on the polynomial level. The first class comprises element-centered and the second face-centered methods. Within both classes we identify methods that achieve superior convergence rates, prove robust with respect to the mesh spacing and the polynomial order and reach a computational complexity of $O(PN)$, where N is the number of unknowns. Further we demonstrate the suitability of the face-centered method for high element aspect ratios.

The paper is organized as follows: In the next section we derive a unified nodal DG formulation of the Poisson problem comprising the symmetric interior penalty method and the local discontinuous Galerkin method. Then we describe the solution methods, i.e. Schwarz, multigrid, and inexact PCG, in Section 3. Section 4 presents the numerical experiments and Section 5 concludes the paper.