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

06/29/2017 ∙ by Gabriel R. Barrenechea, et al. ∙ 0

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.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 8

page 10

page 13

page 15

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.

(a) Taylor-Hood
(b) hdG
Figure 3: Increase of the size of the plateau region for increasing number of subdomains.

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

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

where . If , then is denoted .

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

2.1 Stokes equation

Let us start with -dimensional, , Stokes problem

(2.1)

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)
  • tangential-velocity and normal-flux (TVNF)

    (2.3)
  • normal-velocity and tangential-flux (NVTF)

    (2.4)

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)

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

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)
  • tangential-displacement and normal-normal-stress (TDNNS)

    (2.7)
  • normal-displacement and tangential-normal-stress (NDTNS)

    (2.8)

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 )

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

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)

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

Find

s.t. for all

(3.10)

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)

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)

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

where . If , then is denoted .

The pressure is approximated in the space

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

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

(3.13)

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

(3.15)

and is a stabilisation parameter, and

(3.16)

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)

where

(3.20)

is defined by (3.16), and

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.

(a) L-shaped domain
(b) T-shaped domain
Figure 6: Partition of the domain for 8 subdomains by METIS.

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)

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)

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)

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)

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)

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

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

(4.26)

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)

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.

(a) Domain and boundary
(b) Discrete solution
Figure 9: L-shaped domain problem.
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)

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.

max width= One-level DOF N ORAS SORAS NDTNS-MRAS NDTNS-SMRAS TDNNS-MRAS TDNNS-SMRAS 124 109 4 26 60 26 60 30 59 478 027 16 57 131 69 143 65 140 933 087 32 84 180 109 221 104 211 1 899 125 64 130 293 181 362 161 312 3 750 823 128 209 412 302 568 251 510 Two-level (3 eigenvectors) DOF N ORAS SORAS NDTNS-MRAS NDTNS-SMRAS TDNNS-MRAS TDNNS-SMRAS 124 109 4 18 40 19 36 24 41 478 027 16 37 52 40 57 46 56 933 087 32 49 57 56 67 53 66 1 899 125 64 65 64 70 75 61 74 3 750 823 128 83 64 74 77 75 72 Two-level (5 eigenvectors) DOF N ORAS SORAS NDTNS-MRAS NDTNS-SMRAS TDNNS-MRAS TDNNS-SMRAS 124 109 4 15 32 17 35 24 37 478 027 16 31 41 31 47 42 47 933 087 32 40 48 38 52 53 51 1 899 125 64 49 51 45 53 64 56 3 750 823 128 69 54 49 54 70 53 Two-level (7 eigenvectors) DOF N ORAS SORAS NDTNS-MRAS NDTNS-SMRAS TDNNS-MRAS TDNNS-SMRAS 124 109 4 14 33 16 30 24 35 478 027 16 26 41 25 38 42 44 933 087 32 31 43 25 42 49 46 1 899 125 64 39 47 30 39 59 50 3 750 823 128 58 49 30 43 61 50

Table 1: Comparison of preconditioners for Taylor-Hood discretisation () - the L-shaped domain problem.
(a) Robin interface conditions
(b) NDTNS interface conditions
(c) TDNNS interface conditions
Figure 13: Eigenvalues on one of the floating subdomains in case of uniform decomposition and Taylor-Hood discretisation () - the L-shaped domain problem.

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].

(a) Steel and rubber layers
(b) Discrete solution
Figure 16: Heterogeneous beam.

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

(5.29)

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.

