Analysis of Schwarz methods for a hybridizable discontinuous Galerkin discretization

Analysis of Schwarz methods for a hybridizable discontinuous Galerkin discretization

Martin J. Gander    Soheil Hajian
Abstract

Schwarz methods are attractive parallel solvers for large scale linear systems obtained when partial differential equations are discretized. For hybridizable discontinuous Galerkin (HDG) methods, this is a relatively new field of research, because HDG methods impose continuity across elements using a Robin condition, while classical Schwarz solvers use Dirichlet transmission conditions. Robin conditions are used in optimized Schwarz methods to get faster convergence compared to classical Schwarz methods, and this even without overlap, when the Robin parameter is well chosen. We present in this paper a rigorous convergence analysis of Schwarz methods for the concrete case of hybridizable interior penalty (IPH) method. We show that the penalization parameter needed for convergence of IPH leads to slow convergence of the classical additive Schwarz method, and propose a modified solver which leads to much faster convergence. Our analysis is entirely at the discrete level, and thus holds for arbitrary interfaces between two subdomains. We then generalize the method to the case of many subdomains, including cross points, and obtain a new class of preconditioners for Krylov subspace methods which exhibit better convergence properties than the classical additive Schwarz preconditioner. We illustrate our results with numerical experiments.

Key words. Additive Schwarz, optimized Schwarz, discontinuous Galerkin methods

AMS subject classifications. 65N22, 65F10, 65F08, 65N55, 65H10

1 Introduction

We consider the elliptic model problem

\hb@xt@.01(1.1)

in the weak sense where , and uniformly positive, and is assumed to be a convex polygon for simplicity. Any discretization of this problem, for example by a finite element method (FEM) or a discontinuous Galerkin (DG) method, leads to a large sparse linear system

\hb@xt@.01(1.2)

where is the vector of degrees of freedom representing an approximation of and represents the disretized differential operator. In this paper we consider a hybridizable interior penalty (IPH111We use the acronym IPH for hybridizable interior penalty because this has become the common abbreviation following its introduction in [7] as a member of the family of HDG methods.) discretization which results in a symmetric positive definite (s.p.d.) matrix . An IPH discretization seeks over a triangulation of the domain where is not necessarily continuous across elements. As common to DG methods, IPH imposes the continuity of the solution approximately through penalization techniques, i.e. penalizing jumps of across elements in the bilinear form. The penalization is controlled by a penalty parameter .

Since the matrix of IPH is s.p.d. and sparse, one can use the Conjugate Gradient (CG) method to solve the linear system (LABEL:eq:linsys). The convergence of CG slows down as the condition number grows. It is not hard to show that , where is the maximum diameter of the elements in the triangulation, see for instance [6]. Therefore preconditioning is unavoidable and domain decomposition (DD) preconditioners have been developed and studied for such discretizations, see [2, 12]. IPH as local solvers were also used to precondition classical IP discretizations [1]. One can also design a substructuring preconditioner for a -version of IPH with poly-logarithmic growth in the condition number, see for details [24]. For a similar discretization where the approximation is continuous inside subdomains but discontinuous across subdomains, a substructuring preconditioner was proposed and analyzed for the -version with logarithmic growth in the condition number, see [9].

A favorite preconditioner is the additive Schwarz preconditioner, for which the set of unknowns is partitioned into overlapping or non-overlapping subsets, corresponding to subdomains with maximum diameter . In this paper we only consider the non-overlapping case222There is a subtle difference between overlap at the continuous level of the subdomains, and the discrete level of unknowns, see [14]: no overlap at the level of unknowns means minimal overlap of one mesh size at the continuous level for classical discretizations like finite elements or finite differences. This becomes however even more subtle here with DG discretizations, since the discrete unknowns are coupled through Robin conditions, and no overlap at the level of unknowns really means no overlap at the continuous level, see [15]. and for simplicity study first only two subdomains, a generalization is given in Section LABEL:sec:multi. The non-overlapping two subdomain decomposition results in a natural partitioning of the unknowns . The solution of the linear system by the additive Schwarz method without overlap is equivalent to the block Jacobi iteration

\hb@xt@.01(1.3)

The matrix is also s.p.d. and can be considered as a preconditioner for CG. It can be shown that in this case we have in the absence of a coarse solver; see [12]. Preconditioned CG satisfies then the convergence factor estimate .

On the other hand it has been recently shown in [15] that the block Jacobi iteration in (LABEL:eq:blockJacobi) for an IPH discretization can be viewed as a discretization of a non-overlapping Schwarz method with Robin transmission conditions, i.e.

