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

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:=𝒖⋅𝒕, σn⁢n:=𝝈𝒏⋅𝒏, σn⁢t:=𝝈𝒏⋅𝒕, 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:

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

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):

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

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 g∈L2⁢(Γ) we consider two types of boundary conditions:

1. tangential-velocity and normal-flux (TVNF)

   (2.2){σn⁢n=gon ⁢Γ,ut=0on ⁢Γ,

2. normal-velocity and tangential-flux (NVTF)

   (2.3){σn⁢t=gon ⁢Γ,un=0on ⁢Γ,

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

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∈𝒯h⁡hK. We define the following Sobolev spaces on the triangulation 𝒯h and the set of all edges in ℰh:

with the corresponding broken norms.

The following results will be very useful in what follows.

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 k≥1, we discretise the velocity using the Brezzi–Douglas–Marini space (see [2, Section 2.3.1]) given by

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 1≤m≤k+1 we denote by Πk:[Hm⁢(Ω)]2→BDMhk 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

Furthermore, we introduce for all E∈ℰh 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

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

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

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 w∈L2⁢(K), ΨKk⁢(w) is the unique element of ℙk-1⁢(K) satisfying

We will also use the following results.

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.