Multilevel Preconditioners for DG methods for Jump coefficient problems

# Multilevel Preconditioners for Discontinuous Galerkin Approximations of Elliptic Problems with Jump Coefficients

Blanca Ayuso De Dios Michael Holst Yunrong Zhu  and  Ludmil Zikatanov
July 12, 2019
###### Abstract.

We introduce and analyze two-level and multi-level preconditioners for a family of Interior Penalty (IP) discontinuous Galerkin (DG) discretizations of second order elliptic problems with large jumps in the diffusion coefficient. Our approach to IPDG-type methods is based on a splitting of the DG space into two components that are orthogonal in the energy inner product naturally induced by the methods. As a result, the methods and their analysis depend in a crucial way on the diffusion coefficient of the problem. The analysis of the proposed preconditioners is presented for both symmetric and non-symmetric IP schemes; dealing simultaneously with the jump in the diffusion coefficient and the non-nested character of the relevant discrete spaces presents extra difficulties in the analysis which precludes a simple extension of existing results. However, we are able to establish robustness (with respect to the diffusion coefficient) and nearly-optimality (up to a logarithmic term depending on the mesh size) for both two-level and BPX-type preconditioners. Following the analysis, we present a sequence of detailed numerical results which verify the theory and illustrate the performance of the methods. The paper includes an Appendix with a collection of proofs of several technical results required for the analysis.

###### Key words and phrases:
Multilevel preconditioner, discontinuous Galerkin methods, Crouzeix-Raviart finite elements, space decomposition

## 1. Introduction