\hb@xt@.01(1.4)

where , is the interface between the two subdomains and is precisely the penalty parameter of the IPH discretization. This parameter has to be chosen such that it ensures coercivity and optimal approximation properties. For an IPH discretization, we must have for some constant large enough, independent of , and this scaling cannot be weakened, since otherwise coercivity is lost. On the other hand, optimized Schwarz theory suggests that the iteration in (LABEL:eq:DD-DGH) converges faster if , see [13]. In that case for the contraction factor we have while with the choice for IPH, we have .

The challenge is therefore to design a Schwarz algorithm for IPH with convergence factor , while having the same fixed point as the original additive Schwarz or block Jacobi method for IPH. An idea for doing this can be found for Maxwell’s equation in [10]. This approach was also adopted for IPH in [20], where numerical experiments show that the convergence factor is indeed , while maintaining the same fixed point, but there is no convergence analysis.

We provide in this paper a convergence theory for Schwarz methods applied to IPH discretizations and prove these numerical observations. A similar analysis exists for classical FEM using Schur complement formulations and exploiting eigenvalues of the Dirichlet-to-Neumann (DtN) operator, see [22]. Our analysis uses similar DtN arguments, but is substantially different from [22], since in a DG method continuity conditions are imposed only weakly. We focus in our analysis on the -version with polynomial degree one, and do not study the effect of possible jumps in or higher polynomial degree.

Our paper is organized as follows: in Section LABEL:sec:iph we describe two different but equivalent formulations of IPH, and construct a Schur complement system. In Section LABEL:sec:tech we provide mathematical tools to analyze Schwarz methods formulated using Schur complements. In Section LABEL:sec:schwarz we present the additive Schwarz and a new Schwarz algorithm for IPH in a two subdomain setting and prove their convergence with concrete contraction factor estimates. Section LABEL:sec:multi contains a generalization of the algorithms to the multi-subdomain case. We show in Section LABEL:sec:num numerical experiments to illustrate our analysis, and also verify numerically that the new algorithm provides a better preconditioner for Krylov subspace methods: we observe that the contraction factor is which is much faster than the CG solver preconditioned by one level additive Schwarz.

2 Hybridizable Interior Penalty method

This section is devoted to recall the definition of IPH in two different but equivalent forms, namely the primal and hybridizable formulation. We later in Section LABEL:sec:schwarz design and analyze two Schwarz methods for the hybridizable form and show that the first one is slow and equivalent to a block Jacobi method applied to a primal form, i.e. (LABEL:eq:blockJacobi). However the second Schwarz method takes advantage of hybridizable formulation and achieve faster convergence.

IPH was first introduced in [11] as a stabilized discontinuous finite element method and later was studied as a member of the class of hybridizable DG methods in [7]. It has been shown that it is equivalent to a method called Ultra Weak Variational Formulation (UWVF) for the Helmholtz equation; see [19]. IPH also fits into the framework developed in [3] for a unified analysis of DG methods. IPH is further studied in [21] in the context of incompressible flows.

2.1 Notation

We follow the notation introduced in [3]. Let be a shape-regular and quasi-uniform triangulation of the domain . Let be the diameter of an element of the triangulation defined by and . If is an edge of an element, we denote by the length of that edge. The quasi-uniformity of the mesh implies .

We denote by the set of interior edges shared by two elements in , that is

by the set of boundary edges, and all edges by . We introduce the broken Sobolev space where is the Sobolev space in and is a positive integer. Note that is not necessarily continuous across elements. Therefore the element boundary traces of functions in belong to , where can be double-valued on , but is single-valued on .

We now define two trace operators: let and . Then on we define the average and jump operators

where is the unit outward normal from on . It is clear that these operators are independent of the element enumeration. Similarly for a vector-valued function we define on interior edges

On the boundary, we set the average and jump operators to and . We do not need to define and on .

We define a finite dimensional subspace of by

\hb@xt@.01(2.1)

where is the space of polynomials of degree in the simplex . We denote boundary integrals on an edge by

and similarly for volume terms on an element

If is a subset of , we denote the -norm of along by and . Similarly if is a subset of , we denote the -norm of a by .

For we define functions whose restrictions to each element, , are equal to the gradient of . This operator in the literature is called piecewise gradient and is usually denoted by . For the sake of simplicity we use instead of .

2.2 Primal formulation

To simplify our presentation, we set to be a constant and in the model problem (LABEL:eq:pde). Let , then the IPH bilinear form of the model problem (LABEL:eq:pde) is defined as

