Robust block preconditioners for poroelasticity

In this paper we study the linear systems arising from discretized poroelasticity problems. We formulate one block preconditioner for the two-filed Biot model and several preconditioners for the classical three-filed Biot model under the unified relationship framework between well-posedness and preconditioners. By the unified theory, we show all the considered preconditioners are uniformly optimal with respect to material and discretization parameters. Numerical tests demonstrate the robustness of these preconditioners.



page 1

page 2

page 3

page 4


A nonsymmetric approach and a quasi-optimal and robust discretization for the Biot's model. Part I – Theoretical aspects

We consider the system of partial differential equations stemming from t...

Block preconditioning methods for asymptotic preserving scheme arising in anisotropic elliptic problems

Efficient and robust iterative solvers for strong anisotropic elliptic e...

Algebraic multigrid block preconditioning for multi-group radiation diffusion equations

The paper focuses on developing and studying efficient block preconditio...

Optimal Statistical Hypothesis Testing for Social Choice

We address the following question in this paper: "What are the most robu...

Tug-of-War: Observations on Unified Content Handling

Modern applications and Operating Systems vary greatly with respect to h...

A note on the implementation of the hyperelastic Kilian model in the Abaqus

The main goal of the article is to verify the implementation of the Kili...
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

Poroelasticity, the study of the fluid flow in porous and elastic media, couples the elastic deformation with the fluid flow in porous media. The Biot model has wide applications in geoscience, biomechanics, and many other fields. Numerical simulations of poroelasticity are challenging. A notorious instability of the numerical discretization method for the poroelasticity model is the unphysical oscillation of the pressure under certain conditions Murad.M;Loula.A1994a . There are many possible sources of this instability. One of the most significant sources is the instability of the finite element approximation for the coupled systems Haga.J;Osnes.H;Langtangen.H2012b ; Axelsson.O;Blaheta.R;Byczanski.P2012a . This motivates us to study the well-posedness of the finite element discretization. However, we do not look further into the details of the instability of the Biot model and refer interested readers to Haga.J;Osnes.H;Langtangen.H2012b ; Axelsson.O;Blaheta.R;Byczanski.P2012a ; Aguilar.G;Gaspar.F;Lisbona.F;Rodrigo.C2008a .

Another challenge associated with the Biot model is that of developing efficient linear solvers. Direct solvers have poor performance when the size of problems becomes large. Iterative solvers are good alternatives, as they exhibit better scalability. However, the convergence of iterative solvers is very much problem-dependent such that there is a need for robust preconditioners. For example, the multigrid preconditioned Krylov subspace method usually has optimal convergence rate for the Poission equation and many other symmetric positive definite problems Hackbusch.W1985a ; Xu.J1992a . However, for poroelasticity problems, coupled systems of equations must be solved, which are known to be indefinite and ill-conditioned Ferronato.M;Gambolati.G;Teatini.P2001a . Preconditioning techniques for poroelasticity problems have been the subject of considerable research in the literature Axelsson.O;Blaheta.R;Byczanski.P2012a ; Lipnikov.K2002a ; Bergamaschi.L;Ferronato.M;Gambolati.G2007a ; Phoon.K;Toh.K;Chan.S;Lee.F2002a ; Toh.K;Phoon.K;Chan.S2004a ; Haga.J;Osnes.H;Langtangen.H2011a ; Haga.J;Osnes.H;Langtangen.H2012a ; Turan.E;Arbenz.P2014a and most of the techniques developed are based on the Schur complement approach. In Phoon.K;Toh.K;Chan.S;Lee.F2002a ; Toh.K;Phoon.K;Chan.S2004a , diagonal approximation of the Schur complement preconditioner is used to precondition two-field formulation of the Biot model. In Haga.J;Osnes.H;Langtangen.H2011a ; Haga.J;Osnes.H;Langtangen.H2012a , Schur complement preconditioners are also studied for two-field formulation with the algebraic multigrid (AMG) as the preconditioner for the elasticity block. In Axelsson.O;Blaheta.R;Byczanski.P2012a , Schur complement approaches for three-field formulation are investigated. Recently, robust block diagonal and block triangular preconditioners are developed in adler2017robust for two-field Biot model. And for classical three-filed Biot model, the robust block preconditioners are designed in hong2018parameter ; hong2018conservative

