Accessible Published by De Gruyter March 22, 2018

Hybrid Discontinuous Galerkin Discretisation and Domain Decomposition Preconditioners for the Stokes Problem

Gabriel R. Barrenechea, Michał Bosy ORCID logo, Victorita Dolean ORCID logo, Frédéric Nataf and Pierre-Henri Tournier


Solving the Stokes equation by an optimal domain decomposition method derived algebraically involves the use of nonstandard interface conditions whose discretisation is not trivial. For this reason the use of approximation methods such as hybrid discontinuous Galerkin appears as an appropriate strategy: on the one hand they provide the best compromise in terms of the number of degrees of freedom in between standard continuous and discontinuous Galerkin methods, and on the other hand the degrees of freedom used in the nonstandard interface conditions are naturally defined at the boundary between elements. In this paper, we introduce the coupling between a well chosen discretisation method (hybrid discontinuous Galerkin) and a novel and efficient domain decomposition method to solve the Stokes system. We present the detailed analysis of the hybrid discontinuous Galerkin method for the Stokes problem with nonstandard boundary conditions. The full stability and convergence analysis of the discretisation method is presented, and the results are corroborated by numerical experiments. In addition, the advantage of the new preconditioners over more classical choices is also supported by numerical experiments.

1 Introduction

Discontinuous Galerkin (DG) methods have been first introduced in the early 1970s [25] and they have benefited of a wide interest from the scientific community. The main advantages of these methods are their generality and flexibility as they can be used for a variety of partial differential equations on unstructured meshes. Moreover, they can preserve local properties such as mass and momentum conservation while ensuring a high-order accuracy. However, the cost of these advantages is a larger amount of degrees of freedom in comparison to the continuous Galerkin methods [17] for the same approximation order.

A good compromise between the previous methods, while preserving the high order, are the hybridised versions of DG using divergence conforming spaces such as Raviart–Thomas (RT) and Brezzi–Douglas–Marini (BDM) [2]. These methods are a subset of the hybrid discontinuous Galerkin (HDG) methods, introduced in [8] for second-order elliptic problems, and extended to three-dimensional Stokes equation in [7]. The authors present there the mixed formulation of HDG methods defined locally on each element. They consider many types of boundary conditions that involve normal and tangential velocity, pressure, and tangential stress. The formulations of the methods are similar, the only difference lies in the choice of the numerical traces. A refined analysis of HDG methods for the Stokes equation with Dirichlet boundary conditions was presented in [9], where optimal convergence in the HDG norm is proven, and superconvergent postprocessings are introduced. An alternative strategy to approximate the Stokes problem was used in [16], where a combination of HDG and a symmetric interior penalty (somehow linked to the work [27]) is presented and analysed.

In a further paper [15] this approach is extended to Darcy, and coupled Darcy–Stokes flows. The new formulation includes different degrees of polynomials for finite element spaces associated with different variables.

In all HDG discretisations, one important component is the Lagrange multiplier defined on the inter element edges (faces in 3D). In the work [22] the authors approximate the Navier–Stokes equation using H(div)-conforming spaces (rather than the usual H1-conforming ones) together with an HDG approach. Now, in order to reduce the number of degrees of freedom needed for the Lagrange multiplier, the authors introduce an edge-wise projection into a finite element space of lower degree. This projection was also shown to be of paramount importance to establish the connection between the hybrid high-order [11] and the HDG methods presented in [6].

Even considering the decrease of the number of degrees of freedom provided by the HDG methods with projection, modern applications lead to linear systems whose size is too large to allow the use of direct solvers. Thus, parallel solvers are becoming increasingly important in scientific computing. A natural paradigm to take advantage of modern parallel architectures is the Domain Decomposition method, see e.g. [28, 24, 31, 12]. Domain decomposition methods are iterative solvers based on a decomposition of a global domain into subdomains. At each iteration, one (or two) boundary value problem(s) are solved in each subdomain and the continuity of the solution at the interfaces between subdomains is only satisfied at convergence of the iterative procedure. The partial differential equation is the one of the global problem.

For Additive Schwarz methods and Schur complement methods, the boundary conditions on the interfaces between subdomains, or, interface conditions (IC), are Dirichlet or Neumann boundary conditions. For scalar equations there is a consensus on this choice of IC. On the other hand, for systems of differential equations, such as the Stokes or incompressible elasticity equations, alternatives like normal velocity-tangential flux (NVTF) or tangential velocity-normal flux (TVNF) IC should be superior to the pure velocity (Dirichlet like) or pure stress (Neumann like) IC, see [12, Section 6.6] and the references therein. In [19], this choice was motivated by symmetry considerations, while in [14, 5, 4], they were obtained by an analysis of the systems of partial differential equations by symbolic techniques, mostly linked to the Smith factorization [29]. Similar attempts to derive more intrinsic IC to the nature of the equation to solve were derived [13] for the Euler system.

Due to the difficulty of implementing these IC, previous numerical tests were restricted to decompositions where the boundaries of subdomains are rectilinear so that the normal to the interface is easy to define. The underlying domain decomposition method was a Schur complement method. That is mainly the reason we have considered and analysed a specific HDG method where this kind of degrees of freedom are naturally present.

In this paper, we want to combine appropriate HDG discretisation and the associated domain decomposition methods mentioned above using nonstandard IC. The combination of the two is meant to provide a competitive solving strategy for this kind of partial differential equation system. A different, but somewhat related, approach can be found in [1] where a DG type discretisation is coupled to a discrete Helmholtz decomposition to propose some preconditioners.