\hb@xt@.01(2.2)

where , and . Observe that is symmetric. The definition of the IPH bilinear form is different from the classical Interior Penalty (IP) method only in the last term, i.e. the last term in is not present in IP.

There are two natural energy norms which are equivalent at the discrete level. Let then

\hb@xt@.01(2.3)

One can show that they are equivalent at the discrete level by a local application of the inverse inequality (LABEL:eq:invineq).

Proposition 2.1

Let . Then we have

where and independent of and .

The norm provides a natural norm for boundedness and can be used for showing coercivity. The main ingredients for coercivity are the following inequalities which hold for all :

\hb@xt@.01(2.4)

where and are both independent of and but depend on the polynomial degree. This can be obtained from the trace inequality

\hb@xt@.01(2.5)

where is the polynomial degree, for details see [26, 3].

Proposition 2.2

If , for and sufficiently large, then we have

where , and both constants are independent of .

Note that coercivity holds only for and that has to be big enough to result in a positive . Since and come from the trace inequality, we can choose where is the degree of the polynomials in the simplex. Throughout this paper we assume that is chosen big enough to ensure that any term of type (with , independent of and ) is positive.

Having established that is bounded and coercive, we obtain that the following approximation problem has a unique solution: find such that

\hb@xt@.01(2.6)

Assuming the exact solution is regular enough, it can be shown that

i.e. IPH has optimal approximation order [3, 21]. We emphasize that without setting , the coercivity and optimal approximation properties are lost.

2.3 Hybridizable formulation

In this section we exploit the fact that IPH is a hybridizable method. A method is hybridizable if one can eliminate the degrees of freedom inside each element to obtain a linear system in terms of a single-valued function along the edges, say . Not all DG methods have this property, for example classical IP is not hybridizable. A unified hybridization procedure for DG methods has been introduced and studied in [7] where IPH is also included.

We introduce the general setting by decomposing the domain into two non-overlapping subdomains and . Denoting the interface by , we assume , i.e. the cut does not go through any element of the triangulation. This will result in a natural partitioning of into and which do not overlap but share as a boundary; see for an example Figure LABEL:fig:ddmesh.

Figure 2.1: An unstructured mesh with the interface (thick-dashed).

We denote by the maximum diameter of the subdomains and by the diameter of the mono-domain . We assume .

We introduce local spaces on and by

\hb@xt@.01(2.7)

Note that this domain decomposition setting implies . We define on the interface the space of broken single-valued functions by

\hb@xt@.01(2.8)

For the sake of simplicity we denote the restriction of on by . Observe that the trace of on belongs to .

Let and consider the symmetric bilinear form

\hb@xt@.01(2.9)

where

\hb@xt@.01(2.10)

and

\hb@xt@.01(2.11)

This is an IPH discretization of the model problem in and is treated as a Dirichlet boundary. Therefore inherits coercivity and continuity of the original bilinear form, .

The global bilinear form is also coercive at the discrete level, if is sufficiently large, independent of . To see this we introduce an energy norm for all such that

\hb@xt@.01(2.12)

then by definition of for all we have

\hb@xt@.01(2.13)

We can bound the contribution of each subdomain from below separately:

where we used the inverse inequalities (LABEL:eq:warburton) for terms acting on the interface and (LABEL:eq:coerc1) for terms acting inside subdomains. Here is a constant independent of . Note that we proved the coercivity in a subdomain by subdomain fashion by splitting the terms.

Consider the following discrete problem: find such that

\hb@xt@.01(2.14)

which has a unique solution since is coercive on . One can eliminate the interface variable, , and obtain a variational problem in terms of only. It turns out that this coincides with the variational problem (LABEL:eq:varIPH); for a proof see [21].

The advantage of the variational problem (LABEL:eq:varIPH2) is that each subproblem is communicating through the auxiliary unknown . Therefore we can eliminate the interior unknowns, , and obtain a Schur complement system. If we test (LABEL:eq:varIPH2) with , , and assume that is known, we obtain a local problem: find such that

\hb@xt@.01(2.15)

This is an IPH discretization of the continuous problem

However the boundary condition on is imposed weakly and therefore in the strong sense, see [7, 15, 21].

2.4 Schur complement formulation

We choose nodal basis functions for and denote the space of degrees of freedom (DOFs) of by and similarly for subspaces by . The variational form in (LABEL:eq:varIPH) is equivalent to the linear system . is the system matrix and are the corresponding DOFs of the approximation . We can partition into where corresponds to DOFs of . Then we can arrange the entries of and rewrite the linear system as

\hb@xt@.01(2.16)