based on the uniform stability estimates. Robust preconditioner for a new three-field formulation introducing a total pressure as the third unknown is analyzed in 

lee2017parameter . Robust block diagonal and block triangular preconditioners are also developed in adler2019robust based on the descretization proposed in rodrigo2018new . Other robustness analysis for fixed-stress splitting method and Uzawa-type method for multiple-permeability poroelasticity systems are presented in hong2018fixed and hong2019parameter .

The focus of this paper is on the stability of the linear systems after time discretization and several robust preconditioners for the iterative solvers under the unified relationship framework between well-posedness and preconditioners. The block preconditioners in adler2017robust for two field formulation and in hong2018parameter ; adler2019robust for the three field formulation can be briefly written in this framework. In addition, we analyze the well-posedness of the linear systems and propose other optimal preconditioners for the Biot model Haga.J;Osnes.H;Langtangen.H2012b based on the mapping property Mardal.K;Winther.R2011a . By proposing optimal block preconditioners, we convert the solution of complicated coupled system into that of a few SPD systems on each of the fields.

The rest of this paper is organized as follows. In Section 2, we give a brief introduction of the Biot model. In Section 3, we introduce two theorems in order to prove well-posedness. In Section 4, we address the unified framework indicating the relationship between preconditioning and well-posedness of linear systems. In Section 5 and Section 6, we show the well-posedness and several optimal preconditioners for the Biot model under the unified framework. In Section 7, we present numerical examples to demonstrate the robustness of these preconditioners.

2 The Biot model

The poroelastic phenomenon is usually characterized by the Biot model Biot.M1941a ; Biot.M1955a , which couples structure displacement , fluid flux , and fluid pressure . Consider a bounded and simply connected Lipschitz domain of poroelastic material. As the deformation is assumed to be small, we assume that the deformed configuration coincides with the undeformed reference configuration. Let denote the total stress in this material. From the balance of the forces, we first have

In addition to the elastic stress

the fluid pressure also contributes to the total stress, which results in the following constitutive equation:

Here, and are the Lamé constants and is the Poisson ratio, the symmetric gradient is defined by , and is the Biot-Willis constant. Therefore, we obtain the following momentum equation

Let denote the fluid content. Then, the mass conservation of the fluid phase implies that


where is the source density. The fluid content is assumed to satisfy the following constitutive equation:


where is the fluid storage coefficient. We also have the Biot-Willis constant in this equation, as this poroelastic model is assumed to be a reversible process and the increment of work must be an exact differential Biot.M1941a ; Rice.J2001a ; CHENG.A2014a . Based on (1) and (2), the following equation holds

According to Darcy’s law, we have another equation:

where is the fluid mobility and is the body force for the fluid phase.

We consider all the parameters to be positive. The following boundary conditions are assumed:


where , and , .

The initial conditions are as follows:

where and are given functions.

We use the backward Euler method to discretize the time derivative :

where is the time step size. More sophisticated implicit time discretizations result in similar linear systems. As we are focusing on the properties of the linear systems resulting from the time discretized problem, we consider only the backward Euler method for the sake of brevity. After the implicit time discretization, fast solvers are needed to solve the following three-field system:


Note that the right-hand side of the last equation in (5) includes terms from previous time step due to the time discretization. As the exact form of this right-hand side does not affect the well-posedness of the linear system, we keep using to denote it. We apply this convention to all the right-hand sides in similar situations, throughout the rest of this paper.

To reduce the number of variables, the fluid flux is eliminated to obtain the following two-field system:


In the rest of this paper, we develop block preconditioners for both the two-field and three-field systems.