Concerning the HDG discretisation, in approach similar to the one presented in this work, but using completely discontinuous spaces, is given in [23]. Our analysis is related to the one in that paper, but the method presented herein uses H(div)-conforming spaces, which implies in turn that the Lagrange multipliers are scalar valued. The combination of these two facts reduces the number of degrees of freedom significantly. In addition, the use of nonstandard boundary conditions (motivated by the newly defined domain decomposition preconditioners) makes the analysis somehow more involved.

The rest of the paper is organised as follows. To start with, we introduce the problem and notation in Section 2. In Section 3, we present the hybridisation of a symmetric interior penalty Galerkin method that allows us to impose the TVNF and NVTF boundary conditions in quite a natural way. The formulation is similar to the one from [22] with Dirichlet boundary conditions. In addition to different kinds of boundary conditions, we include a local edge-based projection in order to reduce the dimension of the space of Lagrange multipliers. Our analysis follows the one from [22] (see also [21] for a more detailed version). Thanks to the HDG discretisation, we can consider domain decomposition methods with arbitrary shape of the interfaces and Schwarz type methods. In Section 4, the Additive Schwarz methods are defined at the algebraic level. Section 5 contains the numerical results, including the convergence validation of the HDG method and a comparison of the domain decomposition preconditioners. Finally, some conclusions are drawn.

2 Notation and Preliminary Results

Let Ω be an open polygonal domain in 2 with Lipschitz boundary Γ:=Ω. We use boldface font for tensor or vector variables, e.g. 𝒖 is a velocity vector field. The scalar variables will be italic, e.g. p denotes a pressure scalar value. We define the stress tensor 𝝈:=ν𝒖-p𝑰 and the flux 𝝈𝒏:=𝝈𝒏. In addition, we denote normal and tangential components as follows: un:=𝒖𝒏, ut:=𝒖𝒕, σnn:=𝝈𝒏𝒏, σnt:=𝝈𝒏𝒕, where 𝒏 is the outward unit normal vector to the boundary Γ and 𝒕 is a vector tangential to Γ such that 𝒏𝒕=0.

For DΩ, we use the standard L2(D) space with the following norm:

fD2:=Df2𝑑𝒙for all fL2(D).

Let us define, for m, the following Sobolev spaces:

Hm(D):={vL2(D):𝜶vL2(D) for all |𝜶|m},

where, for 𝜶=(α1,α2)2, |𝜶|=α1+α2, and


In addition, we will use following standard semi-norm and norm for the Sobolev space Hm(D):

|f|Hm(D)2:=|𝜶|=m𝜶fD2,fHm(D)2:=k=0m|f|Hk(D)2for all fHm(D).

In this work, we consider the two-dimensional Stokes problem

