 # Numerical assessment of two-level domain decomposition preconditioners for incompressible Stokes and elasticity equations

Solving the linear elasticity and Stokes equations by an optimal domain decomposition method derived algebraically involves the use of non standard interface conditions. The one-level domain decomposition preconditioners are based on the solution of local problems. This has the undesired consequence that the results are not scalable, it means that the number of iterations needed to reach convergence increases with the number of subdomains. This is the reason why in this work we introduce, and test numerically, two-level preconditioners. Such preconditioners use a coarse space in their construction. We consider the nearly incompressible elasticity problems and Stokes equations, and discretise them by using two finite element methods, namely, the hybrid discontinuous Galerkin and Taylor-Hood discretisations.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In [BBD16] the one-level domain decomposition methods for Stokes equations were introduced in conjunction with the non standard interface conditions. Although it can be observed there the lack of scalability with respect to the number of subdomains. It means that by splitting the problem in a larger number of subdomains leads to the increase of size of the plateau region in the convergence of an iterative method (see Figure 3) when using the one-level domain decomposition methods. This is caused by the lack of global information, as subdomains can only communicate with their neighbours. Hence, when the number of subdomains increases in one direction, the length of the plateau also increases. Even in cases when the local problems are of the same size, the iteration count grows with the increase of the number of subdomain. This can be also observed in all experiments in this manuscript in case of one-level methods.

The remedy for this is the use of a second level in the preconditioner or a coarse space correction that adds the necessary global information. Two-level algorithms have been analysed for several classes of problems in [TW05]. The key point of these kind of methods is to choose the appropriate coarse space. The classical coarse space introduced by Nicolaides in [Nic87]

is defined by vectors that support is in each subdomain. Hence the coarse space has the size equal to a number of subdomains. A coarse space construction, named spectral coarse space, was motivated by the complexity of the problems that classical coarse space performance was not satisfying. This construction allows to enrich a bigger size of the coarse space, but can be also reduced to the classical one. This idea was introduced for the first time in

[BHMV99]

in the case of multigrid methods. It relies on solving local generalised eigenvalue problems allowing to choose suitable vectors for the coarse space.

For overlapping domain decomposition preconditioners, a similar idea was introduced in the case of Darcy equations in [GE10a, GE10b]. The authors of [NXD10] consider also the heterogeneous Darcy equation and presented a different generalised eigenvalue problem based on local Dirichlet-to-Neumann maps. The method has been analysed in [DNSS12] and proved to be very robust in the case of small overlaps. The same idea was extended numerically to the heterogeneous Helmholtz problem in [CDKN14]. The authors of [LNS15]

apply the coarse space associated with low-frequency eigenfunctions of the subdomain Dirichlet-to-Neumann maps for the generalisation of the optimised Schwarz methods, named 2-Lagrange multiplier methods.

The first attempt to extend this spectral approach to general symmetric positive definite problems was made in [EGLW12] as an extension of [GE10a, GE10b]. Since some of the assumptions of the previous framework are hard to fulfil, authors of [SDH14]

proposed slightly different approach for symmetric positive definite problems. Their idea of constructing partition of unity operator associated with degrees of freedom allows to work with various finite element spaces. An overview of different kinds of two-level methods can be found in

[DJN15, Chapters 5 and 7].

Despite the fact that all these approaches provide satisfying results, there is no universal treatment to build efficient coarse spaces in the case of non definite problems such as Stokes equations. The spectral coarse spaces that we use in this work are inspired by those proposed in [HJN15]. The authors introduced and tested numerically symmetrised two-level preconditioners for overlapping algorithms which use Robin interface conditions between the subdomains (see (5.27) for details). They have applied these preconditioners to solve saddle point problems such as nearly incompressible elasticity and Stokes discretised by Taylor-Hood finite elements. In our case, we use non standard interface conditions. Therefore the use of spectral coarse spaces could lead to an important gain.

In this work, we test this improvement in case of nearly incompressible elasticity and Stokes equations that are discussed in Section 2. As the discretisations we use the Taylor-Hood [GR86, Chapter II, Section 4.2] and hybrid discontinuous Galerkin method [CGL09, CG09] presented in Section 3. In Section 4 we introduce the two-level domain decomposition preconditioners. Sections 5 and 6 present the two and three dimensional numerical experiments, respectively. Finally, a summary is outlined in Section 7.