3 Well-posedness of linear systems

In this section, we first introduce several theorems to prove the well-posedness of the following saddle point problem: Find such that , the following equations hold


Here, and are given Hilbert spaces with the inner products and , respectively. The corresponding norms are denoted by and .

Given , the following kernel spaces are important in the analysis:

We consider the orthogonal decomposition of and as follows:

We will use these notation to denote the components of functions in the kernel spaces and their orthogonal complements throughout the rest of this section.

The well-posedness of (7) can be proved provided that , , and satisfy certain properties. In the first case, we assume that is coercive on .

Theorem 3.1

(Xu.J2015a ) Assume that is symmetric, is symmetric positive semi-definite. Let be a semi-norm on such that , . Assume that the following inequalities


hold with the constants , , and independent of parameters. In addition, assume that , . Then, Problem (7) is uniformly well-posed with respect to parameters under the norms and , where and .



To prove well-posedness, we just need to verify the following:

  • the boundedness of under the norms and .

  • the Babuska inf-sup condition: for any ,


Since it is straightforward to verify the boundedness of , we focus on proving the inf-sup condition (12).

According to (11), for , consider its projection . There exists such that

Let , , and . Then, we have

Moreover, we have

Therefore, (12) holds.

Remark 1

It is worth noting that in case , we have and .

It is also possible to only assume that is elliptic in . Then we need to assume that is bounded under the norm .

Theorem 3.2

(Brezzi.F;Fortin.M1991a ) Assume that and are symmetric and positive semi-definite and that (8), (9) and (11) hold. Moreover, assume that


Assume that the constants , , , , and are independent of the parameters. Then, Problem (7) is uniformly well-posed with respect to parameters under the norms and .

Theorems 3.1 and 3.2 will be used to prove the well-posedness in different cases. Note that they are sufficient conditions for the problems to be well-posed. For weaker conditions, we refer to Boffi.D;Brezzi.F;Fortin.M2013a .

In this paper, we are especially interested in the robustness of preconditioners with respect to varying material and discretization parameters guided by the well-posedness of the linear system. Thus we want to emphasize the dependence on these parameters in inequalities. Therefore, we introduce the following notation: , and . Given two quantities and , means that there is a constant independent of these parameters such that . can be similarly defined. if and .

4 Relationship between preconditioning and well-posedness

Given that a variational problem is well-posed, an optimal precondtioner can be developed, in order to speed up Krylov subspace methods, such as Conjugate Gradient Method (CG) and Minimal Residual Method (MINRES). In order to illustrate this fact, we first consider the following variational problem:

Find , such that


where is a given Hilbert space and .

The well-posedness of the variational problem (16) refers to the existence, uniqueness, and the stability of the solution. The necessary and sufficient conditions for (16) to be well-posed are shown in the following theorem. We assume the symmetry in the rest of this section.

Theorem 4.1

(Babuska.I1971a ) Problem (16) is well-posed if and only if the following conditions are satisfied:

  • There exists a constant such that .

  • There exists a constant such that


Consider the operator form of (16):

Define operator such that


Assuming the well-posedness, then the following inequalities hold

Therefore, the condition number of the precondtioned system is proved to be bounded

This type of preconditioners is frequently used in the literature and is characterized as “mapping property” in a recent review paper Mardal.K;Winther.R2011a .

Let be a set of given basis of and be a set of given basis of . Consider the matrix representation of and :

and the vector representation of


Assume is symmetric and is SPD. Denote the mass matrix of by , i.e., , . In fact, . Then


Therefore, .

A more general approach is via norm equivalence matrices Loghin.D;Wathen.A2004a . Given an SPD matrix , inner product and norm can be defined correspondingly:

Nonsingular matrices and are H-norm equivalent, denoted by , if there are constants and independent of the size of the matrices such that

If and is symmetric with respect to , then MINRES preconditioned by has the following convergence estimate Loghin.D;Wathen.A2004a :