(2.1){-νΔ𝒖+p=𝒇in Ω,𝒖=0in Ω,

where 𝒖:Ω¯2 is the unknown velocity field, p:Ω¯ the pressure, ν>0 the viscosity, which is considered to be constant, and 𝒇[L2(Ω)]2 is a given function. For gL2(Γ) we consider two types of boundary conditions:

  • tangential-velocity and normal-flux (TVNF)

    (2.2){σnn=gon Γ,ut=0on Γ,
  • normal-velocity and tangential-flux (NVTF)

    (2.3){σnt=gon Γ,un=0on Γ,

which together with (2.1) define two boundary value problems. We will detail the analysis for the TVNF boundary value problem

(2.4){-νΔ𝒖+p=𝒇in Ω,𝒖=0in Ω,σnn=gon Γ,ut=0on Γ,

since considering the NVTF boundary conditions (2.3) instead is very similar. We will just add a remark when necessary to stress the differences between them. The restriction to homogeneous Dirichlet conditions on ut is made only to simplify the presentation.

Let {𝒯h}h>0 be a regular family of triangulations of Ω¯ made of triangles. For each triangulation 𝒯h, h denotes the set of its edges. In addition, we set hK:=diam(K) for each K𝒯h, and h:=maxK𝒯hhK. We define the following Sobolev spaces on the triangulation 𝒯h and the set of all edges in h:

L2(h):={v:v|EL2(E) for all Eh},
Hm(𝒯h):={vL2(Ω):v|KHm(K) for all K𝒯h}for m,

with the corresponding broken norms.

The following results will be very useful in what follows.

Lemma 2.1 (Inverse and Trace Inequalities).

There exist C,Cmax>0, independent of hK, such that for every KTh and every polynomial function v in K the following inequalities hold:


Moreover, there exists C>0, independent of hK, such that for any vH1(K), the following local trace inequality holds:



For (2.5) see [17, Lemma 1.138] and for (2.6) see [10, Lemma 1.46]. The discrete trace inequality (2.7) follows by standard scaling arguments. ∎

We now introduce the finite element spaces used for discretisation. First, we present them for the particular case in which we consider the TVNF boundary condition. For this case for k1, we discretise the velocity using the Brezzi–Douglas–Marini space (see [2, Section 2.3.1]) given by

BDMhk:={𝒗𝒉H(div,Ω):𝒗𝒉|K[k(K)]2 for all K𝒯h}.

The choice of the above finite element spaces is motivated by two requirements. On the one hand, we have imposed that the finite element space must be conforming in H(div,Ω), while on the other hand, it has to satisfy the inf-sup condition with respect to a broken H1-norm of the velocity. The BDM space appears then as a natural alternative satisfying both these requirements. In fact, other popular alternatives, such as the Raviart–Thomas space, while conforming in H(div,Ω) do not satisfy the inf-sup condition with respect to the norm used in this work.

In addition, for 1mk+1 we denote by Πk:[Hm(Ω)]2BDMhk the BDM projection defined in [2, Section 2.5]. The HDG formulation includes a Lagrange multiplier over the internal edges. In order to propose a discretisation with fewer degrees of freedom, we discretise the Lagrange multiplier u~ using the spaces

Mhk-1:={v~hL2(h):v~h|Ek-1(E) for all Eh},
Mh,0k-1:={v~hMhk-1:v~h=0 on Γ}.

Furthermore, we introduce for all Eh the L2(E)-projection ΦEk-1:L2(E)k-1(E) defined as follows. For every w~L2(E), ΦEk-1(w~) is the unique element of k-1(E) satisfying

EΦEk-1(w~)v~h𝑑s=Ew~v~h𝑑sfor all v~hk-1(E),

and we define Φk-1:L2(h)Mhk-1 by Φk-1|E:=ΦEk-1 for all Eh.

Let us denote 𝑽𝒉:=BDMhk×Mh,0k-1. The pressure is discretised using the following space:

Qhk-1:={qhL2(Ω):qh|Kk-1(K) for all K𝒯h}.

In addition, we define the local L2(K)-projection ΨKk-1:L2(K)k-1(K) for each K𝒯h defined as follows. For every wL2(K), ΨKk(w) is the unique element of k-1(K) satisfying

KΨKk-1(w)vh𝑑x=Kwvh𝑑xfor all vhk-1(K).

We will also use the following results.

Lemma 2.2 (Approximation Results).

There exists C>0, independent of hK, such that for all 𝐯[Hm(K)]2 and vHm(K), 1mk+1, the following interpolation estimates hold:

  • local Brezzi–Douglas–Marini approximation

  • trace L2-projection approximation

  • local L2-projection approximation



For (2.8) see [2, Proposition 2.5.1], for (2.9) see [18, Lemma III.2.10], and for (2.10) see the proof of [17, Theorem 1.103]. ∎

3 Hybrid Discontinuous Galerkin Method

In this section, we introduce the HDG method proposed in this work, study its well-posedness, and analyse its error.

3.1 The Discrete Problem

From now on we will use to denote the element-wise gradient. First, we multiply the first equation from (2.1) by a test function 𝒗𝒉BDMhk and integrate by parts. This gives


Since the normal and tangential vectors are perpendicular (𝒏𝒕=0), we can split (3.1) as


For the solution of the Stokes problem (2.1), 𝝈𝒏 is continuous across all interior edges. Moreover, since 𝒗𝒉BDMhk, we see that (𝒗𝒉)n is continuous across all interior edges, and thus we can rewrite (3.2) as


In addition, since 𝝈𝒏 is continuous across all interior edges, we have

K𝒯hKσntv~h𝑑s=0for all v~hMh,0k-1,

and we can add this to (3.3) to get


Denoting u~=ut on h leads to (ut-u~)=Φk-1(ut-u~)=0 on h. Applying the boundary conditions (2.2), we can rewrite (3.4) as


where τ>0 is a stabilisation parameter. Hence, we define the velocity bilinear form a:𝑽𝒉×𝑽𝒉 as


where ε{-1,1} and τ>0 is a stabilisation parameter, and b:𝑽𝒉×Qhk-1 as


With these definitions we propose the HDG method for the TVNF boundary value problem (2.4):

Find (𝒖𝒉,u~h,ph)𝑽𝒉×Qhk-1 such that for all (𝒗𝒉,v~h,qh)𝑽𝒉×Qhk-1,


Remark 3.1.

The use of H(div)-conforming spaces not only decreases the number of degrees of freedom in comparison to [23], but leads as well to a simpler bilinear form b, as the jump terms appearing in that reference are no longer needed for stability or consistency.

3.2 Well-Posedness of the Discrete Problem

Let us consider the following semi-norm:


Lemma 3.2.

The semi-norm |||||| defined by (3.7) is a norm on 𝐕𝐡.


Since |||||| is a semi-norm, we only need to show that

|||(𝒘𝒉,w~h)|||=0𝒘𝒉=𝟎 and w~h=0.

Let us suppose (𝒘𝒉,w~h)𝑽𝒉 and |||(𝒘𝒉,w~h)|||=0. Then 𝒘𝒉=0 in all K𝒯h, and thus 𝒘𝒉|K=𝑪K for all K𝒯h. Now, since 𝒘𝒉[0(K)]2 in every K, we have

Φk-1((𝒘𝒉)t-w~h)K=0(𝒘𝒉)t=w~h in each Eh.

Since w~h is single valued on all the edges in h, we have that (𝒘𝒉)t is continuous in Ω. Moreover, 𝒘𝒉 belongs to BDMhk, so (𝒘𝒉)n is also continuous in Ω. Then, 𝒘𝒉 is continuous in Ω, and thus 𝒘𝒉=𝑪2 in Ω. Finally,

(𝒘𝒉)t=(𝑪)t=0 on Γ𝒘𝒉=𝟎 in Ω,

which finishes the proof since w~h=(𝒘𝒉)t on every edge. ∎

Lemma 3.3.

There exists C>0 such that, for all (𝐰,w~),(𝐯,v~)[H1(Ω)H2(Th)]2×L2(Eh) and qL2(Ω), we have



Let us start with (3.8). Using the Cauchy–Schwarz inequality, we get


Therefore, using the triangle inequality and the trace L2-projection approximation (2.9), we get


Thus, using the Cauchy–Schwarz inequality gives


Finally, we get (3.8) for C=(2+c1+c2). The continuity (3.9) is analogous. ∎

To show the well-posedness of (3.6) we need the ellipticity of the bilinear form a and an inf-sup condition for the bilinear form b. We start by showing that a is elliptic with respect to ||||||.

Lemma 3.4.

There exists α>0 such that for all (𝐯𝐡,v~h)𝐕𝐡,


If ε=-1 in the definition (3.5), then inequality (3.10) only holds under the additional hypothesis of τ being large enough. If ε=1 in (3.5), then inequality (3.10) holds for arbitrary τ.


First, since 𝒏𝒗𝒉|E[k-1(E)]2 for all Eh, we have


To bound the middle term in terms of the other two, we consider two cases.

Case 1: ε=1. Then (3.11) reduces to


It only remains to show that the right-hand side of (3.12) is an upper bound (up to a constant) for the norm |||||| given by (3.7). Using the discrete trace inequality (2.6), we get


and then


which proves (3.10) with α=1/(1+Cmax2).

Case 2: ε=-1. Then (3.11) becomes


Using the Cauchy–Schwarz inequality, we deduce


Since 𝒗𝒉BDMhk is a piecewise polynomial, we can apply the discrete trace inequality (2.6) to the second term, followed by Young’s inequality to arrive at


Finally, if we suppose τ>2Cmax2 and set


then using (3.13), we get (3.10) for α=C/(1+Cmax2). ∎

The next step towards stability is proving the inf-sup condition for b(,).

Lemma 3.5.

There exists β>0 independent of hK such that

sup(𝒗𝒉,v~h)𝑽𝒉b((𝒗𝒉,v~h),qh)|||(𝒗𝒉,v~h)|||βνqhΩfor all qhQhk-1.


According to the Fortin criterion, see [17, Lemma 4.19], we need to prove that there exists a Fortin operator 𝚷:[H1(Ω)]2𝑽𝒉 such that for every 𝒗[H1(Ω)]2 the following conditions hold:

(3.14)b((𝒗,v~),qh)=b(𝚷(𝒗),qh)for all qhQhk-1,

Let 𝒗[H1(Ω)]2 and let us consider the operator 𝚷(𝒗):=(Πk(𝒗),Φk-1(vt)). It is well known, see [2, Section 2.5], that Πk satisfies (3.14). To prove (3.15) we denote (𝒘𝒉,w~h):=𝚷(𝒗). Then using the discrete trace inequality (2.6) and the fact that the projection is a bounded operator, we get


Applying the triangle inequality for the last term of (3.16), we arrive at


Using the stability of Πk leads to


Using (2.8) and the local trace inequality (2.7), we get


Finally, using the trace L2-projection approximation (2.9) for the third term, we get


Then collecting (3.17), (3.18) and (3.19), we obtain (3.15) with


which finishes the proof. ∎

Using the last two results and the standard theory of variational problems with constraints [2, Section 4.2], we deduce there exists a unique solution of (3.6). In addition, thanks to the derivation carried out in Section 3.1, method (3.6) is consistent as the following result shows.

Lemma 3.6 (Consistency).

Let (𝐮,p)[H1(Ω)H2(Th)]2×L2(Ω) be the solution of problem (2.4) and u~=ut on all edges of Eh. If (𝐮𝐡,u~h,ph)𝐕𝐡×Qhk-1 solves (3.6), then for all (𝐯𝐡,v~h,qh)𝐕𝐡×Qhk-1 the following holds:


3.3 Error Analysis

In this section, we present the error estimates for the method. These estimates are proved using the norm


The first step is the following version of Cea’s lemma.

Lemma 3.7.

Let (𝐮,p)[H1(Ω)H2(Th)]2×L2(Ω) be the solution of (2.4), u~=ut on all edges in Eh, and let (𝐮𝐡,u~h,ph)𝐕𝐡×Qhk-1 be the solution of (3.6). Then there exists C>0, independent of h and ν, such that



Let us denote


Using Lemmas 3.4 and 3.5, and [17, Proposition 2.36], we get the following stability for B: There exists βB>0, independent of h and ν, such that for all (𝒗𝒉,v~h,qh)𝑽𝒉×Qhk-1 there exists (𝒘𝒉,w~h,rh)𝑽𝒉×Qhk-1 such that |||(𝒘𝒉,w~h,rh)|||h=1, and


Now Lemma 3.3 yields the following continuity for B: there exists CB>0 such that


Let (𝒗𝒉,v~h,qh)𝑽𝒉. Then, using Lemma 3.6, the triangle inequality, (3.21) and (3.22), we arrive at


Thus, we get (3.20) with C:=1+CB/βB. ∎

We prove the following error estimate using standard interpolation estimates.

Lemma 3.8 (HDG Error).

Let us assume (𝐮,p)[H1(Ω)Hk+1(Th)]2×Hk(Th) is the solution of (2.4), and u~=ut on all edges in Eh. If (𝐮𝐡,u~h,ph)𝐕𝐡×Qhk-1 solves the discrete problem (3.6), then there exists C>0, independent of h, such that



Let us consider the Fortin operator 𝚷 defined in the proof of Lemma 3.5. If 𝚷(𝒖)=(𝒘𝒉,w~h), then by using the triangle inequality and boundedness of the projection Φk-1, we get


For the first term from (3.24), we use the BDM approximation (2.8) to get


Next we use the local trace inequality (2.7) to get


Let k𝒖 be the usual Lagrange interpolant of degree k of 𝒖 (see [17, Example 1.31]). Using the triangle inequality followed by the local inverse inequality (2.5), the local Lagrange approximation [17, Example 1.106] and (2.8), we see that (3.26) becomes


For the third term in (3.24), we use (2.7) and (2.8) to get


The last term in (3.24) is bounded using (2.9) as follows:


Finally, the local L2-projection approximation (2.10) gives


Thus, putting together (3.24) with (3.25), (3.27)–(3.30), and using the shape regularity of the mesh, we get




and the result (3.23) follows from Lemma 3.7. ∎

3.4 NVTF Boundary Conditions

As we mentioned before, the analysis in case of NVTF boundary conditions (2.3) is similar. Thus, we just highlight the main differences. If we consider NVTF boundary conditions (2.3), then to discretise the velocity we use the BDM space

BDMh,0k:={𝒗𝒉BDMhk:(𝒗𝒉)n=0 on Γ}.

For the Lagrange multiplier we use the polynomial space Mhk-1, while the pressure is discretised using


The product space then becomes 𝑽𝒉:=BDMh,0k×Mhk-1 and we pose the following discrete problem:

Find (𝒖𝒉,u~h,ph)𝑽𝒉×Qh,0k-1 such that for all (𝒗𝒉,v~h,qh)𝑽𝒉×Qh,0k-1,


In obtaining (3.31) the only difference step in the derivation is that now (3.4) is replaced by


Concerning the analysis, the proofs of all the results presented in the last sections remain essentially unchanged.

4 The Domain Decomposition Preconditioner

Let us assume that we have to solve the linear system


where 𝐀 is the matrix arising from discretisation of the Stokes equations on the domain Ω, 𝑼 is the vector of unknowns and 𝑭 is the right-hand side. To accelerate the performance of an iterative Krylov method applied to this system we will consider domain decomposition preconditioners which are naturally parallel [12, Chapter 3]. They are based on an overlapping partition of the computational domain.

Let {𝒯h,i}i=1N be a partition of the triangulation 𝒯h. For an integer value l0, we define an overlapping decomposition {𝒯h,il}i=1N such that 𝒯h,il is a set of all triangles from 𝒯h,il-1 and all triangles from 𝒯h𝒯h,il-1 that have non-empty intersection with 𝒯h,il-1, and 𝒯h,i0=𝒯h,i. With this definition the width of the overlap will be 2l. Furthermore, if Wh stands for the finite element space associated to 𝒯h, then Wh,il is the local finite element space on 𝒯h,il that is a triangulation of Ωi.

Let 𝒩 be the set of indices of degrees of freedom of Wh and let 𝒩il be the set of indices of degrees of freedom of Wh,il for l0. Moreover, we define the restriction operator 𝐑i:WhWh,il as a rectangular matrix |𝒩il|×|𝒩| such that if 𝑽 is the vector of degrees of freedom of vhWh, then 𝐑i𝑽 is the vector of degrees of freedom of Wh,il in Ωi. Abusing notation, we denote by 𝐑i both the operator and its associated matrix. The extension operator from Wh,il to Wh and its associated matrix are both given by 𝐑iT. In addition, we introduce a partition of unity 𝐃i as a diagonal matrix |𝒩il|×|𝒩il| such that


where 𝐈𝐝|𝒩|×|𝒩| is the identity matrix.

We are ready to present the first preconditioner, called Restricted Additive Schwarz (RAS) [3], given by


We also introduce a new preconditioner that is a modification of the above one. The modification is similar to the Optimized RAS [30], however we do not use Robin IC. For this, let 𝐁i be the matrix associated to a discretisation of (2.1) in Ωi where we impose either TVNF (2.2) or NVTF (2.3) boundary conditions in Ωi. Then, the preconditioner reads


Remark 4.1.

The improvement of convergence in the case of Optimized RAS depends on the choice of the parameter. This parameter depends on the problem and discretisation. The big advantage of MRAS preconditioners is that they are parameter-free.

4.1 Partition of Unity

The above definitions of the preconditioners can be associated with any discretisation of the problem. However, each discretisation involves the construction of a relevant partition of unity 𝐃i, i=1,,N. We discuss here the construction of 𝐃i when the problem (2.1) is discretised by the HDG method in case k=1, either with TVNF boundary conditions (3.6), or NVTF boundary conditions (3.31). Let us introduce the piecewise linear functions χ~il of 𝒯h such that

χ~il={1on all nodes of 𝒯h,i0,0on other nodes.

Now we define the piecewise linear functions χil of 𝒯h,il as follows:


Obviously i=1Nχil=1. We define the partition of unity matrix 𝐃i as a block diagonal matrix where the first block 𝐃iBDM is associated with BDMh1, the second block 𝐃iM with Mh0, and the third block 𝐃iQ with Qh0. The degrees of freedom of the BDM elements are associated with the normal components on the edges of the mesh. For these finite elements, the diagonal of 𝐃iBDM is a vector obtained by interpolating χil at the two points of the edges. The degrees of freedom of the Lagrange multiplier finite elements are associated with the edges of the mesh. For these finite elements, the diagonal of 𝐃iM is a vector obtained by interpolating χil at the midpoints of the edges. For pressure finite elements, the diagonal of 𝐃iQ is a vector obtained by interpolating χil at the midpoints of the elements.

5 Numerical Results

In this section, we present a series of numerical experiments aimed at confirming the theory developed in Section 3, and to give a computational comparison of the preconditioners discussed in the previous section. All experiments have been made by using FreeFem++ [20], which is a free software specialised in variational discretisations of partial differential equations.

5.1 Convergence Validation

The computational domain for both test cases considered here is the unit square Ω=(0,1)2. We present the results for k=1, that is, the discrete space is given by BDMh1×Mh,00×Qh0 for TVNF boundary conditions and BDMh,01×Mh0×Qh,00 for NVTF boundary conditions. We test both the symmetric method (ε=-1) and the non-symmetric method (ε=1). For both cases we have followed the recommendation given in [21, Section 2.5.2] and taken τ=6.

The first example aims at verifying the formulation with TVNF boundary conditions (3.6). We choose the right-hand side 𝒇 and the boundary datum g such that the exact solution is given by


In Figures a and b we show the results of the usual convergence order tests for the symmetric case and the non-symmetric case by plotting in log-log scale the error as a function of the size of the mesh. We notice that they validate the theory from Section 3.3. In addition, an optimal h2 convergence rate is observed for 𝒖-𝒖𝒉Ω. The proof of this last fact is lacking. This is due to the lack of regularity results concerning the solution of the continuous problem with these nonstandard boundary conditions, and hence the usual Aubin–Nitsche approach can not be advocated.

Figure 1

Error convergence of the HDG method with TVNF boundary condition – the first example.

(a) Symmetric bilinear form (ε=-1{\varepsilon=-1})


Symmetric bilinear form (ε=-1)

(b) Non-symmetric bilinear form (ε=1{\varepsilon=1})


Non-symmetric bilinear form (ε=1)

In the second example we want to verify the formulation with NVTF boundary conditions (3.31). In this case, the right-hand side 𝒇 and the boundary datum g are chosen for the following exact solution:


The convergence error with respect to the size of the mesh is depicted on the log-log plots for the symmetric case and the non-symmetric case in Figures a and b, respectively. We can see that they not only validate the theory from Section 3.3, but also show an optimal h2 convergence rate for 𝒖-𝒖𝒉Ω.

Figure 2

Error convergence of the HDG method with NVTF boundary condition – the second example.

(a) Symmetric bilinear form (ε=-1{\varepsilon=-1})


Symmetric bilinear form (ε=-1)

(b) Non-symmetric bilinear form (ε=1{\varepsilon=1})


Non-symmetric bilinear form (ε=1)

5.2 Comparison of Different Domain Decomposition Preconditioners

In this section, we compare the standard RAS preconditioner (4.1) with the newly introduced preconditioners, that is, the ones based on nonstandard IC. We call them MRAS preconditioners (4.2) and more precisely TVNF-MRAS for which 𝐁i is the matrix arising from the discretisation of (2.1) in Ωi with IC (2.2) on Ωi, and NVTF-MRAS for which 𝐁i is the matrix arising from the discretisation of (2.1) in Ωi with IC (2.3) on Ωi. As we mentioned before, our preconditioners do not depend on the used discretisation, that is why we add also similar preconditioners but based on a more standard discretisation, that is, the lowest order Taylor–Hood discretisation [18, Section II.4.2]. In all cases, they are used in conjunction with a Krylov iterative solver such as GMRES [26]. In addition, N stands for the number of subdomains. In all tables we present the number of iterations needed to achieve an euclidean norm of the error (with respect to the one domain solution) smaller than 10-6. We have implemented the RAS preconditioner (4.1) and the MRAS (4.2), using both TVNF and NVTF interface conditions.

We start with the second example from the previous section. However, now we consider the symmetric (ε=-1) formulation with TVNF boundary conditions (3.6). The mesh is uniform and contains 125 000 triangles for a total of 565 003 degrees of freedom for the Taylor–Hood discretisation and 689 000 degrees of freedom for the HDG discretisation. We use a random initial guess for the GMRES iterative solver. The overlapping decomposition into subdomains can be uniform (Unif) or generated by METIS (MTS) and it has two layers of size h in the overlap.

Table 1

Preconditioners comparison – the first test case.


The first thing that we can notice from Table 1 is the important convergence improvement in case of RAS applied to a system resulting from a HDG discretisation in comparison to the RAS applied to the system resulting from the Taylor–Hood discretisation despite the fact that the number of degrees of freedom is slightly bigger in the first case. The change in discretisation presumably leads to better conditioned systems to solve. Also the MRAS preconditioner with both discretisations performs better than the standard RAS method which fully justifies the use of the new IC no matter the discretisation method. Moreover, as expected, the number of iterations increases with respect to the number of the subdomains and this behaviour is common to the three preconditioners. It is worth noticing that this increase is slower than the expected linear one.

We also plot the convergence of the error for the different discretisations in Figure 3. We observe that in all cases the MRAS preconditioner (4.2) shortens the plateau region in the convergence curves significantly which leads, automatically, to an important reduction in the number of iterations.

Figure 3

Convergence of error in the 64 subdomains case – the first test case.

(a) Taylor–Hood, uniform decomposition


Taylor–Hood, uniform decomposition

(b) HDG, uniform decomposition


HDG, uniform decomposition

(c) Taylor–Hood, METIS decomposition


Taylor–Hood, METIS decomposition

(d) HDG, METIS decomposition


HDG, METIS decomposition

Now we consider the Poiseuille problem and we choose the right-hand side 𝒇 and the TVNF boundary condition such that the exact solution is given by


We use the same mesh as in the previous case and a random initial guess for the GMRES iterative solver, but this time the overlapping decomposition has three layers of size h in the overlap.

Table 2

Preconditioners comparison – the Poiseuille problem.


The conclusions stay the same as in previous example since the results from Table 2 are similar to the previous ones. We consider a different problem, however on the same mesh. Hence the global matrix is the same in both cases. Thus, we can notice a reduction in the number of iterations caused by the increase of the width of the overlap.

We observe in Figure 4 that once again the MRAS preconditioner (4.2) leads to a reduction of the number of iterations in all cases.

Figure 4

Convergence of error in the 144 subdomains case – the Poiseuille problem.

(a) Taylor–Hood, uniform decomposition


Taylor–Hood, uniform decomposition

(b) HDG, uniform decomposition


HDG, uniform decomposition

(c) Taylor–Hood, METIS decomposition


Taylor–Hood, METIS decomposition

(d) HDG, METIS decomposition


HDG, METIS decomposition

Figure 5

Numerical solution of the T-shaped domain problem.

(a) Velocity field 𝒖{\boldsymbol{u}}


Velocity field 𝒖

(b) Pressure p


Pressure p

The last example is on a T-shaped domain Ω=(0,1.5)×(0,1)(0.5,1)×(-1,1), and we impose mixed boundary conditions given by

{𝒖(x,y)=(4y(1-y),0)Tif x=0,σnn(x,y)=0,ut(x,y)=0if x=1.5,𝒖(x,y)=(0,0)Totherwise.

In Figure 5 we plot the numerical solution obtained with the HDG discretisation using τ=6 on a coarse mesh. In this case, we used a mesh containing 379 402 triangles, which gives linear systems of a size 1 712 352 for the Taylor–Hood discretisation and 2 089 735 for the HDG discretisation. The initial guess in the GMRES iterative solver is zero. The overlapping decomposition into subdomains is generated by METIS and it has two layers of size h in the overlap.

Table 3

Preconditioners comparison – the T-shaped domain problem.


According to Table 3 the conclusions remain the same, that is, the standard RAS method performs far better when applied to a HDG discretisation with respect to a Taylor–Hood one and the MRAS preconditioners are better than the standard RAS preconditioner for both discretisations. Figure 6 again confirms the superiority of the MRAS preconditioner over the RAS preconditioner.

Figure 6

Convergence of error for METIS decomposition in the 800 subdomains case – the T-shaped domain problem.

(a) Taylor–Hood



(b) HDG



6 Conclusion

In this paper, we introduced a HDG method for the Stokes equations that naturally discretises nonstandard boundary value problems such as those with TVNF and NVTF boundary conditions. This approach can be extended naturally to the case of incompressible, or nearly incompressible, elasticity. We proved the well-posedness and convergence with respect to the norm (3.7) of this method and in the numerical experiments from Section 5.1 we validated the theory and observed the optimal convergence.

To solve the discretised problem we introduced two different kinds of preconditioners with nonstandard boundary conditions whose optimality has been proved by algebraic techniques. We compared the newly introduced preconditioners to the more standard RAS preconditioner and numerical tests from Section 5.2 clearly show their superiority for different test cases in two space dimensions. Moreover, our numerical experience seems to hint that the linear systems arising from the HDG discretisation are better conditioned than those obtained using the Taylor–Hood element. This can be seen from the fact that the RAS preconditioner already performs far better when applied to the HDG method than when using Taylor–Hood.

We observed, as expected, that the Schwarz preconditioners are not scalable with respect to the number of subdomains. However, this can be fixed by using an appropriate coarse spaces [12, Chapter 4]. A suitable choice of a coarse space will be a subject of future research.

Funding source: Centre for Numerical Analysis and Intelligent Software

Award Identifier / Grant number: EP/G036136/1

Funding statement: This research has been partially funded by the Engineering and Physical Sciences Research Council of Great Britain under the Centre for Numerical Algorithms and Intelligent Software (NAIS) for the evolving HPC platform grant EP/G036136/1.


We thank Frédéric Hecht from Laboratory J. L. Lions for comments that greatly improved the FreeFem++ codes.


[1] B. Ayuso de Dios, F. Brezzi, L. D. Marini, J. Xu and L. Zikatanov, A simple preconditioner for a discontinuous Galerkin method for the Stokes problem, J. Sci. Comput. 58 (2014), no. 3, 517–547. Search in Google Scholar

[2] D. Boffi, F. Brezzi and M. Fortin, Mixed Finite Element Methods and Applications, Springer Ser. Comput. Math. 44, Springer, Heidelberg, 2013. Search in Google Scholar

[3] X.-C. Cai and M. Sarkis, A restricted additive Schwarz preconditioner for general sparse linear systems, SIAM J. Sci. Comput. 21 (1999), no. 2, 792–797. Search in Google Scholar

[4] T. Cluzeau, V. Dolean, F. Nataf and A. Quadrat, Preconditionning techniques for systems of partial differential equations based on algebraic methods, Technical Report 7953, INRIA, 2012, Search in Google Scholar

[5] T. Cluzeau, V. Dolean, F. Nataf and A. Quadrat, Symbolic techniques for domain decomposition methods, Domain Decomposition Methods in Science and Engineering XX, Springer, Berlin (2013), 27–38. Search in Google Scholar

[6] B. Cockburn, D. A. Di Pietro and A. Ern, Bridging the hybrid high-order and hybridizable discontinuous Galerkin methods, ESAIM Math. Model. Numer. Anal. 50 (2016), no. 3, 635–650. Search in Google Scholar

[7] B. Cockburn and J. Gopalakrishnan, The derivation of hybridizable discontinuous Galerkin methods for Stokes flow, SIAM J. Numer. Anal. 47 (2009), no. 2, 1092–1125. Search in Google Scholar

[8] B. Cockburn, J. Gopalakrishnan and R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2009), no. 2, 1319–1365. Search in Google Scholar

[9] B. Cockburn, J. Gopalakrishnan, N. C. Nguyen, J. Peraire and F. J. Sayas, Analysis of HDG methods for Stokes flow, Math. Comp. 80 (2011), no. 274, 723–760. Search in Google Scholar

[10] D. A. Di Pietro and A. Ern, Mathematical Aspects of Discontinuous Galerkin Methods, Math. Appl. (Berlin) 69, Springer, Heidelberg, 2012. Search in Google Scholar

[11] D. A. Di Pietro and A. Ern, A hybrid high-order locking-free method for linear elasticity on general meshes, Comput. Methods Appl. Mech. Engrg. 283 (2015), 1–21. Search in Google Scholar

[12] V. Dolean, P. Jolivet and F. Nataf, An Introduction to Domain Decomposition Methods. Algorithms, Theory, and Parallel Implementation, Society for Industrial and Applied Mathematics, Philadelphia, 2015. Search in Google Scholar

[13] V. Dolean and F. Nataf, A new domain decomposition method for the compressible Euler equations, ESAIM Math. Model. Numer. Anal. 40 (2006), no. 4, 689–703. Search in Google Scholar

[14] V. Dolean, F. Nataf and F. Rapin, Deriving a new domain decomposition method for the Stokes equations using the Smith factorization, Math. Comp. 78 (2009), no. 266, 789–814. Search in Google Scholar

[15] H. Egger and C. Waluga, A hybrid discontinuous Galerkin method for Darcy–Stokes problems, Domain Decomposition Methods in Science and Engineering XX, Lect. Notes Comput. Sci. Eng. 91, Springer, Berlin (2013), 663–670. Search in Google Scholar

[16] H. Egger and C. Waluga, hp analysis of a hybrid DG method for Stokes flow, IMA J. Numer. Anal. 33 (2013), no. 2, 687–721. Search in Google Scholar

[17] A. Ern and J. L. Guermond, Theory and Practice of Finite Elements, Appl. Math. Sci. 159, Springer, New York, 2004. Search in Google Scholar

[18] V. Girault and P. A. Raviart, Finite Element Methods for Navier–Stokes Equations. Theory and Algorithms, Springer Ser. Comput. Math. 5, Springer, Berlin, 1986. Search in Google Scholar

[19] P. Gosselet and C. Rey, Non-overlapping domain decomposition methods in structural mechanics, Arch. Comput. Methods Eng. 13 (2006), no. 4, 515–572. Search in Google Scholar

[20] F. Hecht, New development in FreeFem++, J. Numer. Math. 20 (2012), no. 3–4, 251–265. Search in Google Scholar

[21] C. Lehrenfeld, Hybrid discontinuous Galerkin methods for solving incompressible flow problems, Dissertation, Rheinisch-Westfälischen Technischen Hochschule Aachen, 2010. Search in Google Scholar

[22] C. Lehrenfeld and J. Schöberl, High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows, Comput. Methods Appl. Mech. Engrg. 307 (2016), 339–361. Search in Google Scholar

[23] I. Oikawa, Analysis of a reduced-order HDG method for the Stokes equations, J. Sci. Comput. 67 (2016), no. 2, 475–492. Search in Google Scholar

[24] A. Quarteroni and A. Valli, Domain Decomposition Methods for Partial Differential Equations, Clarendon Press, Oxford, 1999. Search in Google Scholar

[25] W. H. Reed and T. R. Hill, Triangular mesh methods for the neutron transport for a scalar hyperbolic equation, Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973. Search in Google Scholar

[26] Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Statist. Comput. 7 (1986), no. 3, 856–869. Search in Google Scholar

[27] D. Schötzau, C. Schwab and A. Toselli, Mixed hp-DGFEM for incompressible flows, SIAM J. Numer. Anal. 40 (2002), no. 6, 2171–2194. Search in Google Scholar

[28] B. F. Smith, P. E. Bjørstad and W. Gropp, Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press, Cambridge, 1996. Search in Google Scholar

[29] H. J. S. Smith, On systems of linear indeterminate equations and congruences, Philos. Trans. A 151 (1861), 293–326. Search in Google Scholar

[30] A. St-Cyr, M. J. Gander and S. J. Thomas, Optimized multiplicative, additive, and restricted additive Schwarz preconditioning, SIAM J. Sci. Comput. 29 (2007), no. 6, 2402–2425. Search in Google Scholar

[31] A. Toselli and O. Widlund, Domain Decomposition Methods – Algorithms and Theory, Springer Ser. Comput. Math. 34, Springer, Berlin, 2005. Search in Google Scholar

Received: 2017-09-06
Revised: 2018-02-02
Accepted: 2018-02-28
Published Online: 2018-03-22
Published in Print: 2019-10-01

© 2020 Walter de Gruyter GmbH, Berlin/Boston