## 2 The differential equations

Let be an open polygon in or an open Lipschitz polyhedron in , with Lipschitz boundary . We use

to denote the dimension of the space. We use bold for tensor or vector variables. In addition we denote normal and tangential components as follows

, and , where is the outward unit normal vector to the boundary .

For , we use the standard space and denotes the set of all continuous functions on the closure of a set . Let us define following Sobolev spaces

 Hm(D) :={v∈L2(D): ∀ |α|≤m ∂αv∈L2(D)} for m∈N, H12(∂D) :={~v∈L2(∂D): ∃v∈H1(D)~v=tr(v)}, H(div,D) :={v∈[L2(D)]d: ∇⋅v∈L2(D)},

where, for and we denote , and is the trace operator. In addition, we use the following notation of the space including boundary and average conditions

 L20(D) :={v∈L2(D): ∫Dv dx=0}, H1~Γ(D) :={v∈H1(D): v=0 on ~Γ},

where . If , then is denoted .

Now we present the two differential problems considered in this work.

### 2.1 Stokes equation

 (2.1) {−νΔu +∇p=fin Ω∇⋅u=0in Ω,

where is the velocity field, the pressure, the viscosity which is considered to be constant and is a given function. We define the stress tensor and the flux as . For and we consider three types of boundary conditions:

• Dirichlet (non-slip)

 (2.2) u=uD on Γ;
• tangential-velocity and normal-flux (TVNF)

 (2.3) {ut=0 on Γσnn=g on Γ;
• normal-velocity and tangential-flux (NVTF)

 (2.4) {un=0 on Γσnt=g on Γ.

The third type of boundary condition has already been considered for the Stokes problem in [AdDBM14].

### 2.2 Nearly incompressible elasticity equation

From a mathematical point of view, the nearly incompressible elasticity problem is very similar to the Stokes equations. The difference is that instead of considering the gradient , the symmetric gradient is used. We want to solve the following -dimensional, , problem

 (2.5) {−2μ∇⋅ε(u)+∇p=fin Ω−∇⋅u=1λpin Ω

where is the displacement field, the pressure, is a given function, and are the Lamé coefficients defined by

 λ=Eν(1+ν)(1−2ν), μ=E2(1+ν),

where is the Young modulus and the Poisson ratio. We define the stress tensor as and its normal component as . For we consider three types of boundary conditions:

• mixed such that and

 (2.6) {u=0 on ΓDσsymn=0 on ΓN;
• tangential-displacement and normal-normal-stress (TDNNS)

 (2.7) {ut=0 on Γσsymnn=g on Γ;
• normal-displacement and tangential-normal-stress (NDTNS)

 (2.8) {un=0 on Γσsymnt=g on Γ.

The second type of boundary condition has already been considered for linear elasticity equation in [PS11].

## 3 The numerical methods

Let be a regular family of triangulations of made of simplices. For each triangulation , denotes the set of its facets (edges for , faces for ). In addition, for each element , , and we denote . We define the following broken Sobolev spaces on the set of all edges in (for )

 L2(Eh) :={v: v|E∈L2(E) ∀ E∈Eh}.

Moreover, for , denotes the space of polynomials of total degree smaller than, or equal to, on the set .

We now present the two discretisations to be used in the numerical experiments.

### 3.1 Taylor-Hood discretisation

We first consider the Taylor-Hood discretisation using the following approximation spaces

 THkh :={vh∈[H1(Ω)]d:vh|K∈[Pk(K)]d∀ K∈Th}, Rk−1h :={qh∈C0(¯Ω):qh|K∈Pk−1(K)∀ K∈Th}

where (see [GR86, Chapter II, Section 4.2]).

If (2.1) is supplied with the homogeneous boundary conditions (2.2), then the discrete problem reads:

Find

s.t. for all

 (3.9) ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∫Ων∇uh:∇vh dx−∫Ωph∇⋅vh dx=∫Ωfvh dx−∫Ω∇⋅uhqh dx=0.

In case of TVNF boundary conditions (2.3), we define , and the discrete problem reads:

Find

s.t. for all

 (3.10) ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∫Ων∇uh:∇vh dx−∫Ωph∇⋅vh dx=∫Ωfvh dx+∫Γg(vh)t ds−∫Ω∇⋅uhqh dx=0.

If NVTF boundary conditions (2.4) are used, then we define the following space
, and the discrete problem reads:

Find

s.t. for all

 (3.11) ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∫Ων∇uh:∇vh dx−∫Ωph∇⋅vh dx=∫Ωfvh dx+∫Γg(vh)n ds−∫Ω∇⋅uhqh dx=0.

In similar way, if the problem (2.5) is supplied with the boundary conditions (2.6), then the discrete problem reads

Find

s.t. for all

 (3.12) ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∫Ω2με(uh):ε(vh) dx−∫Ωph∇⋅vh dx=∫Ωfvh dx−∫Ω∇⋅uhqh dx−1λ∫Ωphqh dx=0.

The rest of the discrete problems associated with (2.5) that is supplied with TDNNS boundary conditions (2.7) or NDTNS boundary conditions (2.8) are similar to (3.10) or (3.11), respectively.

### 3.2 Hybrid discontinuous Galerkin discretisation

We restrict the discussion of this to two dimensional case . This method has been presented and analysed in [BBD16]. The velocity is approximated using the Brezzi-Douglas-Marini spaces (see [BBF13, Section 2.3.1]) of degree given by

 BDMkh :={vh∈H(div,Ω): vh|K∈[Pk(K)]2 ∀ K∈Th}, BDMkh,~Γ :={vh∈H(div,Ω): vh|K∈[Pk(K)]2 ∀ K∈Th∧ (vh)n=0 on ~Γ},

where . If , then is denoted .

The pressure is approximated in the space

 Qk−1h

Finally, a Lagrange multiplier, aimed at approximating the tangential component of the velocity is introduced. The space swhere this multipliers is sought are given by

 Mk−1h :={~vh∈L2(Eh): ~vh|E∈Pk−1(E) ∀ E∈Eh}, Mk−1h,~Γ :={~vh∈Mk−1h:~vh=0 on ~Γ},

where . If , then is denoted . Furthermore, we introduce for all the -projection defined by

 (3.13) ∫EΦk−1E(~w)~vh ds=∫E~w~vh ds∀ ~vh∈Pk−1(E),

and we denote defined as for all .

If (2.1) is supplied with the homogeneous boundary conditions (2.2), then the discrete problem reads:

Find

s.t. for all ,

 (3.14)

where

 a((wh,~wh),(vh,~vh)):= ∑K∈Th(∫Kν∇wh:∇vh dx (3.15) −∫∂Kν((wh)t−~wh)(∂nvh)t ds

and is a stabilisation parameter, and

 (3.16) b((vh,~vh),qh):=−∑K∈Th∫Kqh∇⋅vh dx.

If TVNF boundary conditions (2.3) are used, then the discrete problem reads:

Find

s.t. for all ,

 (3.17)

In case of NVTF boundary conditions (2.4), the discrete problem reads:

Find

s.t. for all ,

 (3.18)

In similar way, if the problem (2.5) is supplied with the mixed boundary conditions (2.6), then the discrete problem reads:

Find

s.t. for all

 (3.19) ⎧⎪⎨⎪⎩as((uh,~uh),(vh,~vh))+b((vh,~vh),ph)=∫Ωfvh dxb((uh,~uh),qh)+c(ph,qh)=0,

where

 ∑K∈Th(∫K2με(wh):ε(vh) dx −∫∂K2μ(εn(wh))t((vh)t−~vh) ds (3.20) −∫∂K2μ((wh)t−~wh)(εn(vh))t ds

is defined by (3.16), and

 c(rh,qh):=−1λ∫Ωrhqh ds.

The rest of the discrete problems associated with (2.5) that is supplied with TDNNS boundary conditions (2.7) or NDTNS boundary conditions (2.8) are similar to (3.17) or (3.18), respectively.

## 4 The domain decomposition preconditioners

Let us assume that we have to solve the following linear system where is the matrix arising from discretisation of the Stokes or linear elasticity equation on the domain , is the vector of unknowns, and is the right hand side. To accelerate the performance of an iterative Krylov method [DJN15, Chapter 3] applied to this system we will consider domain decomposition preconditioners which are naturally parallel. They are based on an overlapping decomposition of the computational domain.

Let be a partition of the triangulation (see examples in Figure 6). For an integer value , we define an overlapping decomposition such that is a set of all triangles from and all triangles from that have non-empty intersection with , and . With this definition the width of the overlap will be of . Furthermore, if stands for the finite element space associated to , is the local finite element spaces on , which is a triangulation of .

Let be the set of indices of degrees of freedom of and the set of indices of degrees of freedom of for . Moreover, we define the restriction operator as a rectangular matrix such that if is the vector of degrees of freedom of , then is the vector of degrees of freedom of in . The extension operator from to and its associated matrix are both given by . In addition we introduce a partition of unity as a diagonal matrix such that

 (4.21) Id=N∑i=1RTiDiRi,

where

is the identity matrix.

We first recall the Modified Restricted Additive Schwarz (MRAS) preconditioner introduced in [BBD16] for the Stokes equation. This preconditioner is given by

 (4.22) M−1MRAS=N∑i=1RTiDiB−1iRi,

where is the matrix associated to a discretisation of Stokes equation (2.1) in where we impose either TVNF (2.3) or NVTF (2.4) boundary conditions in . In case of a discretisation of elasticity equation (2.5) in , we impose either TDNNS (2.7) or NDTNS (2.8) boundary conditions in .

We new introduce a symmetrised variant of (4.22) called Symmetrised Modified Restricted Additive Schwarz (SMRAS), given by

 (4.23) M−1SMRAS=N∑i=1RTiDiB−1iDiRi.

### 4.1 Two-level methods

A two-level version of the SMRAS and MRAS preconditioners will be based on a spectral coarse space obtained by solving the following local generalised eigenvalue problems

Find s.t.

 (4.24) ~AjVjk=λjkBjVjk,

where are local matrices associated to a discretisation of local Neumann boundary value problem in . Let be a user-defined threshold. We define as the vector space spanned by the family of vectors , , corresponding to eigenvalues smaller than . The value of is chosen such that for a given problem and preconditioner, the behaviour of the method should be robust in the sense that, its convergence should not depend, or depends very weakly, on the number of subdomains.

We are now ready to introduce the two-level method with coarse space . Let be the -orthogonal projection onto the coarse space . The two-level SMRAS preconditioner is defined as

 (4.25) M−1SMRAS,2=P0A−1+(Id−P0)M−1SMRAS(Id−PT0).

Furthermore, if is a matrix whose rows are a basis of the coarse space , then

 P0A−1=RT0(R0ART0)−1R0.

In similar way, we can introduce the two-level MRAS preconditioner

 (4.26) M−1MRAS,2=P0A−1+(Id−P0)M−1MRAS(Id−PT0).

## 5 Numerical results for two dimensional problems

In this section we assess the performance of the preconditioners defined in Section 4.1. We will compare the newly introduced ones with that of ORAS and SORAS introduced in [HJN15]. These kind of preconditioners are associated with the Robin interface conditions and require an optimised parameter as it can be seen in (5.27

) below. The big advantage of SMRAS and MRAS preconditioners from the previous section is that they are parameter-free. We consider the partial differential equation model for nearly incompressible elasticity and Stokes flow as problems of similar mixed formulation. Each of these problems is discretised by using the Taylor-Hood methods from Section

3.1 and the hdG discretisation from Section 3.2.

Our experiments will be based on the classical weak scaling test. This test is built as follows. A domain is split into a triangulation . For each of element , , and we denote the mesh size by . Then, this triangulation is split into overlapping subdomains of size H, in such a way remains constant. In the absence of a second level in the preconditioner, if the number of subdomains grows, then the convergence gets slower. A coarse space provides a global information and leads to a more robust behaviour.

The simplest way to build a coarse space is to consider the zero energy modes. More precisely they are the eigenvectors associated with the zero eigenvalues of (

4.24) on a floating subdomains. Hence, by a floating subdomain we mean a subdomain without Dirichlet boundary condition on any part of the boundary. Then the matrix on the left hand side of (4.24) is singular and there are zero eigenvalues. These zero energy modes are the rigid body motions (three in two dimensions, six in three dimensions) for the elasticity problem, and the constants (two in two dimensions, three in three dimensions) for the Stokes equations. Unfortunately for some cases, this choice is not sufficient, so we have collected the smallest eigenvalues for each subdomain and build a coarse space by including the eigenvectors associated to them. The different values of are presented in the table in brackets.

All experiments have been made by using FreeFem++ [Hec12], which is a free software specialised in variational discretisations of partial differential equations. We use GMRES [SS86] as an iterative solver. Generalized eigenvalue problems to generate the coarse space are solved using ARPACK [LSY98]. The overlapping decomposition into subdomains can be uniform (Unif) or generated by METIS (MTS) [KK98]. In each of the examples in this section we consider decomposition with two layers of mesh size in the overlap. Tables show the number of iterations needed to achieve a relative norm of the error smaller than , , where is the solution of global problem given by direct solver and denotes -th iteration of the iterative solver. In addition, stands for number of degrees of freedom and for the number of subdomains in all tables.

### 5.1 Taylor-Hood discretisation

In this section we consider the Taylor-Hood discretisation from Section 3.1 with different values of for nearly incompressible elasticity and Stokes equations.

#### 5.1.1 Nearly incompressible elasticity

Since we consider the preconditioners with various interface conditions we need to comment the way of imposing them. ORAS and SORAS preconditioners follow [HJN15] and use Robin interface conditions. This means, the weak formulation of the linear elasticity problem contains the following term

 (5.27) ∫∂Ωi∖Γσsymn(vh)n ds+∫∂Ωi∖Γ2αμ(2μ+λ)λ+3μuhvh ds

where again, following [HJN15] we choose . Fortunately, the MRAS and SMRAS preconditioners are parameter-free. For all numerical experiments associated in this section we use zero as an initial guess for the GMRES iterative solver. Moreover, the overlapping decomposition into subdomains is generated by METIS.

###### Test case 1 (The L-shaped domain problem).

We consider the L-shaped domain clamped on the left side and partly from a top and bottom as it is depicted in Figure (a)a. This example is similar to the one in [CS15]. The associated boundary value problem is

 (5.28) ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩−2μ∇⋅ε(u)+∇p=(0,−1)Tin Ω−∇⋅u=1λpin Ωu(x,y)=(0,0)Ton ΓDσsymn(x,y)=(0,0)Ton ΓN.

The physical parameters are and (nearly incompressible). In Figure (b)b we plot the mesh of the bent domain.

We choose for the Taylor-Hood discretisation. In Figure 13 we plot the eigenvalues of one floating subdomain. The clustering of small eigenvalues of the generalised eigenvalue problem defined in (4.24) suggests the number of eigenvectors to be added to the coarse space. The three zero eigenvalues correspond to the zero energy modes.

The results of Table 1 show a clear improvement in the scalability of the two-level preconditioners over the one-level ones. In fact, using five eigenvectors per subdomain, the number of iterations is virtually unaffected by the number of subdomains. All two-level preconditioners show a comparable performance. For this case, increasing the dimension of the coarse space beyond eigenvectors does not seem to improve the results dramatically.

###### Test case 2 (The heterogeneous beam problem).

We consider a heterogeneous beam with ten layers of steel and rubber. Five layers are made from steel with the physical parameters and , and other five are made from rubber with the physical parameters and as it is depicted in Figure (a)a. A similar example was considered in [HJN15].

The computational domain is the rectangle . The beam is clamped on its left side, hence we consider the following problem

 (5.29) ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩−2μ∇⋅ε(u)+∇p=(0,−1)Tin Ω−∇⋅u=1λpin Ωu(x,y)=(0,0)Ton ∂Ω∩{x=0}σsymn(x,y)=(0,0)Ton ∂Ω∖{x=0}.

In Figure (b)b we plot the mesh of the bent beam. Because of the heterogeneity of the problem, we do not notice a clear clustering of the eigenvalues (see Figure 20). In such case it is well known that coarse space including only three zero energy modes is not sufficient [DNSS12]. That is why we consider a coarse space built using 5 or 7 eigenvectors per subdomain.

As in the previous example, the introduction of a coarse space provides a significant improvement in the number of iterations needed for convergence. Due to the high heterogeneity of this problem, more eigenvectors per subdomain are needed to obtain scalable results. We notice an important improvement of the convergence when using two-level methods (see Table 2). Although we get a stable number of iterations only when considering a coarse space which is sufficiently big.

#### 5.1.2 Stokes equation

We now turn to the Stokes discrete problem given in Sections 3.1. Once again in case of ORAS and SORAS we choose as in [HJN15] for the Robin interface conditions (5.27). In the first case we consider a random initial guess for the GMRES iterative solver. Later with the second example we use zero as an initial guess.