Consider the preconditioner defined as the matrix representation of in (18). It is easy to see that . Note that .

This can help in the design of preconditioners for CG and MINRES. Preconditioning GMRES differs in that it usually depends on the field of value analysis Loghin.D;Wathen.A2004a .

In the rest of the paper, we will use Theorem 3.1 and Theorem 3.2 to prove the well-posedness of the different formulations of the Biot model under different choices of . Then, based on the well-posedness, we show the corresponding optimal block preconditioners.

5 A two-field formulation

The preconditioning for the two-field system (6) has been studied extensively in the literature Phoon.K;Toh.K;Chan.S;Lee.F2002a ; Toh.K;Phoon.K;Chan.S2004a ; Haga.J;Osnes.H;Langtangen.H2012a ; Haga.J;Osnes.H;Langtangen.H2011a , where the Schur complement approach is usually used to develop preconditioners. In this paper, similar to adler2017robust , we briefly formulate a preconditioner based on the well-posedness of the linear systems for the two-field Biot model.

We first study the well-posedness of (6), beginning by changing the variable in order to symmetrize (6). With an abuse of notation, we still use the notation for pressure after the change of variable. Next, we introduce the function space for the displacement and the pressure. Due to the boundary conditions (3), we consider

for the displacement and

for the pressure. Here, we use the subscript “c” to suggest the continuity of the functions in . We assume in the rest of this paper so that the elasticity operator is nonsingular on . We also assume that such that the divergence operator is surjective on the pressure space.

Then, we define the following bilinear forms:

where and .

Now, we introduce the notation for the kernel spaces:

The variational formulation of (6) is as follows:

Find such that , the following equations hold


We define the norms as follows:


where .

This variational formulation (19) is proved to be well-posed under the norms and provided that the following inf-sup condition holds


It is well known that (21) holds for , and , on a bounded domain with Lipschitz boundary Bramble.J;Lazarov.R;Pasciak.J2001a ; Bramble.J2003a . Moreover, (21) holds for stable Stokes FEM pairs Boffi.D;Brezzi.F;Fortin.M2013a .

Theorem 5.1

Assume that the inf-sup conditions (21) holds and . The system (19) is uniformly well-posed with respect to parameters under the norms and defined in (20).


To prove the well-posedness, we just need to verify the assumptions of Theorem 3.1.

As we assume that , we know that and then (14) is trivial. By definition, (8), (9), and (10) are straightforward to verify.

Based on (21), the following inf-sup condition is implied


Then (11) is verifed.

Therefore, the proof is finished by applying Theorem 3.1.

With the well-posedness of (19) proved, an optimal preconditioner can be formulated. We first introduce some matrix notation. Given finite element basis functions and for and , respectively, define the following stiffness matrices: , , and .

The matrix form of the system and preconditioner are


Remark 2

In case , the kernel space contains constant functions. We can similarly prove the well-posedness, but the norm has a term , which results in a dense matrix in the preconditioner. We refer to Lee.J;Mardal.K;Winther.R2015a for constructing preconditioners related to .

In the literature, the preconditioners for two-field formulation are mostly based on Schur complement approaches. The exact Schur complement preconditioner of , i.e.,

is known to be an optimal preconditioner Elman.H;Silvester.D;Wathen.A2005a , although is dense and cannot be obtained. Practical approximations of , such as

have also been investigated Phoon.K;Toh.K;Chan.S;Lee.F2002a ; Toh.K;Phoon.K;Chan.S2004a ; Haga.J;Osnes.H;Langtangen.H2012a ; Haga.J;Osnes.H;Langtangen.H2011a .

The two-field formulation is usually considered computationally efficient, as it involves the fewest variables and, therefore, has smaller linear systems to solve than the three-field formulation (5). However, the two-field formulation (with continuous pressure elements) exhibits oscillations in the pressure field, and more expanded systems such as the three-field formulation, are shown to be more stable Ferronato.M;Castelletto.N;Gambolati.G2010a ; Haga.J;Osnes.H;Langtangen.H2012b . Motivated by this fact, we study a three-field formulation Haga.J;Osnes.H;Langtangen.H2012b in the next section.