We use nodal basis functions for and denote by the corresponding DOFs for . Then the variational form (LABEL:eq:varIPH2) can be written as

\hb@xt@.01(2.17)

where . Since this matrix is s.p.d. and the same holds also for its diagonal blocks, we can form a Schur complement system. We define and Then the Schur complement system reads

\hb@xt@.01(2.18)
Definition 2.3 (discrete harmonic extension)

For all , we denote by the discrete harmonic extension into ,

\hb@xt@.01(2.19)

The corresponding is called generator. In other words is an approximation obtained from the IPH discretization in using as Dirichlet data; i.e. .

The following result shows that an application of can be viewed as finding the harmonic extension, , and then evaluating a “Robin-like trace” on the interface.

Proposition 2.4

Let and define its harmonic extension by . Then for all .

Proof. Let . Then by definition of and we have

for all , which completes the proof, since .       

3 Properties of the Schur complement and technical tools

The main goal of this section is to provide estimates for the minimum and maximum eigenvalues of the and for . We use the estimate for the operators to prove convergence of the Schwarz method and provide the contraction factor later in Section LABEL:sec:schwarz. In particular we prove in this section that the following estimates hold for all :

\hb@xt@.01(3.1)
\hb@xt@.01(3.2)

where all constants are positive and independent of , and . Since and are symmetric, we can use Rayleigh quotient arguments and obtain an estimate for the minimum and maximum eigenvalues. One can also obtain an estimate with polynomial degree dependency using the techniques of this section.

The only constraint on the shape of the subdomains is a star-shape assumption. To prove the above estimates we need trace and Poincaré inequalities for totally discontinuous functions. The following trace estimate is due to Feng and Karakashian [12, Lemma 3.1]. The Poincaré inequality is due to Brenner, see [5].

Lemma 3.1 (Trace inequality)

Let be a star-shape domain with diameter , and triangulation . Then, for any , we have

Lemma 3.2 (Poincaré inequality)

Let be an open connected polygonal domain with diameter , and triangulation . Then, for any we have

where is a measurable subset of with nonzero measure.

3.1 Eigenvalue estimates for

In order to obtain estimates for the eigenvalues of the operator, we first recall Definition LABEL:HarmExtDef of a harmonic extension: is called harmonic extension of if it satisfies . Now multiplying this relation by from left we get

where we used , and the definition of . Hence if then we have

\hb@xt@.01(3.3)

Now recall that is coercive and bounded over , therefore . Thus if we relate the energy norm of the harmonic extension, , to the -norm of we obtain the desired estimate (LABEL:eq:Best). More precisely we can show that the estimate

\hb@xt@.01(3.4)

holds, where and are constants independent of . Observe that while the upper bound estimate in (LABEL:eq:Best) is less than one. We show later how one can obtain a sharp upper bound estimate as in (LABEL:eq:Best).

Let us start with the lower bound of inequality (LABEL:eq:harmIneq). First we introduce an extension by zero operator which is defined for all as

For a graphical illustration see Figure LABEL:fig:extzero.

Figure 3.1: Illustration of the extension by zero, , for elements which share an edge with the interface, e.g. , and those which do not, e.g. .

Note that there are elements like which physically share a node and not an edge with the interface, but we leave in to be zero. More precisely, only those elements which share an edge with the interface are non-zero.

We show in the Appendix, see also [23], that in an element, , with an edge we have

\hb@xt@.01(3.5)

where , and and all are independent of . This yields the following result which relates the energy of the extension by zero to its -norm on the interface.

Lemma 3.3

Let and be its extension by zero into . We have

where .

Proof. First note that by definition and are non-zero only on those elements which share an edge with the interface. We call them . Then we have

which completes the proof with .       

Now we are able to relate the energy of a harmonic extension, , to the -norm of on the interface.

Lemma 3.4

Let and be its harmonic extension into . Then we have

where .

Proof. Since is the harmonic extension of , it satisfies (LABEL:eq:harmonicsat) (with ). Let . Then by definition of we have

Note that . We can bound the right-hand side from below, therefore

which is positive if and sufficiently large. By continuity of we have

Note that we are able to use instead of since we work with discrete spaces. An application of Lemma LABEL:lemma:Zoptphi completes the proof with .       

The upper bound in (LABEL:eq:harmIneq) can be obtained much easier using coercivity of the .

Lemma 3.5

Let and be its harmonic extension into . Then we have

where .

Proof. Since is the harmonic extension of , it satisfies (LABEL:eq:harmonicsat) (with ). Using the fact that is coercive we have