The hybridizable discontinuous Galerkin (HDG) method was introduced in [Cockburn:2009]
with the purpose of reducing the computational cost of discontinuous Galerkin (DG) methods while retaining the conservation and stability properties of DG methods. This is achieved by introducing facet variables and eliminating local (element-wise) degrees-of-freedom. This static condensation can significantly reduce the size of the global problem. Indeed, it was shown in[Kirby:2012, Yakovlev:2016] that the HDG method either outperforms or demonstrates comparable performance when compared to the CG method. This is in part also due to the local postprocessing which allows one to obtain a superconverged solution. However, they mention that these results hold true only when a direct solver is used; when an iterative solver is used the HDG method falls behind performance-wise.
The literature on iterative solvers and preconditioners for CG discretizations is vast. In contrast, there are only few studies on solvers for HDG and hybridized discretizations. We mention, for example, [Gopalakrishnan:2009] which presents a convergence analysis of multigrid for a hybridized Raviart–Thomas discretization, [Cockburn:2014]
which analyzes an auxiliary space multigrid method for HDG discretizations of elliptic partial differential equations,[Fabien:2019] which considers parallel geometric multigrid for HDG methods, [Wildey:2019] which presents a unified geometric multigrid method for hybridized finite element methods, and [Dobrev:2019] which considers the solution of hybridized systems by algebraic multigrid. Furthermore, a performance comparison of a variation of the multigrid method proposed in [Cockburn:2014] applied to CG, HDG and DG discretizations for the Poisson problem is conducted in [Kronbichler:2018]. They show that high-order continuous finite elements give the best time to solution for smooth solutions followed by matrix-free solvers for DG and HDG (using algebraic multigrid).
An alternative to HDG methods is the embedded discontinuous Galerkin (EDG) method which was introduced and analyzed in [Cockburn:2009b, Guzey:2007, Wells:2011]. The difference between an EDG and HDG method is that the facet variables in an EDG method are continuous between facets; in the HDG method they are discontinuous between facets. For the EDG method this means that after static condensation it has the same number of global degrees-of-freedom (DOFs) as a continuous Galerkin (CG) method and less DOFs than an HDG method on a given mesh. The EDG method, however, does not have the superconvergent properties of the postprocessed HDG solution and has therefore been studied less in the literature. However, we will see that the algebraic structure of the linear system resulting from an EDG method is better suited to fast iterative solvers than the linear system resulting from an HDG discretization. Indeed, this paper is motivated by the observation that multigrid methods applied to EDG discretizations of the Laplacian outperform multigrid methods applied to HDG discretizations of the Laplacian. This was observed, for example, in the context of the Stokes problem in [Rhebergen:2020] using the block-preconditioners developed in [Rhebergen:2018b].
The first goal of this paper is to present a geometric multigrid method for the hybridized and embedded discontinuous Galerkin discretizations of the Laplacian. The challenge in designing efficient multigrid methods for these discretizations is twofold. Firstly, since the facet function spaces do not form a nested hierarchy on refined grids, the design of intergrid transfer operators is not trivial. We therefore use the recently introduced Dirichlet-to-Neumann (DtN) maps proposed in [Wildey:2019] for hybridized finite element methods. The second challenge is the design of an efficient relaxation scheme. The smoothers used in multigrid methods applied to discontinuous Galerkin (DG) discretizations of an elliptic PDE are usually the classical (block) Jacobi and (block) Gauss–Seidel smoother (see for example [Gopalakrishnan:2003, Hemker:2003, Hemker:2004, Raalte:2005]). We, however, use additive Vanka-type relaxation [farrell2019local, Frommer:1999Vanka, Saad2003book] since Vanka-type relaxation is more suitable for parallel computing on general meshes than (block) Gauss–Seidel and, as we will show, results in a multigrid method that requires less iterations than when using Jacobi relaxation. Using Vanka-type relaxation requires the definition of Vanka-patches. We will study two different types of patches, namely, element-wise and vertex-wise patches. We remark that Vanka-type relaxation has been studied also in the context of discontinuous Galerkin methods for the Stokes problem in [Adler:2017].
The second and main goal of this paper is to present a two-dimensional local Fourier analysis (LFA) of the geometric two-grid error propagation operator of HDG and EDG discretizations of the Laplacian. Local Fourier analysis is used to predict the efficiency of multigrid methods. We will show that the performance of multigrid applied to an embedded discontinuous Galerkin discretization is similarly efficient and scalable to when multigrid is applied to a continuous Galerkin discretization. We furthermore show that multigrid applied to EDG outperforms multigrid applied to HDG, confirming what was previously observed but not explained in [Rhebergen:2020]. We remark that local Fourier analysis has previously been used to study the performance of multigrid applied to discontinuous Galerkin discretizations (see for example [Fidkowski:2005, Hemker:2003, Hemker:2004, Vegt:2012a, Vegt:2012b, Raalte:2005]). One-dimensional LFA has been used also to study a multilevel method for the HDG discretization of the Helmholtz equation [Chen:2014]. However, to the best of our knowledge, a two-dimensional LFA has not been applied in the context of hybridizable and embedded discontinuous Galerkin discretizations, and the performance of geometric multigrid for these types of methods has not been analyzed before.
The remainder of this work is organized as follows. In Section 2 we discuss the hybridized and embedded discontinuous Galerkin discretization of the Poisson problem. We present geometric multigrid with additive Vanka relaxation for the hybridized and embedded trace system in Section 3. A local Fourier analysis of the corresponding two-grid method is presented in Section 4. Our theory is applied and verified by numerical examples in Section 5 and we draw conclusions in Section 6.
2 The EDG and HDG methods
In this section we present EDG and HDG methods for the Poisson problem: equationparentequation
where is a bounded polygonal domain with boundary , is a given source term, and is the unknown.
2.1 The discretization
We will discretize Eq. 1 by an EDG and an HDG method. For this, denote by a tesselation of into non-overlapping quadrilateral elements . We will denote the diameter of an element by and the maximum diameter over all elements by . The boundary of an element is denoted by
and the outward unit normal vector onis denoted by . An interior face is shared by two adjacent elements and while a boundary face is a part of that lies on . We denote the set of all faces by and the union of all faces by .
Denote by and
the spaces of tensor product polynomials of degreeon, respectively, element and face . We consider the following discontinuous finite element function spaces:
For the EDG and HDG methods we then define
We note that the HDG method uses discontinuous facet function spaces and that the EDG method uses continuous facet function spaces.
For notational purposes we denote function pairs in by . For functions we write and define . Similarly, for functions where we write and .
The interior penalty EDG and HDG methods are given by [Cockburn:2009, Cockburn:2009b, Wells:2011]: find such that
Here is a penalty parameter that needs to be chosen sufficiently large [Wells:2011].
2.2 Static condensation
A feature of the EDG and HDG methods is that local (element) degrees-of-freedom can be eliminated from the discretization. For higher-order accurate discretizations this static condensation can significantly reduce the size of the problem. To obtain the reduced problem we define the function such that its restriction to the element satisfies: given and ,
for all . If satisfies Eq. 5, then where and . Furthermore, satisfies [Cockburn:2009, Rhebergen:2018b]:
We remark that Eq. 8 is the EDG or HDG method after eliminating the element degrees-of-freedom.
It will be useful to consider also the matrix representation of the EDG and HDG methods. For this, let be the vector of the discrete solution with respect to the basis for and let be the vector of the discrete solution with respect to the basis for . We can write Eq. 5 as
where , , are the matrices obtained from , , and , respectively. Since is a block diagonal matrix it is cheap to compute its inverse. Then, using we eliminate from Eq. 10 and find:
3 Geometric multigrid method
The geometric multigrid algorithm consists of: (1) applying pre-relaxation on the fine grid; (2) a coarse-grid correction step in which the residual is restricted to a coarse grid, a coarse-grid problem is solved (either exactly or by applying multigrid recursively), interpolating the resulting solution as an error correction to the fine grid approximation; and (3) applying post-relaxation.
In this section we present the different operators in a geometric multigrid method for the solution of the trace system Eq. 11. To set up notation, let be a finite sequence of increasingly coarser meshes with . For we denote the restriction operator by , the prolongation operator by , and the coarse-grid operator by .
3.1 Relaxation scheme
Many different relaxation methods may be used in multigrid algorithms. In our analysis we consider additive Vanka type relaxation (block Jacobi relaxation defined by Vanka patches) and compare its performance to the classical relaxation iterations of pointwise Jacobi and Gauss–Seidel. In this section we introduce additive Vanka relaxation relaxation [Frommer:1999Vanka, Saad2003book] following the description in [farrell2019local].
Let denote the set of DOFs of and let , , be subsets of unknowns with . Let be the restriction operator mapping from vectors over the set of all unknowns, , to vectors whose unknowns consist of the DOFs in . Then is the restriction of to the -th block of DOFs. Moreover, let for be a diagonal weight matrix for each block , where is the dimension of . Then, for a given approximation , we solve in each Vanka block the linear system
and update according to
where is a tunable parameter. For the rest of this paper it is useful to note that the error-propagation operator of the additive Vanka relaxation scheme is given by
Depending on the discretization method (EDG, HDG, and CG), the sets are chosen differently. However, for all discretization methods we consider two classes of determining , namely, via vertex-wise patches and via element-wise patches. Vertex-wise patches consist of the DOFs on the vertex , the DOFs on the interior of all edges that share vertex , and any DOFs on the interiors of the elements that contains vertex . On element-wise patches consists of all DOFs on the element and its boundary. As an example we plot the different patches for CG for , EDG for , and HDG for in Fig. 1. Note that these two Vanka-type patches are applicable also on unstructured meshes.
Given the set we next describe the weight matrix . Here is the reciprocal of the number of patches that contain DOF . For example, consider the case of Continuous Galerkin with and a vertex-wise Vanka block (see Fig. (a)a). The vertex DOF is not shared by other patches, therefore its weight is 1. The edge degrees-of-freedom are shared by two patches, therefore their weight is 1/2. The DOFs on the interior of an element are shared by four patches, therefore their weight is 1/4. For this example we furthermore note that is a matrix and, if is an matrix, is a matrix.
We end this section by noting that the size of increases with increasing degree of the polynomial approximation used in the discretization. To apply the additive Vanka smoother Eq. 14 it is necessary to compute the inverse of . To reduce the cost of computing we therefore also consider a lower-triangular approximation to . We refer to this as Lower-Triangular-Vertex-Wise and Lower-Triangular-Element-Wise Vanka relaxation.
3.2 Grid-transfer operators and the coarse-grid operator
While standard finite-element interpolation as in [Sampath2010FEM] can be used for a Continuous Galerkin finite element method, this approach is not possible for hybridized methods since with . We therefore use the Dirichlet-to-Neumann (DtN) interpolation from [Wildey:2019] which we describe next for the HDG and EDG methods.
In what follows, we assume . First we decompose the set of all fine-level faces in as , with the set of all faces in that are in the coarse-level face set (they form the boundaries of the coarse elements), and with the set of all faces in that are not in (they lie in the interior of the coarse elements). The idea is then to split the set of degrees-of-freedom of into two groups , where is the set of DOFs located on the edges in and is the set of DOFs located on the edges in . We illustrate this for the HDG method with on a rectangular mesh in Fig. 2. Then denoting by the part of corresponding to , it is clear that .
We adapt [Sampath2010FEM] to define the first part of the prolongation operator mapping, from to . In particular, we define , with , as
For the multigrid method, however, we require a prolongation operator mapping from to the whole of , which we discuss next.
We can then define the DtN prolongation operator as
We remark that this operator can be computed locally. We furthermore remark that it was shown in [Wildey:2019] that this DtN prolongation operator preserves energy when transferring information between two levels, i.e., for ,
Given the prolongation operator we then define the restriction operator as and use the Galerkin approximation of as our coarse-grid operator, i.e., .
4 A local Fourier analysis framework for HDG, EDG, and CG
Let denote a fine mesh of and denote a coarse mesh of such that . From now on, all meshes we consider will be Cartesian. The two-grid error-propagation operator is given by
where , are the number of pre- and post relaxation sweeps, respectively, and denotes the identity operator. To analyze the two-grid error-propagation operator, and hence to obtain a measure of the efficiency of the two-grid method applied to HDG, EDG, and CG discretizations of the Poisson problem, we use Local Fourier Analysis (LFA) [Trottenberg:book, Wienands:2005].
LFA was introduced in [Brandt1977LFA] to study the convergence behavior of multigrid methods for boundary value problems. Assume is a discrete operator obtained by discretizing a PDE on an infinite two dimensional domain.
can be thought of as a matrix of infinite size, but we represent it by operators that operate on the DOFs near a generic grid point and that are specified by two-dimensional stencils that we assume have constant stencil coefficients. The eigenfunctions ofcan be expressed by discrete Fourier modes, resulting in a representation of by an matrix, where and is small and depends on the discretization. This is called the symbol of . Specifically, for a scalar PDE, when there is a single degree of freedom located on the mesh of a cartesian grid, then and the symbol is a scalar. However, when there are different degrees of freedom that may be located at different locations on the grid (e.g. nodes, edges, centers, etc.), then equals the number of different degrees of freedom, and the symbol is matrix-valued.
Fundamental to LFA is that the properties of can be described by the small matrix . For example, the efficiency of multigrid on a finite grid is measured by the spectral radius of the two-grid error-propagation matrix Eq. 17. However, since is typically very large, we will find instead the spectral radius of the symbol of the two-grid error-propagation operator corresponding to the extension of
to the infinite domain. In contrast to node-based discretization problems with one DOF per node, we consider here discretizations with multiple DOFs per node, edge, and element. The key to LFA is the identification of the eigenspace ofand the symbol of , and finding the LFA representation of defined by Eq. 17. This will be done in this section. The methodology used in our LFA analysis for high-order hybridized and embedded discontinuous Galerkin discretizations of the Laplacian is similar to the methods used in the recent papers [Boonen2008curl, de2020convergence, farrell2019local, He:2019, hetwo, Maclachlan:2011, Rodrigo:2015]. The theory for the invariant space of edge-based operators was first given in [Boonen2008curl], for the curl-curl equations. While describing the general principles of the approach in the next sections, we extend the theory of [Boonen2008curl] to the case where we have multiple degrees of freedom located on multiple different grid locations, including vertical and horizontal edges, nodes and cell centers of the grid.
4.1 Infinite grid
For our analysis we follow a similar approach as [Boonen2008curl] by first defining an appropriate infinite grid and subgrids. To define these subgrids, we first lump all the -type DOFs on a horizontal edge and -type DOFs on a vertical edge to the midpoint of that edge. All -type DOFs will be lumped to the center of an element and -type DOFs are located at a node. We then consider the following two-dimensional infinite uniform grid , with subgrids
We will refer to as -type points on the grid , associated with -type DOFs. Furthermore, , , and will be referred to as -, -, and -type points on the grid . These are associated, respectively, with -, -, and -type DOFs. We remark that it is possible that there is more than one DOF at a particular location. We will use subscripts to distinguish between them. For example, and are two -type DOFs on a horizontal edge. The coarse grids are defined similarly.
The CG method for consists of -, -, -, and -type DOFs and therefore requires all four subgrids , . We refer to [hetwo] for the case . For CG only has -type DOFs and therefore requires only the grid. EDG is identical to CG for . For the EDG method for the -type DOFs have been eliminated by static condensation and therefore requires only the subgrids , . For the HDG method only has - and -type DOFs and so requires the subgrids and .
4.2 Partitioning of the discrete operator
Let be the EDG/HDG matrix given by Eq. 11 or the matrix obtained from a continuous Galerkin discretization of the Laplacian defined on the mesh . We will treat as an operator on the infinite mesh . To take into account that the DOFs on , , , and are different, we partition the operator according to the different groups of DOFs:
When thinking about as an infinite matrix, this corresponds to ordering the grid points in in the order of the subgrids , , and, for example, operator is a mapping from a grid function on -type points to a grid function on -type points.
Note furthermore that the EDG method does not include the -type DOFs and the HDG method does not include the - and -type DOFs. As an example, consider the HDG method. Then we write
If, furthermore, in HDG (so that there are two DOFs per edge), then we can again partition the submatrix , into a -block matrix, which is given by
where, for example, the operator maps from the second degree of freedom on the -edges to the first degree of freedom on the -edges.
In the general case we use the short-hand notation
where is the number of DOFs on a single edge (if or ), in a single node (if ), or on a single element (if ). Let be the total number of DOFs per element.
Let be a grid function on , and by we mean the grid function restricted to grid points . Consider operator from the degree of freedom on grid to the degree of freedom on grid . The action of linear operator on grid function is given by
where the (constant) coefficients define the stencil representation of as
Here we assume that only a finite number of coefficients are nonzero, i.e., is a finite set of offset vectors such that . Note that our notation in Eq. 23 emphasizes that the grid function is defined on grid , and its function values are linear combinations of values of grid function restricted to grid , as expressed in the right-hand side of Eq. 23.
We now describe the four different possible stencil types for . For this, let with and denote by a grid point in on which the operator acts. Furthermore, let be a fixed value in .
Stencil type 1. This considers the case where and in Eq. 24 share the same grid locations, i.e., . Let be a finite set such that for . Then
where depends on the discretization. Note that is the value of the stencil at .
Stencil type 2. Next, let be a finite set such that with and . Then
Stencil type 3. Now let be a finite set such that with and . Then
Stencil type 4. Finally, let be a finite set such that with and . Then
4.3 Properties of the symbol of
In this section we will determine the symbol of . Note that is a block operator with blocks that are each characterized by a stencil as defined in Eqs. 23 and 24. We will follow a similar approach as presented in [Boonen2008curl] to account for acting on different groups of DOFs, see Eq. Eq. 19. Our aim is to characterize the eigenfunctions of in terms of the Fourier modes
where . Each of the operators is a block Toeplitz operator with Toeplitz blocks, and as such their eigenfunctions are given by Eq. 25, but we need to determine how these eigenfunctions combine to make up invariant subspaces of , taking into account the different degrees of freedom and the different grid locations. To account for this we redefine the grid.
Let and let there be DOFs, denoted by , located at location . We define copies of , i.e., . We then define the ‘extended’ grid in which we order the grid points as follows: all grid points in followed by the grid points in , then the copies of the grid points in , the copies of the grid points in , and the copies of the grid points in (note that this ordering is consistent with the ordering of in Eq. 22). We denote the definition of the extended grid using the symbol:
We note that is defined similarly.
For example, for HDG with we write
Consider now the application of the discrete operator acting on a function defined on . Considering the restriction of the grid function to the grid, i.e., evaluating in , we obtain
Note that, in Eq. 29, the grid function is evaluated on the grid, with and .
In particular, if we take to be the Fourier mode from Eq. 25 defined on all grid points , we obtain, for the grid function evaluated on the grid, i.e., in ,