max width= One-level DOF N ORAS SORAS NDTNS-MRAS NDTNS-SMRAS TDNNS-MRAS TDNNS-SMRAS 44 963 8 168 301 160 267 177 264 87 587 16 226 490 245 462 229 424 177 923 32 373 711 447 684 440 672 347 651 64 615 1000 728 1000 746 1000 707 843 128 973 1000 1000 1000 1000 1000 1 385 219 256 1000 1000 1000 1000 1000 1000 Two-level (5 eigenvectors) DOF N ORAS SORAS NDTNS-MRAS NDTNS-SMRAS TDNNS-MRAS TDNNS-SMRAS 44 963 8 109 160 136 147 148 136 87 587 16 136 204 192 200 181 184 177 923 32 193 291 296 275 326 276 347 651 64 260 304 363 282 491 299 707 843 128 412 356 420 369 601 346 1 385 219 256 379 414 448 400 711 317 Two-level (7 eigenvectors) DOF N ORAS SORAS NDTNS-MRAS NDTNS-SMRAS TDNNS-MRAS TDNNS-SMRAS 44 963 8 76 118 124 115 133 103 87 587 16 106 146 166 138 159 123 177 923 32 157 202 203 185 302 214 347 651 64 178 191 225 170 326 182 707 843 128 140 114 153 112 266 122 1 385 219 256 119 86 118 77 259 94

Table 2: Comparison of preconditioners for Taylor-Hood discretisation () - the heterogeneous beam.
(a) Robin interface conditions
(b) NDTNS interface conditions
(c) TDNNS interface conditions
Figure 20: Eigenvalues on one of the floating subdomains in case of METIS decomposition and Taylor-Hood discretisation () - the heterogeneous beam.

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.

(a) Velocity field
(b) Pressure
Figure 23: Numerical solution of the driven cavity problem - the driven cavity problem.
Test case 3 (The driven cavity problem).

The test case is the driven cavity. We consider the following problem on the unit square

(5.30)

In Figure 23 we plot the vector field and pressure, after solving numerically the problem.

We start with the two energy modes only (see Figure 27). This already provides some improvement. Then, we add more eigenvectors to see if they bring improvement.

(a) Robin interface conditions
(b) NVTF interface conditions
(c) TVNF interface conditions
Figure 27: Eigenvalues on one of the floating subdomains in case of uniform decomposition and Taylor-Hood discretisation () - the driven cavity problem.

max width= One-level DOF N ORAS SORAS NVTF-MRAS NVTF-SMRAS TVNF-MRAS TVNF-SMRAS Unif MTS Unif MTS Unif MTS Unif MTS Unif MTS Unif MTS 91 003 4 12 17 24 34 22 22 34 40 22 25 30 40 362 003 16 28 35 56 67 52 53 90 106 54 53 70 84 813 003 36 39 75 92 103 85 91 165 185 91 88 118 136 1 444 003 64 53 91 120 144 120 135 254 283 132 132 169 206 2 728 003 121 80 278 180 212 182 280 412 580 199 213 251 439 5 768 003 256 1000 1000 271 317 303 452 917 955 322 319 397 695 Two-level (2 eigenvectors) DOF N ORAS SORAS NVTF-MRAS NVTF-SMRAS TVNF-MRAS TVNF-SMRAS Unif MTS Unif MTS Unif MTS Unif MTS Unif MTS Unif MTS 91 003 4 10 14 18 22 19 17 26 30 27 20 21 26 362 003 16 20 25 32 37 33 34 50 62 60 40 42 51 813 003 36 27 33 36 44 47 49 62 86 79 53 59 63 1 444 003 64 31 42 38 53 104 66 85 114 85 52 62 79 2 728 003 121 39 103 39 51 74 81 85 133 92 86 62 93 5 768 003 256 300 849 46 54 109 108 146 132 91 78 63 90 Two-level (5 eigenvectors) DOF N ORAS SORAS NVTF-MRAS NVTF-SMRAS TVNF-MRAS TVNF-SMRAS Unif MTS Unif MTS Unif MTS Unif MTS Unif MTS Unif MTS 91 003 4 9 12 13 16 16 15 18 20 25 20 16 18 362 003 16 16 20 21 24 27 22 28 37 56 37 26 35 813 003 36 23 27 25 26 33 30 39 40 65 41 28 37 1 444 003 64 26 36 27 29 40 34 35 45 77 45 28 42 2 728 003 121 35 41 29 32 43 38 34 48 84 72 29 47 5 768 003 256 66 60