Let be a bounded polygon (for ) or polyhedron (for ) and . We consider the following second order elliptic equation with strongly discontinuous coefficients:

 {−∇⋅(κ∇u)=f in Ω,u=0 on ∂Ω. (1.1)

The scalar function denotes the diffusion coefficient which is assumed to be piecewise constant with respect to an initial non-overlapping (open) subdomain partition of the domain , denoted with and for . Although the (polygonal or polyhedral) regions might have complicated geometry, we will always assume that there is an initial shape-regular triangulation such that is a constant for all . Problem (1.1) belongs to the class of interface or transmission problems, which are relevant to many applications such as groundwater flow, electromagnetics and semiconductor device modeling. The coefficients in these applications might have large discontinuities across the interfaces between different regions with different material properties. Finite element discretizations of (1.1) lead to linear systems with badly conditioned stiffness matrices. The condition numbers of these matrices depend not only on the mesh size, but also on the largest jump in the coefficients.

Much research has been devoted to developing efficient and robust preconditioners for conforming finite element discretizations of (1.1). Nonoverlapping domain decomposition preconditioners, such as Balancing Neumann-Neumann , FETI-DP  and Bramble-Pasciak-Schatz Preconditioners  have been shown to be robust with respect to coefficient variations and mesh size (up to a logarithmic factor), in theory and in practice, but only if special exotic coarse solvers (such as those based on discrete harmonic extensions [37, 48, 38]) are used (see also ). The construction and use of such exotic coarse spaces is avoided in other multilevel methods, such as the Bramble-Pasciak-Xu (BPX) or multigrid preconditioners, for which it has always been observed that when used with conjugate gradient (CG) iteration, result in robust and efficient algorithms with respect to jumps in the coefficients, independently of the problem dimension. However, their analysis (based on the standard CG theory) predict a deterioration in the rate of convergence with respect to both the coefficients and the mesh size, By resorting to more sophisticated CG theory (see [6, Section 13.2], ) which accounts for and exploits the particular spectral structure of the preconditioned systems111Namely, that there are a few small eigenvalues due to the jump coefficient distribution that have no influence in the (observed) overall convergence of the iteration, the authors in [58, 61] show that standard multilevel and overlapping domain decomposition methods lead to nearly optimal preconditioners for CG algorithms. (See also ). Much less attention has been devoted to nonconforming approximations. Overlapping preconditioners for the lowest order Crouzeix-Raviart approximation of (1.1) are found in [52, 51], where the analysis depends on the assumption that the coefficient is quasi-monotone.

In this article, we consider the construction and analysis of preconditioners for the Interior Penalty (IP) Discontinuous Galerkin (DG) approximation of (1.1). Based on discontinuous finite element spaces, DG methods can deal robustly with partial differential equations of almost any kind, as well as with equations whose type changes within the computational domain. They are naturally suited for multi-physics applications, and for problems with highly varying material properties, such as (1.1). The design of efficient solvers for DG discretizations has been pursued only in the last ten years; and, while classical approaches have been successfully extended to second order elliptic problems, the discontinuous nature of the underlying finite element spaces has motivated the creation of new techniques to develop solvers. Additive Schwarz methods (of overlapping and non-overlapping type) are considered and analyzed in [39, 34, 2, 3, 4, 11]. Multigrid methods are studied in [41, 20, 18, 17, 50, 29]. Two-level methods are presented in [31, 22, 23]. More general multi-level methods based on algebraic techniques are considered in [47, 46]. However, all the analysis in these works consider only the case of a smoothly or slowly varying diffusivity coefficient. For problem (1.1), only in [34, 35, 36] have the authors introduced and analyzed non-overlapping BBDC and FETI-DP domain decomposition preconditioners for a Nitsche type method where a Symmetric Interior Penalty DG discretization is used (only) on the skeleton of the subdomain partition, while a standard conforming approximation is used in the interior of the subdomains. Robustness and quasi-optimality is shown in for the Additive and Hybrid BBDC  and FETI-DP  preconditioners, even for the case of non-matching grids. As it happens for conforming discretizations, the construction and analysis of these preconditioners rely on the use of exotic coarse solvers, which might complicate the actual implementation of the method.

The goal of this article is to design, and provide a rigorous analysis of, a simple multilevel solver for the lowest order (i.e. piecewise linear discontinuous) approximation of a family of Interior Penalty (IPDG) methods. To ease the presentation, we focus on a minor variant of the classical IP methods, penalizing only the mean value of the jumps: the “weakly penalized” or IPDG-0 methods (called Type-0 in ). Our approach follows the ideas in , and it is based on a splitting of the DG space into two components that are orthogonal in the energy inner product naturally induced by the IPDG-0 methods.

Roughly speaking, the construction amounts to identifying a “low frequency” space (the Crouzeix-Raviart elements) and then defining a complementary space. However, a notable difference takes place in the DG space decomposition introduced for the Laplace equation [24, 10]. For problem (1.1), the subspaces depend on the coefficient , and this is certainly related to the splittings used in algebraic multigrid (AMG ). With the orthogonal splitting of the DG space at hand, the solution of problem (1.1) reduces to solving two sub-problems: a non-conforming approximation to (1.1), and a problem in the complementary space containing high oscillatory error components. We show the latter subproblem is easy to solve, since it is spectrally equivalent to its diagonal form, and so CG with a diagonal preconditioner is a uniform and robust solver.

For the former subproblem, following [58, 61], we develop and analyze (in the standard and asymptotic convergence regimes) a two-level method and a BPX preconditioner. Nevertheless, dealing simultaneously with the jump in the coefficient and the non-nested character of the Crouziex-Raviart (CR) spaces presents extra difficulties in the analysis which precludes a simple extension of [58, 61]. We are able to establish nearly optimal convergence and robustness (with respect to both the mesh size and the coefficient ) for the two-level method and for the BPX preconditioner (up to a logarithmic term depending on the mesh size). The resulting algorithms involve the use of a solver in the CR space that is reduced to a smoothing step followed by a conforming solver. Therefore, in particular one can argue that any of the robust and efficient solvers designed for conforming approximations of problem (1.1) could be used as a preconditioner here. Finally we mention that, although the two-level and multilevel methods we propose are based on the piecewise linear IP-0 methods, they could be used as preconditioners for the solution of the linear systems arising from high order DG methods.

### Outline of the paper

The rest of the paper is organized as follows. We introduce the IPDG-1 and IPDG-0 methods for approximating (1.1) in §2 and revise some of their properties. The space decomposition of DG finite element space is introduced in §3. Consequences of the space splitting are described in §4. The two-level and multi-level methods for the Crouzeix-Raviart approximation of (1.1) are constructed and analyzed in §5. Numerical experiments are included in §6, to verify the theory and assess the performance and robustness of the proposed preconditioners. In §7 we briefly comment on how the developed solvers and theory can be extended for the classical IPDG-1 family. The paper is completed with an Appendix where we have collected proofs of several technical results required in our analysis.

Throughout the paper we shall use the standard notation for Sobolev spaces and their norms. We will use the notation , and , whenever there exist constants independent of the mesh size and the coefficient or other parameters that , , and may depend on, and such that and , respectively. We also use the notation for .

## 2. Discontinuous Galerkin Methods

In this section, we introduce the basic notation and describe the DG methods we consider for approximating the problem (1.1).

Let be a shape-regular family of partitions of into -simplices (triangles in or tetrahedra in ). We denote by the diameter of and we set . We also assume that the decomposition is conforming in the sense that it does not contain hanging nodes and that , with being quasi-uniform initial triangulation that resolves the coefficient . We denote by the set of all edges/faces and by and the collection of all interior and boundary edges/faces, respectively. The space is the set of element-wise functions, and refers to the set of functions whose traces on the elements of are square integrable.

Following , we recall the usual DG analysis tools. Let and be two neighboring elements, and let , be their outward normal unit vectors, respectively (). Let and be the restriction of and to . We set:

 2{{ζ}} =(ζ++ζ−),[[ζ]]=ζ+n++ζ−n− on e∈Eoh, 2{{τ}} =(τ++τ−),[[τ]]=τ+⋅n++τ−⋅n− on e∈Eoh.

We also define the weighted average, , for any with :

 {{ζ}}δe=δeζ++(1−δe)ζ−,{{τ}}δe=δeτ++(1−δe)τ−,on e∈Eoh. (2.1)

For , we set

 [[ζ]]=ζn,{{τ}}={{τ}}δe=τon e∈E∂h. (2.2)

We will also use the notation

 (u,w)Th=∑T∈Th∫Tuwdx∀u,w∈L2(Ω),⟨u,w⟩Eh=∑e∈Eh∫euwds∀u,w,∈L2(Eh).

The DG approximation to the model problem (1.1) can be written as

where is the piecewise linear discontinuous finite element space, and is the bilinear form defining the method.

In this paper, we focus on a family of weighted Interior Penalty methods (see ), with special attention given to a variant (weakly penalized) of them. The bilinear form defining the classical family of weighted IP methods , here called IP()-1 methods, is given by , with

 A(v,w) =(κ∇hv,∇w)Th−⟨{{κ∇v}}βe,[[w]]⟩Eh+θ⟨[[v]],{{κ∇w}}βe⟩Eh (2.3) +⟨αh−1eκe[[v]],[[w]]⟩Eh,∀v,w∈VDGh.

where gives the SIPG()-1 methods; leads to NIPG()-1 methods; and gives the IIPG()-1 methods. Here, denotes the dimensional Lebesgue measure of .The penalty parameter is set to be a positive constant; and it has to be taken large enough to ensure coercivity of the corresponding bilinear forms when . The symmetric method was first considered in  and later in [33, Section 4] for jump coefficient problems (although there it was written using a slightly different notation and DG was only used in the skeleton of the partition). It was later extended to advection-diffusion problems in  and .

We also introduce the corresponding family of IP()-0 methods, which use the mid-point quadrature rule for computing the integrals in the last term in (2.3) above. That is, we set with

 A0(v,w) =(κ∇v,∇w)Th−⟨{{κ∇v}}βe,[[w]]⟩Eh+θ⟨[[v]],{{κ∇w}}βe⟩Eh (2.4) +⟨αh−1eκeP0e([[v]]),[[w]]⟩Eh,∀v,w∈VDGh,

where is the -projection onto the piecewise constants on . We note that this projection satisfies . In (2.3) and (2.4), for any with , the coefficient and the weight are defined as follows:

 κT=κ|T,βe=κ−κ++κ−,whereκ±=κ|T±, (2.5)

The coefficient as the harmonic mean of and :

 κe:=2κ+κ−κ++κ−. (2.6)

The weight depends on the coefficient and therefore it might vary over all interior edges/faces (of the subdomain partition resolving the coefficient ).

###### Remark 2.1.

We note that one could take as , since both are equivalent:

 min{κ+,κ−}≤κe=2κ+κ−κ++κ−≤2min{κ+,κ−}≤2κ±. (2.7)

The equivalence relations in (2.7) show that the results on spectral equivalence and uniform preconditioning given later for (2.3) with defined in (2.6) (the harmonic mean) will automatically hold for method (2.3) with . To fix the notation and simplify the presentation, we stick to definition (2.6) for .

### Weighted Residual Formulation

Following  we can rewrite the two families of IP methods in the weighted residual framework: For all ,

 A(v,w) =(−∇⋅(κ∇v),w)Th+⟨[[κ∇v]],{{w}}1−βe⟩Eoh+⟨[[v]],B1(w)⟩Eh, (2.8) A0(v,w) =(−∇⋅(κ∇v),w)Th+⟨[[κ∇v]],{{w}}1−βe⟩Eoh+⟨[[v]],P0e(B1(w))⟩Eh, (2.9)

where is defined as:

 B1(w)=θ{{κ∇w}}βe+αh−1eκe[[w]],∀e∈Eh. (2.10)

Throughout the paper both the weighted residual formulation (2.8)-(2.9) and the standard one (2.3)-(2.4) will be used interchangeably.

We now establish a result that guarantees the spectral equivalence between and .

###### Lemma 2.2.

Let be a bilinear form corresponding to a IP()-1 method (2.3) and let be the corresponding IP()-0 bilinear form as defined in (2.4). Then there exists a positive constant , depending only on the shape regularity of the mesh and the penalty parameter (but independent of the coefficient and the mesh size ) such that,

 A0(v,v)≤A(v,v)≤c0(α)A0(v,v)∀v∈VDGh. (2.11)
###### Proof.

The lower bound follows immediately from the fact that the projection is an -orthogonal projection and therefore has unit norm. The upper bound would follow if we show

 ∑e∈Ehαh−1eκe∥[[v]]∥20,e≤C(∑T∈ThκT∥∇v∥20,T+∑e∈Ehαh−1eκe∥P0e[[v]]∥20,e),

which can be proved by arguing exactly as in [10, 19, 8] and taking into account (2.7). ∎

By virtue of Lemma 2.2, it will be enough throughout the rest of the paper to focus on the design and analysis of multilevel preconditioners for the IP()-0 methods. At least in the symmetric case, the preconditioners proposed for SIPG()-0 will exhibit the same convergence (asymptotically) when applied to SIPG()-1.

### Continuity and Coercivity of IP(β)-0 methods

The family of methods (2.4) can be shown to provide an accurate and robust approximation to the solution of (1.1). We define the energy norm:

 |||v|||2DG0:=∑T∈ThκT∥∇v∥20,T+∑e∈Ehκeh−1e∥P0e([[v]])∥20,e. (2.12)

Then, is continuous and coercive in the above norm, with constants independent of the mesh size and the coefficient :

 Continuity: |A0(v,w)| ∀v,w∈VDGh, (2.13) Coercivity: A0(v,v) ∀v∈VDG0h. (2.14)

Although the proof of (2.14) and (2.13) is standard, we sketch it here for completeness. Note first that for each such that , the weighted average can be rewritten as:

 {{κ∇v}}βe =βe(κ+(∇v)+)+(1−βe)(κ−(∇v)−) =κ−κ++κ−κ+(∇v)++κ+κ++κ−κ−(∇v)− =κ+κ−κ++κ−[(∇v)++(∇v)−]=κe{{∇v}}. (2.15)

Trace inequality , inverse inequality  and (2.7) imply the following bounds

 he∥{{κ∇v}}βe∥20,e ≤Ct(κe)2(∥∇v∥20,T+∪T−+h2|∇v|21,T+∪T−) ≤2(κe)Ct(1+C2inv)(κ+∥∇v∥20,T++κ−∥∇v∥20,T−).

This inequality, combined with Cauchy-Schwarz inequality and (2.7), gives

 ∣∣⟨{{κ∇v}}βe,[[w]]⟩Eh∣∣ =∣∣ ∣∣∑e∈Eh∫eκe{{∇v}}P0e([[w]])ds∣∣ ∣∣ ≤⎛⎝∑e∈Eh1αheκe∥{{∇v}}∥20,e⎞⎠1/2⎛⎝∑e∈Ehαh−1eκe∥P0e([[w]])∥20,e⎞⎠1/2 ≤8Ct(1+C2inv)α∑T∈ThκT∥∇v∥20,T+14∑e∈Ehαh−1eκe∥P0e([[w]])∥20,e.

Now (2.13) follows from Cauchy-Schwarz inequality. The inequality (2.14) is proved by setting in (2.3) and taking into account the above estimate. We have then

 A0(v,v) =∑T∈ThκT∥∇v∥20,T+α∑e∈Ehκeh−1e∥P0e([[v]])∥20,e−(1−θ)⟨{{κ∇v}}βe,[[v]]⟩Eh ≥|||v|||2DG−|1−θ|∣∣⟨{{κ∇u}}βe,P0e([[v]])⟩Eh∣∣ ≥(1−8Ct(1+C2inv)α)∑T∈ThκT∥∇v∥20,T+4−|1−θ|4α∑e∈Ehκeh−1e∥P0e([[v]])∥20,e,

and (2.14) follows immediately by taking large enough (if ). Moreover, notice that both constants in (2.13) and (2.14) depend on the shape regularity of the mesh partition but are independent of the coefficient .

Obviously, continuity and coercivity also hold for the IP()-1 methods (2.3) if the norm (2.12) is replaced by

 |||v|||2DG:=∑T∈ThκT∥∇v∥20,T+∑e∈Ehκeh−1e∥[[v]]∥20,e. (2.16)

See  or  for a detailed proof. For both families of methods, optimal error estimates in the energy norms (2.12) and (2.16) can be shown, arguing as in . See also  for further discussion on the -error analysis of these methods.

## 3. Space decomposition of the VDGh space

In this section, we introduce a decomposition of the -space that will play a key role in the design of the solvers for the DG discretizations (2.3) and (2.4). In [10, 24], it is shown that the discontinuous piecewise linear finite element space admits the decomposition: , where denotes the standard Crouzeix-Raviart space defined as

 VCRh={v∈L2(Ω):v|T∈P1(T)∀T∈Th and P0e([[v]]⋅n)=0∀e∈Eoh}, (3.1)

and the complementary space is a space of piece-wise linear functions with average zero at the mass centers of the internal edges/faces:

 Z={z∈L2(Ω):z|T∈P1(T)∀T∈Th and P0e({{v}})=0,∀e∈Eoh}.

In , it was shown that this decomposition satisfies when , for all and . We now modify the definition of above in order to account for the presence of a coefficient in the problem (1.1). Let

 Zβ={z∈L2(Ω):z|T∈P1(T)∀T∈Th and P0e({{z}}1−βe)=0,∀e∈Eoh}, (3.2)

where the weight was defined earlier in (2.5). Note that the weight depends on the coefficient , and, as a consequence, the space is also coefficient dependent. In what follows, we shall show that is a space complementary to in and the corresponding decomposition has properties analogous to the properties of the decomposition given in  for the Poisson problem.

For any with , let be the canonical Crouzeix-Raviart basis function on , which is defined by

 φe,T|T∈P1(T),φe,T(me′)=δe,e′∀e′∈Eh(T), % and φe,T(x)=0∀x∉T,

where is the mass center of . We will denote by and the number of simplices and faces (or edges when ) respectively. We also denote by the number of boundary faces.

###### Proposition 3.1.

For any there exists a unique and a unique such that , that is

 VDGh=VCRh⊕Zβ. (3.3)
###### Proof.

For simplicity, throughout the proof we will set , , and for any with We also denote for any with . Since the mesh is made of -simplices

 dimVDGh=(d+1)nT=2nE−nBE,

and it is also obvious that form a basis for Notice that , we can therefore express any as

 u(x) = ∑e∈Eohu+(me)φ+e(x)+∑e∈Eohu−(me)φ−e(x)+∑e∈E∂hu(me)φe(x) = ∑e∈Eoh(β−u+(me)+β+u−(me))(φ+e(x)+φ−e(x)) +∑e∈Eoh(u+(me)−u−(me))(β+φ+e(x)−β−φ−e(x))+∑e∈E∂hu(me)φe(x) = ∑e∈Eoh(1|e|∫e{{u}}1−βeds)(φ+e(x)+φ−e(x)) +∑e∈Eoh(1|e|∫e[[u]]n+ds)(β+φ+e(x)−β−φ−e(x))+∑e∈E∂h(1|e|∫e[[u]]nds)φe(x) = v(x)+zβ(x).

Then for each , we set

 φCRe(x):=φ+e(x)+φ−e(x), (3.4)
 (3.5)

and for all In the definition (3.5) of , we have used for and for Finally, when with for some , we set

 ψze(x)=φe(x),∀x∈T. (3.6)

It is then straightforward to check that

 VCRh=span{φCRe}e∈Eoh,andZβ=span{ψze}e∈Eh.

Hence, for all there exist unique and defined by

 v = ∑e∈Eoh(1|e|∫e{{u}}1−βeds)φCRe(x)∈VCRh, zβ = ∑e∈Eh(1|e|∫e[[u]]n+ds)ψze(x)∈Zβ,

such that . This shows (3.3) and concludes the proof. ∎

###### Remark 3.2.

As we pointed out in the introduction, the definition of the subspace clearly depends on the coefficient , since depends on . Such dependence is often also seen in algebraic multigrid analysis, where the coarse spaces depend on the operator at hand. They are in fact explicitly constructed in this way, the aim being to increase robustness of the methods.

In the proof of Proposition 3.1 above, we have introduced the basis in both and . The canonical Crouzeix-Raviart basis functions are continuous at the mass centers of the faces . The basis in consists of piecewise functions, which are discontinuous across the faces in In fact, for any such that with , we have

 ([[z]]n+)(me′)=ze′,∀e′∈Eh.

To see this, evaluating the jump of at gives

 ([[z]]n+)(me′) = ∑e∈Ehze([[ψze]]n+)(me′)=ze′([[ψze′]]n+)(me′) = {ze′(βe′−(βe′−1))=ze′,e′∈Eoh,ze′,e′∈E∂h.

This relation will also be used later to obtain uniform diagonal preconditioners for the restrictions of and on .

###### Remark 3.3.

For mixed boundary value problems, that is, contains both Neumann boundary and Dirichlet boundary with , the definition of the basis functions on the boundary faces [see (3.6)] needs to be changed as:

 ϕCRe(x)=φe,T(x),e=∂T∩ΓN,for allx∈T,ψze(x)=φe,T(x),e=∂T∩ΓD,for allx∈T. (3.7)

Thus, in case the dimension of is increased (by adding to it functions that correspond to degrees of freedom on ) and the dimension of is decreased accordingly. Clearly things balance out correctly: the identity holds, and also the analysis carries over with very little modification.

Next lemma is a simple but key observation used in the design of efficient solvers.

###### Lemma 3.4.

Let be the bilinear form defined in (2.4). Then,

 A0(v,z)=0∀v∈VCRh,∀z∈Zβ. (3.8)

Furthermore if is symmetric (and positive definite) then the decomposition (3.3) is -orthogonal, namely, .

###### Proof.

From the weighted-residual form of given in (2.9), for all , and all we easily obtain

 A0(v,z) =(−∇⋅(κ∇v),z)Th+⟨[[κ∇v]],{{z}}1−βe⟩Eoh+⟨[[v]],P0e(B1(z))⟩Eh=0.

In the equation above, the first term is zero due to the fact that , so is linear in each , and the coefficient . Last term vanishes (independently of the choice of , or equivalently the choice of ), because , thanks to the definition (3.1) of . The second term vanishes from the definition of (since is constant on each ). Moreover, in the case when is symmetric and positive definite we have that , for all and for all . Thus, for the symmetric method , the spaces and are indeed -orthogonal. The proof is complete. ∎

## 4. Solvers for IP(β)-0 methods

In this section we show how Proposition 3.1 and Lemma 3.4 can be used in the design and analysis of uniformly convergent iterative methods for the IP()-0 methods. We follow the ideas and analysis introduced in  and point out the differences. We first consider the approximation to problem (1.1) with . To begin, let be the discrete operator defined by and let be its matrix representation in the new basis (3.4) and (3.5). We denote by , the vector representation of the unknown function and of the right hand side , respectively, in this new basis. A simple consequence of Lemma 3.4 is that the matrix (in this basis) has block lower triangular structure:

 A0=[Azz00Avz0Avv0], (4.1)

where are the matrix representation of restricted to the subspaces and , respectively, and is the matrix representation of the term that accounts for the coupling (or non-symmetry) . As remarked earlier, for SIPG-0, the stiffness matrix is block-diagonal.

Figure 4.1 gives a 2D example, with two squares and inside the domain We set the coefficients for all and for .