6 A three-field formulation

In this section, we will show the well-posedness of the three field formulation (5), briefly formulate the diagonal block robust preconditioners of hong2018parameter ; adler2019robust as special cases, and propose some new preconditioners for the three field formulation guided by the well-posedness.

6.1 The three-field formulation

We can write (5) as a symmetric problem by rescaling. Introduce

The three-field system (5) can be rewritten as


With an abuse of notation, we still use and to denote the scaled velocity , the scaled pressure , respectively. Then, we introduce the function spaces:

and bilinear forms

We define the corresponding kernel spaces related to

Note that due to the assumption , we have .

Then, the weak formulation is as follows:

Find and such that and , the following equations hold


The additional term corresponds to different versions of the Biot models Axelsson.O;Blaheta.R;Byczanski.P2012a .

The well-posedness of this saddle point problem can be proved with different choices of norms for and . We discuss some of these options in the rest of this section.

6.2 Augmented Lagrangian preconditioners

The stability of the three-field system (23) is closely related to the stability of the pair - and -. In particular, it is considered stable if - satisfies (21) and - satisfies



(25) holds for and and, in discrete cases, there are many stable pairs, such as Raviart-Thomas elements Raviart.P;Thomas.J1977a for and piecewise polynomials for .

The augmented Lagrangian (AL) method Benzi.M;Olshanskii.M2006a ; Xu.J;Yang.K2015a incorporates the constraint into the norm. The constraint here is

Therefore, it is natural to consider the following norms for the AL method.

Let be the projection from to . We define the norms for spaces and as follows:


where is the coefficient in bilinear form , and is an undetermined parameter.

To prove the well-posedness of (24), we just need to verify the assumptions of Theorem 3.1.

Given (21), we have


For the case in which , the right-hand side of (27) is equal to .

Given (25), we can prove another inequality:


Similarly, if we further assume that , the right-hand side of (28) is equal to . Note that this approach is used in Lipnikov.K2002a , where the displacement is set to be zero and the inf-sup condition of the - pair is assumed.

The boundedness of is easy to verify due to the additional term in the norm :


The coercivity of is straightforward to prove, as


Because is uniformly coercive only on , we resort to Theorem 3.2 to prove the well-posedness.

Theorem 6.1

Assume , is uniformly bounded and the inf-sup conditions (21) and (25) hold. Then the system (23) is uniformly well-posed with respect to parameters under the norms and defined in (26).


As , (14) is trivial to prove. Consider for the inf-sup condition of . Due to , the right-hand side of (27) or (28) is equal to . Therefore, the inf-sup condition of is proved.

As , we can prove that . Therefore, the assumptions of Theorem 3.2 hold. Then the proof is finished by applying Theorem 3.2.

It is obvious that we only need to assume either (21) or (25) to prove the well-posedness of (24).

Corollary 1

Assume , is uniformly bounded, and that the inf-sup condition (21) holds. The system (23) is uniformly well-posed with respect to parameters under the norms defined in (26).


The proof follows from (27), (30), (9), and Theorem 3.2.

Corollary 2

Assume that , is uniformly bounded, and the inf-sup condition (25) holds. The system (23) is uniformly well-posed with respect to parameters under the norms defined in (26).


The proof follows from (28), (30), (9), and Theorem 3.2.

Remark 3

The assumption that both (21) and (25) hold results in a smaller parameter than the cases where only one of (21) and (25) holds.

Based on the well-posed formulation, we derive the corresponding optimal block diagonal preconditioner.

6.2.1 Matrix form

We introduce some additional matrix notation. Also, we introduce the FEM basis for . Define the stiffness matrices , , , and .

Then the system matrix of the three-field formulation is

The block preconditioner is