Parameter-Robust Preconditioners For Biot's Model

SeMA Journal (2024) 81:51–80
https://doi.org/10.1007/s40324-023-00336-2
Parameter-robust preconditioners for Biot’s model
Carmen Rodrigo1 · Francisco J. Gaspar1 · James Adler2 · Xiaozhe Hu2 ·

Peter Ohm3 · Ludmil Zikatanov4
Received: 15 April 2023 / Accepted: 19 July 2023 / Published online: 5 August 2023
© The Author(s) 2023
Abstract
This work presents an overview of the most relevant results obtained by the authors regarding
the numerical solution of the Biot’s consolidation problem by preconditioning techniques.
The emphasis here is on the design of parameter-robust preconditioners for the efficient
solution of the algebraic system of equations resulting after proper discretization of such
poroelastic problems. The classical two- and three-field formulations of the problem are con-
sidered, and block preconditioners are presented for some of the discretization schemes that
have been proposed by the authors for these formulations. These discretizations have been
proved to be well-posed with respect to the physical and discretization parameters, what pro-
vides a framework to develop preconditioners that are robust with respect to such parameters
as well. In particular, we construct both norm-equivalent (block diagonal) and field-of-value-
equivalent (block triangular) preconditioners, which are proved to be parameter-robust. The
theoretical results on this parameter-robustness are demonstrated by considering typical
benchmark problems in the literature for Biot’s model.
B Carmen Rodrigo
[email protected]
Francisco J. Gaspar
[email protected]
James Adler
[email protected]
Xiaozhe Hu
[email protected]
Peter Ohm
[email protected]
Ludmil Zikatanov
[email protected]
1 Department of Applied Mathematics, University of Zaragoza, Pedro Cerbuna, 12, 50009 Zaragoza,
Zaragoza, Spain
2 Department of Mathematics, Tufts University, 177 College Avenue, Medford, MA 02155, USA
3 Department of Mathematics, RIKEN Center for Computational Science, 7-1-26
Minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan
4 Department of Mathematics, The Pennsylvania State University, 107 McAllister Building, University
Park, State College, PA 16802, USA
123
52 C. Rodrigo et al.
Keywords Biot’s model · Poroelasticity · Well-posedness · Preconditioning ·

Parameter-robust preconditioners
Mathematics Subject Classification 65F08 · 65F10 · 65N12 · 65N22 · 65N30
1 Introduction
Coupling of fluid flow and mechanical deformation within a porous media is an important
multi-physics problem appearing in many applications in several branches of technology and
natural sciences. Such coupling was already modeled in the early one-dimensional work of
Terzaghi [1], whereas the general three-dimensional mathematical formulation was estab-
lished by Maurice Biot in several pioneering publications (see [2, 3]). Nowadays, geothermal
energy extraction, CO2 storage, hydraulic fracturing or cancer research are among typical
societal relevant applications of this model. In particular, here we consider the quasi–static
Biot’s model for soil consolidation, assuming the porous medium to be linearly elastic,
homogeneous, isotropic, and saturated by a Newtonian fluid. According to Biot’s theory [2],
the consolidation process satisfies a system of partial differential equations (PDEs) which
combines the equations for fluid flow and the elastic deformation of the solid matrix. If we
consider a bounded open subset ⊂ Rd , d = 2, 3 with regular boundary , the governing
equations are given as follows,
−div σ + α∇ p = ρ g, in , (1)

∂ 1
p + α div u + div w = f , in . (2)
∂t Mb
Equation (1) corresponds to the equilibrium equation for the porous medium, which is given
in terms of the effective stress σ , and where p is the pore pressure which is multiplied by the
Biot–Willis constant α. The effective stress is related to the strain tensor by the constitutive
equation σ = 2με(u) + λ div(u)I being λ and μ the so-called Lamé coefficients, which
can be computed in terms of the Young modulus, E, and the Poisson ratio, ν, as follows:
Eν E
λ= and μ = .
(1 − 2ν)(1 + ν) 1 + 2ν
Under the assumptions considered here, the strain tensor can be written in terms of the
displacement vector u as follows ε(u) = 21 (∇u + ∇ut ). Finally, the right-hand side in (1)
is given as the vector-valued function g, which represents the gravitational force, multiplied
by the bulk density ρ = φρ f + (1 − φ)ρs , where ρs and ρ f are the densities of solid and
fluid phases and φ is the porosity.
The fluid mass conservation equation (2) is written in terms of the variation of fluid content,
being Mb the bulk modulus whose inverse is called storage coefficient and is related to the
compressibility of the fluid, and the percolation velocity of the fluid relative to the soil w. The
source term f represents a forced fluid extraction or injection process. Finally, Darcy’s law
1
relates the fluid’s flux w to the fluid’s pressure gradient as follows w = − K (∇ p − ρ f g),
μf
where K is the absolute permeability tensor of the porous medium and μ f is the viscosity
of the fluid.
All these equations can be combined to obtain different formulations of the poroelasticity
problem, depending on which variables one is interested in. The most widely used are the
two field formulation in which the main variables are the displacement vector and the fluid
123
Parameter-robust preconditioners for Biot’s model 53
pressure, and the classical three-field formulation which includes also Darcy velocity in order
to guarantee fluid mass conservation. But, of course, there are other different formulations, for
example the three-field formulation in which the solid pressure is included as third variable as
a way to deal with low-compressible or incompressible materials, or the recently introduced
three-field formulation with the total pressure [4] which is a weighted sum of fluid and solid
pressures. In this work, we will focus on the classical two- and three-field formulations, but
the techniques used here can be applied to any of the formulations of the problem.
As we can observe in the system of equations presented above, many physical parameters
are involved in the model, whose values vary for different real application problems. In fact,
it is important to note that the values of some of these parameters may differ over several
orders of magnitude depending on the application under consideration. As an example, the
permeability can typically range from 10−9 to 10−21 m2 in geophysical applications [4, 5],
whereas in biophysical applications such as in the modeling of soft tissue or bone, the values
for this parameter can range from 10−14 to 10−16 m2 [6–8]. Due to this fact, the design
of discretization techniques and solution algorithms that behave well with respect to this
variation of the parameters is a big challenge and it has received a lot of attention during the
last years.
Although finite difference and finite volume methods have been considered for the dis-
cretization of Biot’s model in the literature [9–12], here we focus on the application of finite
element methods. Stable finite-element schemes have been developed for the different for-
mulations of Biot’s model. For the classical two-field formulation, Taylor–Hood elements
which satisfy an appropriate inf-sup condition have been used [13–15], and also appropriate
stabilization techniques have been proposed for unstable finite-element pairs, such as the com-
bination of piecewise linear elements for both variables [16–18]. For the classical three-field
formulation, a non-conforming stable discretization based on Crouzeix–Raviart elements
for displacement, lowest order Raviart–Thomas–Nédélec elements for Darcy’s velocity, and
piecewise constants for the pressure was proposed in [19]. In [20] a stabilization by bub-
ble functions was presented to obtain a parameter-robust version of a typically considered
discretization based on piecewise linear elements for displacements, Raviart–Thomas for
Darcy’s velocity and piecewise constants for the pressure. Another parameter-independent
approach can be found in [21], where mixed finite element discretizations that fulfills mass
conservation strongly at a discrete level are considered. In addition, stable discretization
schemes have been proposed for other formulations, as for example, for the three-field for-
mulation including the total pressure [4, 22], or the four-field formulation considering stress
tensor, fluid flux, displacement, and pore-pressure as main variables [23]. More recently, in
[24] the authors deal with conservative discretizations of the four-field formulation of Biot’s
equations utilizing total pressure.
Typical discretizations of Biot’s model yield a large-scale sparse linear system of equations
that has to be solved at each time step. Such linear systems are usually ill-conditioned and
difficult to solve in practice. In fact, the solution of these systems of equations constitutes the
main bottleneck in the numerical simulation of coupled poromechanics problems, which is
mandatory for real applications, and therefore the design of efficient and fast solvers for Biot’s
model has become increasingly popular. Due to the size and the properties of the resulting
systems, the application of iterative solution techniques, instead of direct methods, is usual.
In particular, for the poroelastic problems considered here, there are two typical solution
approaches: iterative coupling methods and fully-coupled or monolithic methods. Whereas
monolithic techniques solve the resulting linear system simultaneously for all the involved
unknowns, iterative coupling methods [25, 26], in contrast, consist of a sequential approach
in which either the fluid flow problem or the geomechanics part is solved first, followed by
123
the solution of the other system, repeating this process until a converged solution within a
prescribed tolerance is obtained. These iterative coupling methods have the main advantage
of being able to combine existing software for simulating fluid flow and geomechanics
in order to obtain the solution for the coupled problem. For this reason, the analysis and
application of this type of schemes have received a lot of attention [27–31]. In particular,
in [32, 33] a re-interpretation of the four commonly used sequential splitting methods as
preconditioned-Richardson iterations with block-triangular preconditioning is presented, and
such analysis indicates that a fully-implicit method outperforms the convergence rate of the
sequential-implicit methods. In the context of monolithic methods, one can develop efficient
preconditioners to accelerate the convergence of Krylov subspace methods or to design special
smoothers within a multigrid framework. Examples of this latter approach can be found in
[34–38], whereas the design of preconditioners for poromechanics has been addressed, for
example, in [4, 9, 39–44], where many of the developed techniques are based on the Schur
complement approach. More recently, the focus of this latter approach has been on the design
of parameter-robust preconditioners [4, 45–48].
Thus, in this work we want to present some of the results obtained by the authors within
this context. In particular, we address the construction of preconditioners for Biot’s model,
which are robust with respect to typical scales in material constants spanning over many
orders of magnitude. This is done for different finite element discretizations of the two clas-
sical formulations of Biot’s model, namely the displacement–pressure two-field formulation
and the three-field formulation including the Darcy’s velocity. For the first one, stabilized
finite element pairs proposed and analyzed by the authors in [16, 17] are considered. For the
three-field formulation, we propose preconditioners for two different discretizations: the non-
conforming scheme given in [19] which is based on Crouzeix–Raviart finite elements for the
displacements, lowest order Raviart–Thomas–Nédélec elements for the Darcy velocity, and
piecewise constant approximation for the pressure, and the conforming stabilized scheme pro-
posed in [20] which assumes piecewise linear elements for displacements, Raviart–Thomas
finite elements for Darcy’s velocity and piecewise constants for the pressure. With this pur-
pose, we develop robust block preconditioners (e.g. [49, 50]) to accelerate the convergence
of Krylov subspace methods. These preconditioners take advantage of the block structure of
the discrete problem, decoupling the different fields at the preconditioning stage. Such block
preconditioning is primarily attractive due to its simplicity, which allows us to focus on the
character of the diagonal blocks, and to leverage extensive work on solving simpler problems.
For example, algebraic multigrid [51] or auxiliary space decomposition methods [52, 53] can
be considered for efficiently solving such blocks. Since the considered discretizations upon
we build the preconditioners are well-posed with respect to the physical and discretization
parameters, we are able to develop robust block preconditioners that efficiently solve the
linear systems, independently of such parameters as well. Following the general framework
developed in [33, 43, 46, 54, 55], we first consider block diagonal preconditioners (also
known as norm-equivalent preconditioners). Then, we discuss block triangular (upper and
lower) preconditioners following the framework developed in [54, 56, 57] for Field-of-Value
(FOV) equivalent preconditioners. For both cases, we show that the theoretical bounds on
their performance remain independent of the discretization and physical parameters of the
problem.
We present theoretical results that guarantee the parameter robustness of the proposed
preconditioners, but also we present appropriate numerical experiments to illustrate their
robust performance. In particular, we will consider typical benchmark problems that are
used in the literature to analyze the behavior of numerical methods for poromechanics. All
these test problems were implemented in the HAZmath library [58], developed by some
123
of the authors, which contains routines for finite elements, multilevel solvers, and graph
algorithms. The numerical tests were performed on a workstation with an 8-core 3 GHz Intel
Xeon “Sandy Bridge” CPU and 32 GB of RAM per core.
The rest of the paper is organized as follows. Section 2 presents the general theory in
which is based the design of the proposed block preconditioners. Next, Sect. 3 is divided
into Sects. 3.1 and 3.2 which focus on the two- and three-field formulations of Biot’s model,
respectively. Section 3.1 is composed of four parts: first the formulation of the problem is
presented, followed for two subsections with the proposed block diagonal and block triangular
preconditioners respectively, and last a numerical experiment showing the performance of the
proposed solvers for a benchmark problem in the literature. Section 3.2 first introduces the
considered three-field formulation, and then it is divided into Sects. 3.2.1 and 3.2.2, dealing
with two different discretizations for such a formulation, and whose structure is the same as
that of Sect. 3.1. Finally, concluding remarks are made in Sect. 4.
2 Block preconditioning
Our aim is to derive uniform block preconditioners for the Biot’s model based on the well-
posedness of the considered discretizations. In order to do this, we will follow the general
theory presented in [54, 55], and both block diagonal and triangular preconditioners will be
developed. In the following, we introduce the general framework in which we will construct
such preconditioners.
Let X be a real, separable Hilbert space equipped the inner product (·, ·) X , which induces
the norm · X . Let us denote X as a dual space to X, being ·, · the duality pairing between
them. Then, we can consider an operator A : X → X induced by a symmetric and bounded
bilinear form L(·, ·), that is, A x, y = L(x, y). We assume that such a bilinear form satisfies
the following continuity and inf-sup conditions:
|L(x, y)| ≤ βx X y X , ∀x, y ∈ X, (3)
L(x, y)
inf sup ≥ γ, (4)
x∈X y∈X x X y X
where β, γ > 0. The properties of the bilinear form ensure that A is a bounded and symmet-
ric operator and that the arising system of equations A x = b is well-posed. A competitive
way for solving such a system of equations is to employ an iterative method coupled with
a preconditioning technique. More precisely, we look for a preconditioner M such that
MA x = M b or AM x = b (x = M−1 x) can be solved in a number of iterations indepen-
dent of the system’s size. We can manage this by choosing a preconditioning matrix which
is equivalent to the system matrix in some sense. Thus, in the next subsections we will con-
sider norm and field-of-values equivalence of matrices, which will provide us convergence
bounds for Krylov subspace methods such as the Minimum Residual method (MINRES) or
the Generalized Minimum Residual method (GMRES). In addition, if such equivalences hold
independently of the parameters involved in the problem, we will obtain a robust convergence
with respect to such parameters. With these assumptions, we will be able to design block
diagonal and block triangular parameter-robust preconditioners for solving system A x = b .
123
2.1 Norm-equivalent preconditioner
Let M : X → X be a symmetric positive definite (SPD) operator, which induces an

inner product (x, y)M−1 := M−1 x, y on X and the corresponding norm x2M−1 :=
(x, x)M−1 . Since MA : X → X is symmetric with respect to (·, ·)M−1 , we can use M as
preconditioner for the MINRES algorithm, and it is well-known [59] that the convergence
rate of the preconditioned MINRES satisfies that
κ(MA) − 1
A(x − x m )M ≤ 2ρ m A(x − x 0 )M , with ρ = , (5)
κ(MA) + 1
being x m the m-th iteration of the preconditioned MINRES, x the exact solution, and κ(MA)
denotes the condition number of MA. Then, we are interested in estimating the condition
number κ(MA) in order to know about the behavior of preconditioner M.
Mardal and Winther [55] showed that, if the well-posedness conditions (3)–(4) hold, and
M satisfies
c1 x2X ≤ x2M−1 ≤ c2 x2X , (6)
then, A and M are norm-equivalent and we have the following bound for the condition
number κ(MA) ≤ cc12 γβ .
This implies that the convergence rate of the MINRES method preconditioned with M
satisfies ρ ≤ cc22 β−c 1γ
β+c1 γ , and therefore, if the original problem is well-posed and the norm-
equivalent inequality (6) is satisfied with constants c1 and c2 independent of the physical and
discretization parameters, then the convergence rate of preconditioned MINRES is uniform,
and M is a parameter-robust preconditioner. Later we will use this approach to design norm-
equivalent preconditioners for the Biot’s model and, as we will see, this will lead to block
diagonal preconditioners.
In practice, a natural choice of norm-equivalent preconditioner is the Riesz operator B :
X → X corresponding to the inner product (·, ·) X , that is,
(B f , x) X = f , x, ∀ f ∈ X , x ∈ X. (7)
It can be easily seen that B satisfies the norm-equivalent inequality (6) with constants c1 and
c2 equal to one, and then κ(BA) ≤ γβ . Thus, if the well-posedness constants β and γ are inde-
pendent on the physical and discretization parameters, B will provide a robust preconditioner.
2.2 FOV-equivalent preconditioner
In this section, we study a more general class of preconditioners as the field-of-values-

equivalent (FOV-equivalent) preconditioners are. Let us consider a general operator M L :
X → X. M L and A are said to be FOV-equivalent if there exist constants and ϒ such
that, for any x ∈ X,
(M L A x, x)M−1 M L A xM−1
≤ , ≤ ϒ. (8)
(x, x)M−1 xM−1
In general, M : X → X can be any SPD operator, but we will consider M to be a SPD

norm-equivalent preconditioner, and and ϒ are positive constants such that ≤ ϒ.
Such an operator M L can be used as a preconditioner for the Generalized Minimum Residual
123
(GMRES) algorithm, in the way that, from [60, 61], the convergence rate of the preconditioned
GMRES method is given by
2 m
M L A(x − x m )2M−1 ≤ 1 − 2 M L A(x − x 0 )2M−1 , (9)
ϒ
where x m is the m-th iteration of the GMRES method preconditioned with M L and x is the
exact solution.
If constants and ϒ are independent of the physical and discretization parameters, then
M L is a uniform left preconditioner for GMRES and is referred to as an FOV-equivalent
preconditioner. This approach leads to block lower triangular preconditioners.
Similar arguments also apply to right preconditioners for GMRES, so that MU : X → X
and A, are FOV equivalent if, for any x ∈ X ,
(AMU x , x )M AMU x M
≤ , ≤ ϒ. (10)
(x , x )M x M
In this case, if and ϒ are independent of the physical and discretization parameters, MU is
a uniform right preconditioner for GMRES. Such an approach leads to block upper triangular
preconditioners.
Naturally, FOV-equivalent preconditioners can be built based on the Riesz operator and
the FOV-equivalence which can be derived from the well-posedness conditions (3)–(4), as
done for example in [54].
In the next section, we apply the theoretical framework presented here to design appropri-
ate block preconditioners for the chosen discretizations of Biot’s model. We have to take into
account that proper weighted norms have to be chosen so that the well-posedness constants are
robust with respect to the physical and discretization parameters and thus the preconditioners
also maintain this robustness property.
3 Parameter-robust preconditioners for Biot’s model
3.1 Two-field formulation
First, we deal with the classical two-field formulation of Biot’s model (1)–(2), where the
main unknowns are the displacement vector u and the pore pressure p.
Let L 2 () be the Hilbert space of square integrable scalar-valued functions defined on
, and let H 1 () denote the space of square integrable scalar-valued functions whose first
derivatives are in L 2 (), and H 1 () its vector counterpart. We consider appropriate Sobolev
spaces V ⊂ H 1 () and Q ⊂ H 1 (), which would include the information about the
imposed boundary conditions, and by applying integration by parts, we obtain the following
variational form:
For each t ∈ (0, T ], find (u(t), p(t)) ∈ V × Q, such that
a(u, v) − α(div v, p) = (ρ g, v), ∀v ∈ V , (11)

1 ∂p ∂u
− q − α div , q − a p ( p, q) = ( f , q), ∀q ∈ Q, (12)
Mb ∂t ∂t
where bilinear forms a(·, ·) and a p (·, ·) are given as follows:

a(u, v) = 2μ ε(u) : ε(v) + λ div u div v,

123

K
a p ( p, q) = ∇ p · ∇q,
μf
(·, ·) denotes the standard inner product on L 2 (), and the problem is completed with suitable
initial conditions.
In the rest of the section, we neglect the first term in the flow equation, assuming that the
fluid is incompressible, since this results in the most challenging case.
Discretization
In order to introduce the finite element spatial discretization, we consider finite dimensional
spaces V h ⊂ V and Q h ⊂ Q. In particular, here we consider two stabilized discretizations
for the two-field formulation of Biot’s model proposed in [16, 17]. For the first one, V h = V l
is the space of piece-wise (with respect to a triangulation Th ) linear continuous vector-valued
functions on and Q h consists of piece-wise linear continuous scalar valued functions. This
finite element pair is referred as the P1 -P1 scheme. The second considered discretization is
the so-called MINI element, introduced in [62], where Q h is as before and V h = V l ⊕ V b ,
being V b the so-called space of bubble functions defined as
V b = span{ϕb,T e1 , . . . , ϕb,T ed }T ∈Th , ϕb,T = αT λ1,T . . . λd+1,T ,
where λm,T are the barycentric coordinates on T , e j are the canonical Euclidean basis vectors
in Rd and αT is a normalizing constant for ϕb,T . Both finite element pairs have to be stabilized.
The P1 –P1 scheme needs a stabilization term to satisfy a uniform inf-sup condition, whereas
both discretizations need such an stabilization for eliminating non-physical oscillations of the
pressure approximation seen in practice [16]. Taking this into account, and using an implicit
Euler scheme for time discretization on a time interval (0, tmax ] with constant time-step size
τ , the fully discretized scheme at time tn , n = 1, 2, . . . is as follows:
Find unh ∈ Vh and phn ∈ Q h , such that,
a(unh , v h ) − α(div v h , phn ) = (g(tn ), v h ), ∀v h ∈ V h , (13)

− α(div ∂¯t unh , qh ) − a p ( phn , qh ) − ηh 2 (∇ ∂¯t phn , ∇qh ) = ( f (tn ), qh ), ∀qh ∈ Q h , (14)
where ∂¯t unh := (unh − un−1 ¯ n n−1

h )/τ , ∂t ph := ( ph − ph )/τ , and η represents the stabilization
n
parameter.
At each time step, the resulting discrete problem can be written as the following two-by-
two block linear system:

Au α BT u fu
A x = b, A = , x = , and b = , (15)
α B −τ A p − ηh 2 L p p fp
where a(u, v) → A u , −(div u, q) → B, a p ( p, q) → A p , and (∇ p, ∇q) → L p represent

the discrete versions of the variational forms.
Our target is to design block diagonal and triangular preconditioners for the solution of
this discrete problem. In order to do that, as previously mentioned, we follow the frame-
work proposed in [54, 55], which is based on the well-posedness of the discretized linear
system at each time step. This well-posedness was already proved in [17]. In order to obtain
parameter-robust preconditioners, however, we need to prove this well-posedness of the lin-
ear system (15), making special emphasis on the independence of the constant arising from
the analysis with respect to any physical and discretization parameters.
123
With this purpose, the key is to define a proper weighted norm on X = V h × Q h as

follows. For x = (u, p)T , let us define,
α2
x2X := u2A u + τ p2A p + ηh 2 p2L p + p2 , (16)
ζ2

2μ
where u2A u := a(u, u), p2A p := a p ( p, p), p2L p := (∇ p, ∇ p), ζ = λ+ d , · is
the standard L 2 −norm,
and d = 2 or 3 is the dimension of the problem.
Next, we define the following composite bilinear form,
L(x, y) = (A u u, v) + α(Bv, p) + α(Bu, q) − τ (μ−1

f K∇ p, ∇q) − ηh (∇ p, ∇q),
2
(17)
for x = (u, p) and y = (v, q), and the well-posedness of problem (13)–(14) is demonstrated
by the following theorem.
Theorem 1 The bilinear form L(x, y) given in (17) satisfies the continuity and inf-sup
conditions (3)–(4) with constants γ and β independent of the physical and discretization
parameters, and A defined in (15) is therefore an isomorphism from X to X provided that
the stabilization parameter, η, satisfies η = δ αζ 2 with δ > 0.
2
Proof Due to our choice of finite element spaces X = V h × Q h , we have the following
inf-sup condition (see [63]),
(div v, q)
sup ≥ γ B0 q − ξ 0 h∇q, ∀ q ∈ Q h , (18)
v∈V h v1
where γ B0 > 0 and ξ 0 ≥ 0 are constants that do not depend on the mesh size (notice that if we
√
use the MINI element, then ξ 0 = 0). From the definition of ζ , it holds that v A u ≤ dζ v1 ,
and therefore we can reformulate the inf-sup condition, (18), as follows,
(Bv, q) (Bv, q) γ0 ξ0 γB ξ
sup ≥ sup √ ≥ √ B q − √ h∇q =: q − h∇q,
v∈V h v A u v∈V h dζ v1 dζ dζ ζ ζ
(19)
√ √
where γ B := γ B0 / d and ξ = ξ 0 / d.
(18) and (19), for any p, there exists w ∈ V h such that
Based on the inf-sup conditions
(Bw, p) ≥ γζB p − ζξ h∇ p w A u and, without loss of generality, we may assume
w A u = p, and then we have,

γB ξ
(Bw, p) ≥ p − h∇ p p. (20)
ζ ζ
For a given pair (u, p) ∈ V h × Q h and w defined as above, we choose v = u+θ w, θ = ϑ γ Bζ α
and q = − p and using the facts above, we have,
L(x, y) = (A u u, u+θ w)+α(B(u + θ w), p)−α(Bu, p)+τ (K ∇ p, ∇ p)+ηh 2 (∇ p, ∇ p)

γB α γ 2 α2 γ B α2
≥ u2A u − ϑu A u p + ϑ B 2 p2 − ϑ 2 ξ h∇ p p
ζ ζ ζ
δ α2 2 2
+ τ p2A p + ξ h ∇ p2
ξ2 ζ2
123
⎛ ⎞T ⎛ ⎞⎛ ⎞
u A u 1 −ϑ/2 0 0 u A u
⎜ γ B α p ⎟ ⎜ ⎟⎜ γ α ⎟
⎜ ζ ⎟ ⎜−ϑ/2 ϑ −ϑ/2 0⎟ ⎜ Bζ p ⎟
≥⎜
⎜ α ξ h∇ p⎟
⎟ ⎜
⎜ 0
⎟⎜ ⎟.
⎝ζ ⎠ ⎝ −ϑ/2 δ/ξ 2 0⎟ ⎜α
⎠ ⎝ ζ ξ h∇ p⎠
⎟
√ √
τ p A p 0 0 0 1 τ p A p
If 0 < ϑ < min{2, ξ2δ2 }, the matrix in the middle is SPD and there exists γ0 such that

γ 2 α2 α2
L(x, y) ≥ γ0 u2A u + B 2 p2 + 2 ξ 2 h 2 ∇ p2 + τ p2A p ≥ γ̃ x2X ,
ζ ζ
where γ̃ = γ0 min{γ B2 , ξ 2 /δ}. Since, it is straightforward to verify also that (v, q)2X ≤
γ̄ 2 (u, p)2X , then L satisfies the inf-sup condition (4) with γ = γ̃ /γ̄ . The boundedness of
L, that is (3), is obtained by the continuity of each term and the Cauchy–Schwarz inequality.

Norm-equivalent preconditioners
Once we have proved that the well-posedness condition (4) holds, we look for SPD operators
that satisfy the norm-equivalence inequality (6). As mentioned in Sect. 2.1, one natural choice
is the Riesz operator, which for the two-field stabilized discretization and the norm · X
defined in (16) is given as follows:
−1
Au 0
BD = α2 , (21)
0 τ A p + ηh 2 L p + ζ2
M
where M is the mass matrix of the pressure block.

In practice, applying preconditioner B D implies inverting the diagonal blocks exactly,
which can be expensive and sometimes infeasible. Thus, we consider the following inexact
preconditioner:

D = Su 0 ,
B
0 Sp
where the diagonal blocks Su and S p are SPD and spectrally equivalent to the diagonal blocks
in B D , that is,
c1,u (Su u, u) ≤ (A−1

u u, u) ≤ c2,u (Su u, u), (22)
−1
α2
c1, p (S p p, p) ≤ τ A p + ηh L p + 2 M
2
p, p ≤ c2, p (S p p, p), (23)
ζ
being the involved constants independent of discretization and physical parameters.

D and A are norm-equivalent and we can estimate the condition number κ(B
B D A) = O(1)
by Theorem 1.
FOV-equivalent preconditioners
For a more general approach, we propose now block triangular preconditioners for the sta-
bilized system given in (15), following the FOV-equivalent framework presented in Sect. 2.2.
For simplicity in the analysis, we slightly modify the system matrix A by negating the second
equation, but still denote it by A.
123
First, we consider the following block lower triangular preconditioner:

−1
Au 0
BL = , (24)
−α B τ A p + ηh 2 L p + αζ 2 M
2
and we need to show that the FOV-equivalence in (8) is satisfied.
Theorem 2 There exist constants and ϒ, independent of discretization or physical param-

eters, such that, for any x = (u, p)T = 0,
(B L A x, x)(B D )−1 B L A x(B D )−1
≤ , ≤ ϒ,
(x, x)(B D )−1 x(B D )−1
provided that η = δ αζ 2 with δ > 0.

2
Proof By direct computation first, and applying the reformulated inf-sup condition given
in (19) secondly,
(B L A x, x)(B D )−1 = (u, u) A u + α(B T p, u) + τ ( p, p) A p
+ ηh 2 (L p p, p) + α 2 (B A−1 T
u B p, p)

≥ 0 u2A u + τ p2A p + ηh 2 p2L p + α 2 B T p2A−1

u
2
γ B ξ
≥ 0 u2A u + τ p2A p + ηh 2 p2L p + α 2 p − h∇ p .
ζ ζ
α 2
2 < θ < 1, and taking into account that η = δ ζ 2 with δ > 0, we can
1
By choosing 1+δ/ξ
obtain the lower bound in the theorem:

(B L A x, x)(B D )−1 ≥ 0 u2A u + τ p2A p
2
γ B2 α 2 δ 1 α
+(1 − θ ) 2 p2 + 1 + 2 − ξ 2 h 2 ∇ p2
ζ ξ θ ζ2

α2
≥ 0 1 u2A u + τ p2A p + ηh 2 p2L p + 2 p2
ζ
=:
(x, x)(B D )−1 ,
2
where 1 := min{1, (1 − θ )γ B2 , 1 + ξδ2 − θ1 ξδ }. The upper bound ϒ can be obtained
directly from the continuity of each term, the Cauchy-Schwarz inequality, and the fact that
(y,z)
B T p A−1
u
≤ ζ1 p. To show the last inequality we use the identity z A u = y AAu . We
u
then have

−1 T
B T p A−1
u
= A u B p, B Tp = A−1 T −1 T
u B p, A u B p
Au
−1

v, A u B p A
T
= A−1u B p A u = sup
T u
v v A u
(Bv, p) Bv L 2
= sup ≤ p sup . (25)
v v A u v v A u
123
Since by definition we have that v2A u = 2με(v)2L 2 + λBv2L 2 and also the spatial
dimension satisfies d ≥ 1, it is evident that
2μ 2μ
Bv2L 2 = ε11 (v) + · · · + εdd (v)2L 2
d d (26)

≤ 2μ ε11 (v)2L 2 + · · · + εdd (v)2L 2 ≤ 2μ ε(v)2L 2 .
Therefore, we obtain

2μ
ζ 2 Bv2L 2 = + λ Bv2L 2
d (27)
≤ 2μ ε(v)2L 2 + λBv2L 2 = v2A u .
Bv L 2
Taking square root and a supremum over v shows that supv v A u ≤ ζ1 .

The inexact counterpart of (24) is also considered,

−1 −1
L = Su
B
0
−1 , (28)
−α B S p
where Su and S p satisfy the spectrally equivalences (22) and (23). Then, similarly, we can
L is also FOV-equivalent with A and, therefore, it can be used as preconditioner
show that B
for GMRES. Due to the fact that the proofs are similar, we only state the results here. The
next theorem is easily seen to be analogous to [64, Theorem 4.4].
Theorem 3 If I −Su A u A u ≤ ρ with 0 ≤ ρ < 1, there exist constants and ϒ, independent

of discretization and physical parameters, such that, for any x = (u, p)T = 0, it holds that
L A x, x)(B )−1
(B L A x(B )−1
B
≤ D
, D
≤ ϒ,
(x, x)(BD )−1 x(BD )−1
provided that η = δ αζ 2 with δ > 0.

2
In a similar way, we can construct block upper triangular preconditioners, both exact and
inexact versions,
−1 −1 −1
Au α BT Su α BT
BU =
and BU = , (29)
0 τ A p + ηh 2 L p + αζ 2 M S −1
2
0 p
and we can prove the FOV-equivalence (10), so that they give rise to uniform right precon-
ditioners for GMRES.
Numerical experiment: 3D footing problem
In order to demonstrate the robustness of the proposed preconditioners, we consider a clas-
sical benchmark in the literature, that is, a 3D footing problem (see for example [34]).
The computational domain, depicted in the left side of Fig. 1, is a three dimensional cube
= (−32, 32) × (−32, 32) × (0, 64) that represents a block of porous soil, on which a
uniform load of intensity 0.1 N/m2 is applied in a square of size 32 × 32 m2 at the middle
of the top of the domain. The base of the domain is fixed in all directions, and the remaining
boundary is considered to be stress free. A no-flow boundary condition is applied at the
compression zone, and the remaining boundary is assumed to be free to drain.
123
Fig. 1 Computational domain and boundary conditions on the left picture, and example of the deformation
obtained after applying the uniform load on the right picture
For the simulations we consider the stabilized P1 -P1 scheme, which has been imple-
mented in the HAZmath library [58]. As an example of the solution of this benchmark,
we show in the right side of Fig. 1 the deformation obtained after applying the uniform
load.
Regarding the considered material properties, in this numerical experiment, we fix the
Young’s modulus E = 3 × 104 N/m2 , whereas the Poisson’s ratio ν and the permeability K
will be varied in some of the experiments. We will consider a diagonal permeability tensor
K = K I, with K constant.
We start by analyzing the performance of the proposed preconditioners with respect to
the discretization parameters, so that we consider different values of the mesh size h and
the time step size τ . In particular, h will vary from 1/4 to 1/32, whereas τ will range
from 0.0001 to 0.1. The physical parameters are considered fixed as K = 10−6 m2 and
ν = 0.2.
We consider flexible GMRES as the outer iteration with a relative residual stopping criteria
of 10−6 , which is preconditioned with the block diagonal and triangular proposed precon-
D , B
ditioners B D , B L , BU , as well as, with their inexact counterparts B L , B
U . When these
latter are considered, the diagonal blocks are solved inexactly by preconditioned GMRES
with a tolerance of 10−2 . In Table 1, we show the number of GMRES iterations needed to
fulfill the stopping criterion, and we can observe that the proposed preconditioners are robust
with respect to the discretization parameters h and τ .
Next, we investigate the robustness of the block preconditioners with respect to the physical
parameters K and ν. Thus, we fix the mesh size as h = 1/16 and the time step size as τ = 0.01,
and vary the values of the permeability K and the Poisson’s ratio ν. In Table 2, we show the
number of GMRES iterations needed to achieve the stopping criterion, for the different exact
and inexact preconditioners. The results corresponding to the variation of the permeability
are shown on the left part of the table, where K varies from 1 to 10−10 while ν is fixed as 0.2.
On the right part of the table, we fix K = 10−6 and we consider different values of ν from
0.1 to 0.499, which already is close to the incompressible material case. From the iteration
counts in the table, we can see that the proposed preconditioners are also quite robust with
respect to the physical parameters.
123
Table 1 Two-field stabilized BD BL BU

P1 -P1 scheme: iteration counts
τ /h 1 1 1 1 1 1 1 1 1 1 1 1
for the block preconditioners for 4 8 16 32 4 8 16 32 4 8 16 32
different values of the
discretization parameters (∗ 0.1 7 7 8 * 5 5 6 * 4 4 4 *
means the direct method for 0.01 7 7 8 * 5 5 6 * 4 4 5 *
solving diagonal blocks is out of
0.001 7 7 8 * 5 5 6 * 5 5 6 *
memory)
0.0001 7 7 8 * 5 5 6 * 5 5 6 *
D
B L
B U
B
τ /h 1 1 1 1 1 1 1 1 1 1 1 1
4 8 16 32 4 8 16 32 4 8 16 32
0.1 8 8 9 9 6 6 8 8 6 6 8 8
0.01 8 8 9 9 6 6 8 8 6 6 8 8
0.001 8 8 9 9 6 6 8 8 6 6 8 8
0.0001 8 8 9 9 7 6 8 8 6 7 8 8
Table 2 Two-field stabilized P1 -P1 scheme: iteration counts when varying physical parameters K or ν
ν = 0.2 and varying K K = 10−6 and varying ν

1 10−2 10−4 10−6 10−8 10−10 0.1 0.2 0.4 0.45 0.49 0.499
BD 4 7 8 8 8 8 7 8 11 11 12 12
BL 2 5 6 6 6 6 5 6 8 8 8 9
BU 3 4 5 5 5 5 4 5 6 6 5 4
D
B 5 8 9 9 9 9 8 9 12 13 14 13
L
B 5 7 8 8 8 8 7 8 11 11 12 12
U
B 5 7 8 8 9 8 7 8 7 8 17 11
3.2 Three-field formulation
In many of the applications where the coupling of fluid flow and mechanical deformation
within a porous media is important, the flow of the fluid through the medium is of primary
interest. Although from the reduced displacement–pressure formulation the fluid flux can be
recovered, a natural approach is to introduce this value as a primary variable instead, so that we
deal with the classical three-field formulation of the Biot’s model. The incorporation of this
extra unknown can be seen as a disadvantage against the two-field formulation, regarding the
computational cost, but there appear other advantages of this approach. The direct calculation
of the fluid flux, instead of its calculation by post-processing allows a higher order of accuracy
in its computation, and in addition, the mass conservation for the fluid phase is ensured by
using continuous elements for the fluid flux variable.
Since the three-field formulation considered here has the displacement u, the pore pressure
p and the Darcy’s velocity w as primary unknowns, let us define proper Sobolev spaces
V ⊂ H 1 (), Q ⊂ H 1 () and W ⊂ H(div, ), where H(div, ) contains the square
integrable vector-valued functions with square integrable divergence. Then, the variational
formulation for Biot’s three-field consolidation model reads:
For each t ∈ (0, T ], find (u(t), p(t), w(t)) ∈ V × Q × W such that
a(u, v) − α( p, div v) = (ρ g, v), ∀ v ∈ V , (30)
123

1 ∂p ∂u
, q + α div , q + (div w, q) = ( f , q), ∀ q ∈ Q, (31)
Mb ∂t ∂t
(K −1 μ f w, r) − ( p, div r) = (ρ f g, r), ∀ r ∈ W , (32)
where a(u, v) corresponds to the bilinear form representing linear elasticity, already intro-
duced in Sect. 3.1.
In order to consider the finite element approximation of problem (30)–(32), given a par-
tition of the domain into d-dimensional simplices, Th , we associate a triple of piecewise
polynomial, finite-dimensional spaces,
Vh ⊂ V, Q h ⊂ Q, W h ⊂ W . (33)
By using backward Euler as a time discretization on a time interval (0, T ] with constant
time-step size τ , the discrete scheme corresponding to the three-field formulation (30)–(32)
reads:
Find (umh , ph , w h ) ∈ V h × Q h × W h such that
m m
a(umh , v h ) − α( ph , div v h ) = (ρ g, v h ), ∀ v h ∈ V h ,
m
(34)
1 m

ph , qh + α div um h , qh + τ (div w h , qh ) = ( f , qh ), ∀ qh ∈ Q h ,
m
(35)
Mb
τ (K −1 μ f wm h , r h ) − τ ( ph , div r h ) = τ (ρ f g, r h ), ∀ r h ∈ W h ,
m
(36)

where ( f , qh ) = τ ( f , qh ) + M1b phm−1 , qh + α div um−1 h , qh , and (um
h , ph , w h ) is an
m m
approximation to (u(·, tm ), p(·, tm ), w(·, tm )) , at time tm = mτ, m = 1, 2, . . . . Note that

the last equation has been scaled by τ for symmetry.
Moreover, this discrete variational form can be represented in block matrix form,
⎛ ⎞
⎛ ⎞ Au α BuT 0
uh ⎜ ⎟
A ⎝ ph ⎠ = b, with A = ⎜ ⎝ −α Bu Mb M p −τ Bw ⎠ ,
1 ⎟ (37)
wh
0 τ BwT τ Mw
where u, p, and w are the unknown vectors for the displacement, the pressure, and the Darcy
velocity, respectively. The blocks in the definition of matrix A correspond to the following
bilinear forms:
a(uh , v h ) → A u , −(div uh , qh ) → Bu ,
−(div w h , qh ) → Bw , (K −1 μ f w h , r h ) → Mw , ( ph , qh ) → M p .
3.2.1 Stabilized P1 −RT−P0
Due to its application to existing reservoir engineering simulators, one of the most frequently
considered schemes for Biot’s model is a three-field formulation based on piecewise linear
elements for displacements, Raviart–Thomas–Nédélec elements for the fluid flux, and piece-
wise constants for the pressure, that is, the so-called P1 –RT–P0 scheme. This finite element
scheme, however, does not satisfy an inf-sup condition uniformly with respect to the physical
parameters of the problem, and therefore, in [20] we proposed a stabilization of this popular
element which gives rise to uniform error bounds, which is the discretization that we consider
here.
123
More specifically, to discretize displacements we choose a piecewise linear continuous

finite-element space, V h,1 , enriched with edge/face (2D/3D) bubble functions, V b , that is
V h = V h,1 ⊕V b (see [65, pp. 145–149]), for Darcy’s velocity W h is the lowest order Raviart–
Thomas–Nédélec space (RT0), that is, W h = {w h ∈ W | w h |T = a+ηx, a ∈ Rd , ∀T ∈ Th },
and we use a piecewise constant space (P0 ) for Q h in order to approximate the pressure field.
Thus, any uh ∈ V h can be written as uh = ulh + ubh , with ulh ∈ V h,1 and ubh ∈ V b , and we
can rewrite system (37) as follows,
⎛ b⎞ ⎛ ⎞
uh Abb Abl α BbT 0
⎜ ul ⎟ ⎜ T ⎟
⎜ ⎟ ⎜ Abl All α BlT 0 ⎟
A ⎜ h ⎟ = b, with A = ⎜ ⎜ −α B
⎟, (38)
⎝ ph ⎠ ⎝ b −α Bl M1b M p −τ Bw ⎟ ⎠
wh 0 0 τ BwT τ Mw
where ub , ul , p, and w are the unknown vectors for the bubble components of the displace-
ment, the piecewise linear components of the displacement, the pressure, and the Darcy
velocity, respectively. New blocks in the definition of matrix A correspond to the following
bilinear forms:
a(ubh , v bh ) → Abb , a(ulh , v bh ) → Abl , a(ulh , vlh ) → All ,
−(div ubh , qh ) → Bb , −(div ulh , qh ) → Bl ,

Abb Abl
and are related to those in (37) by A u = and Bu = Bb Bl .
Alb All
The well-posedness of the discretized system provides a convenient framework with which
to construct block preconditioners, and in order to obtain parameter-robust preconditioners,
such well-posedness property must hold independently to the discretization and physical
parameters. In order to demonstrate such a uniform stability of the discrete system, we
introduce the following weighted norm, for x = (u, p, w)T ∈ X = V h × Q h × W h ,
1/2
|||x||| := u2A u + c−1 p p 2
Mp + τ w 2
Mw + τ 2
c p Bw w 2
−1 , (39)
M p
2 −1 √
where c p := αζ 2 + M1b with ζ := λ + 2μ/d as before, and d = 2 or 3 is the dimension
of the problem. Next, for x = (u, p, w), y = (v, q, r), we introduce the following composite
bilinear form
L(x, y) := a(u, v) − (α p, div v) + τ (K −1 μ f w, r) − τ ( p, div r)

1
− p, q − (α div u, q) − τ (div w, q). (40)
Mb
Then, in [20], it was shown that under certain conditions (referred to as Stokes–Biot stability
[20, Def. 3.1]) on the space X, the block matrix form A defined in (38) is well-posed with
respect to the weighted norm (39).
Theorem 4 If the triple (V h , W h , Q h ) is Stokes–Biot stable, that is,

• a(u, v) ≤ C V u1 v1 , for all u, v ∈ V h ,
• a(u, v) ≥ αV u21 , for all u ∈ V h ,
• The pair of spaces (W h , Q h ) is Poisson stable, i.e., it satisfies stability and continuity
conditions required by the mixed discretization of the Poisson equation,
• The pair of spaces (V h , Q h ) is Stokes stable,
123
then, it is fulfilled that L(·, ·) is continuous with respect to |||·|||, and the following inf-sup
condition holds,
L(x, y)
inf sup ≥ γ, (41)
0= x∈X 0 = y∈X |||x|||||| y|||
where the constants involved are independent of mesh size h, time step size τ , and the physical
parameters.
Proof The proof follows directly from Case I in the proof found in [21, Theorem 6], although
the weighted norm used here is slightly different from the norm used in [21].

Next, we use the well-posedness property to develop block preconditioners for A. Fol-
lowing the general framework previously introduced, we first consider block diagonal
preconditioners (norm-equivalent preconditioners) and then we discuss block triangular
(upper and lower) preconditioners following the framework developed for Field-of-Value
(FOV) equivalent preconditioners. For both cases, the theoretical bounds on their perfor-
mance remain independent of the discretization and physical parameters of the problem.
Since system (38) satisfies inf-sup condition (41), based on the framework proposed in
Sect. 2.1, a natural choice for a norm-equivalent preconditioner is the Riesz operator with
respect to the inner product corresponding to weighted norm (39), that is,
⎛ ⎞−1
Au 0 0
⎜ 0 α2
+ M1b M p ⎟
BD = ⎜ ⎝ ζ2
0 ⎟ ,
⎠ (42)
2 −1
α
0 0 τ M w + τ 2 ζ 2 + Mb
1
Aw
where Aw := BwT M −1p Bw .

In order to avoid the need of inverting exactly the diagonal blocks, we also propose the
following inexact preconditioner,
⎛ ⎞
Su 0 0
D = ⎝ 0 S p 0 ⎠ .
B (43)
0 0 Sw
which is obtained by replacing the diagonal blocks in (42) by their spectrally-equivalent
symmetric and positive definite approximations, Su , S p and Sw , satisfying that
c1,u (Su u, u) ≤ (A−1 u u, u) ≤ c2,u (Su u, u), (44)

−1
α2 1 −1
c1, p (S p p, p) ≤ + M p p, p ≤ c2, p (S p p, p), (45)
ζ2 Mb
⎛ ⎞
2 −1 −1
α 1
c1,w (Sw w, w) ≤ ⎝ τ Mw + τ 2 + Aw w, w⎠≤ c2,w (Sw w, w), (46)
ζ2 Mb
where the involved constants c1,u , c1, p , c1,w , c2,u , c2, p , and c2,w are independent of dis-
cretization and physical parameters.
In practice, naturally Su can be defined by standard multigrid methods, and S p is obtained
2 −1
by a diagonal scaling such that S p = αζ 2 + M1b M −1 p since when using P0 elements
123
M p results to be diagonal. The choice of Sw , however, may require special attention since
2 −1
for large values of τ matrix τ Mw + τ 2 αζ 2 + M1b Aw is numerically close to singular
and then special preconditioners are needed. In particular, Sw can be defined by either an
HX-preconditioner (Auxiliary Space Preconditioner) [53, 66] or multigrid with special block
smoothers [52].
Following the FOV-equivalent framework previously introduced in Sect. 2.2, next, we
consider more general preconditioners. In particular, block triangular preconditioners are
designed for the linear system given in (38).
We first consider the block lower triangular preconditioner,
⎛ ⎞−1
Au 0 0
⎜−α Bu α2
+ 1 ⎟
BL = ⎜
⎝ ζ2 Mb Mp 0 ⎟
⎠ , (47)
−1
α2
0 τ BwT τ Mw + τ2 ζ2
+ 1
Mb Aw
for which we can prove the FOV-equivalence property in (8), that is,
Theorem 5 Assuming a shape regular mesh and the discretization described above, there
exist constants and ϒ, independent of discretization and physical parameters, such that,
for any x = (u, p, w)T = 0,
(B L A x, x)(B D )−1 B L A x(B D )−1

≤ , ≤ ϒ.
(x, x)(B D )− 1 x(B D )−1
Analogously to the two-field formulation case, we also consider an inexact approach to

the lower triangular preconditioner given as follows,
⎛ ⎞−1
Su−1 0 0
L = ⎝ −α Bu
B S −1
p 0 ⎠ , (48)
0 τ BwT Sw−1
where Su , S p and Sw , satisfy the spectral equivalence relations (44)–(46) and can be chosen
in practice as we have specified in the previous section. Thus, B L also satisfies the FOV-
equivalence property, as stated in the following theorem.
2 −1
Theorem 6 If I − Su A u A u ≤ ρ ≤ 0.2228, and S p = αζ 2 + M1b M −1
p , then there exists
constants and ϒ, independent of discretization and physical parameters, such that, for
any x = (u, p, w)T = 0,

L A x, x −1
B L A x(B )−1
B
(B D )
≤ , D
≤ ϒ.
(x, x)(BD )− 1 x(BD )−1
The proofs of the above two theorems turn out to be a special case of some of the proofs
included in [47], and thus are omitted here. With these results, we can use the theoreti-
cal framework introduced in Sect. 2.2 and then B L and B L result to be parameter-robust
preconditioners for the GMRES method.
123
Fig. 2 2D physical and

computational domain for
Mandel’s problem
Similar arguments can also be applied for the design of block upper triangular precondi-
tioners, so that we can consider
⎛ ⎞−1
Au α BuT 0
⎜ 0 α2
+ M1b M p −τ Bw ⎟
BU = ⎜ ⎝ ζ2 ⎟ ,
⎠ (49)
2 −1
α
0 v0 τ M w + τ ζ 2 + Mb
2 1
Aw
and the corresponding inexact preconditioner,

⎛ −1 ⎞−1
Su α BuT 0
U = ⎝ 0
B S −1
p −τ Bw ⎠ . (50)
0 0 Sw−1
Analogous theorems to Theorems 5 and 6 can be stated, so that the FOV-equivalence property
U can be proved, giving rise to parameter robustness for the block upper triangular
for BU and B
preconditioners.
Numerical experiment: Mandel’s problem
We consider Mandel’s problem in two-dimensions, which models an infinitely long saturated
porous slab sandwiched between a top and bottom rigid frictionless plate. At time t = 0,
each plate is loaded with a constant vertical force of magnitude 2F per unit length as shown
in Fig. 2.
This problem is an important benchmarking tool as the analytical solution is known [67],
and therefore, it has been widely used in the literature to test numerical methods for poroelastic
problems. The analytical solution for pressure is given by
∞
sin αn αn x −αn2 ct
p(x, y, t) = 2 p0 cos − cos αn exp , (51)
αn − sin αn cos αn a a2
n=1
where p0 = B(1 + νu )F, B = 1 is the Skempton’s coefficient, νu = 3ν+B(1−2ν)

1
3a 3−B(1−2ν) is the
undrained Poisson ratio, c is the consolidation coefficient given by c = K (λ + 2μ), and αn
are the positive roots to the nonlinear equation,
1−ν
tan αn = αn .
νu − ν
Due to symmetry of the problem we only need to solve in the top right quadrant, defined
as = (0, 1) × (0, 1). We cover with a uniform triangular grid by dividing an N × N
uniform square grid into right triangles, where the mesh spacing is defined by h = N1 . For the
material properties, μ f = 1, α = 1, and Mb = 106 , the Lamé coefficients are computed in
terms of the Young’s modulus, which is fixed to be E = 104 , and the Poisson’s ratio, ν. This
123
Table 3 Three-field stabilized P1 −RT−P0 scheme: iteration counts for the block preconditioners with varying
discretization parameters
BD BL BU
τ /h 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
8 16 32 64 128 8 16 32 64 128 8 16 32 64 128
0.1 39 40 40 40 38 19 19 18 17 17 19 19 19 18 17
0.01 26 34 39 39 38 15 18 19 18 17 14 17 18 18 17
0.001 23 23 28 34 37 11 12 15 17 18 10 11 14 17 17
0.0001 21 21 21 21 21 11 10 10 13 15 8 9 9 12 14
D
B L
B U
B
τ /h 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
8 16 32 64 128 8 16 32 64 128 8 16 32 64 128
0.1 39 40 40 40 36 19 20 19 19 18 19 19 19 18 20
0.01 26 34 39 39 38 15 18 19 19 18 14 17 18 18 17
0.001 23 23 23 34 37 11 13 15 17 18 10 12 15 17 17
0.0001 21 22 21 23 29 11 11 11 13 15 9 9 10 12 15
latter, as well as the permeability tensor, K , will be varied in order to show the robustness
of the proposed preconditioners with respect to the physical parameters. In all test cases, we
consider a diagonal permeability tensor K = K I, with constant K .
Next, we test the performance of the proposed block preconditioners, B D , B L , BU , B D ,

B L , and BU , with respect to discretization and physical parameters. For each test we use
flexible GMRES to solve the linear system obtained from the bubble-enriched P1 −RT0−P0
discretization. A stopping tolerance of 10−8 was used for the relative residual of the linear
system, measured relative to the norm of the right hand side. The exact solves for the blocks
in B D , B L , and BU are done using the UMFPACK library [68–71]. For the inexact precon-
ditioners, Su is inverted using GMRES preconditioned with unsmoothed aggregation AMG
in a V-cycle, solved to a relative residual tolerance of 10−3 , the Sw block is solved using an
auxiliary space preconditioned GMRES to a relative residual tolerance of 10−3 [52, 53, 66],
and the action of S p is directly computed, since using a piecewise constant finite-element
space for pressure results in a diagonal matrix for M p .
First, we study the robustness of the preconditioners with respect to the discretization
parameters. Thus, next test covers different values for the mesh size and the time step size,
while the physical parameters are fixed as ν = 0.0 and K = 10−6 . Then, in Table 3,
we show the number of preconditioned GMRES iterations needed to solve the full bubble
system up the required tolerance, by using the proposed exact and inexact block precondi-
tioners and for different mesh sizes and time-step sizes. We can observe from the relatively
consistent iteration counts the robustness of the preconditioners with respect to the dis-
cretization parameters. Since the block upper and lower triangular preconditioners contain
more coupling information than the block diagonal preconditioners, as a result we see that
they preform better than the block diagonal preconditioners. Finally, it is also important to
notice that the performance impact of using the inexact block preconditioners is negligible
versus using the exact block preconditioners, which implies that the inexact precondition-
ers could potentially be solved with less strict tolerance, resulting in more computational
efficiency.
Next, we analyze the robustness of the proposed preconditioners with respect to the phys-
ical parameters. In order to do that, the mesh size is fixed to h = 128
1
and the time-step size is
123
Table 4 Three-field stabilized P1 −RT−P0 scheme: iteration counts for the block preconditioners with varying
physical parameters K and ν

1 10−2 10−4 10−6 10−8 10−10 0.1 0.2 0.4 0.45 0.49 0.499
BD 23 25 35 38 29 19 45 52 39 36 28 20
BL 7 11 15 17 15 9 16 19 11 11 9 10
BU 13 16 17 16 15 7 20 22 16 14 11 16
D
B 35 33 36 38 29 19 45 52 39 26 23 17
L
B 14 15 16 18 15 10 17 20 14 12 11 12
U
B 27 22 17 17 15 8 21 24 17 16 10 16
assumed to be τ = 0.01, whereas the physical values of ν and K are varied. Then, in Table 4
we show the corresponding iteration counts for the block preconditioned GMRES.
On the left part of the table, we present the results obtained when the Poisson’s ratio
is fixed as ν = 0.0 and the permeability varies from 1 to 10−10 , and on the right part
of the table the permeability is fixed as K = 10−6 whereas different values for the Pois-
son’s ratio are considered, even approaching the incompressible material case (ν = 0.5).
Again, we observe robustness of the proposed preconditioners, this time with respect to
the physical parameters. In addition, we can see that the use of the inexact precondition-
ers has minimal impact on their performance. Finally, notice that an interesting result is
the better performance when the system is approaching the two considered limit cases,
namely in the limit of impermeability (K → 0) and in the limit of incompressibility
(ν → 0.5).
Remark. Bubble eliminated system
A noteworthy result of [20] is that the degrees of freedom added to the faces to obtain the
stabilized scheme can be eliminated resulting in a stable scheme with the same number
of unknowns as in the initial P1 −RT0−P0 discretization. In particular, one can replace
the enrichment bubble block, Abb , in (37) with a spectrally equivalent diagonal matrix,
Dbb := (d + 1)diag(Abb ), obtaining the following linear operator,
⎛ ⎞
Dbb Abl α BbT 0
⎜ T ⎟
⎜ Abl All α BlT 0 ⎟
⎜ ⎟
AD = ⎜ ⎟, (52)
⎜ −α Bb −α Bl 1
−τ Bw ⎟
⎝ Mb M p ⎠
0 0 τ BwT τ Mw
where the stabilization term can be eliminated in a straightforward way by static condensation,
yielding,
⎛ T D −1 A T D −1 B T ⎞
All − Abl bb bl α BlT − α Abl bb b 0
⎜ ⎟
AE = ⎜ −1 −1 T
⎝ −α Bl + α Bb Dbb Abl Mb M p + α Bb Dbb Bb −τ Bw ⎠ .
1 2 ⎟ (53)
0 τ BwT τ Mw
In such a way, an optimal stable discretization with the lowest possible number of degrees
of freedom, equivalent to a discretization with P1 −RT0−P0 elements, is obtained.
123
Due to the spectral equivalence between Abb and Dbb , the formulation (52) is still well-
posed (see [20]), and remains well-posed independently of the physical and discretization
parameters. Since (37) is indefinite, the well-posedness of (53) does not simply follow directly,
but it was proved in [47], where this well-posedness allowed the design of block precondi-
tioners for the “bubble-eliminated system” given in (53).
3.2.2 Non conforming CR−RT−P0
Finally, we consider a stable finite element method for the three-field formulation of Biot’s
model that was proposed in [19] to avoid non-physical oscillations in the pressure field
by using the lowest possible approximation order. This scheme considers a nonconforming
finite element space for the displacements u and conforming finite element spaces for both
the Darcy velocity w and the pressure p. More concretely, as in Sect. 3.2.1, W h is the lowest
order Raviart–Thomas–Nédélec space [72–74], and Q h is the space of piecewise constant
functions with respect to the triangulation Th , that is, P0 . Regarding the displacements, V h
is chosen as the Crouzeix–Raviart finite element space [75], which consists of vector valued
functions which are linear on every element T ∈ Th and satisfy the following continuity
conditions

2
V h = v h ∈ L () [[v h ]]e = 0, for all e ∈ E , o
e
where [[v h ]]e denotes the jump across an interior face e ∈ E o for a given function v h .
Then (um h , ph , w h ) ∈ V h × Q h × W h must satisfy the discrete formulation in (34)–(36).
m m
However, if we simply take ah (·, ·) = a(·, ·), as done before, Korn’s inequality may fail for
the standard discretization by Crouzeix–Raviart elements without additional stabilization,
and, therefore, ah (·, ·) is not coercive. It is also possible that Korn’s inequality hold, but
the coercivity constant will approach infinity, blowing up, as the mesh size h approaches
zero. For more details on nonconforming linear elements for elasticity problems and discrete
Korn’s inequality, we refer to [76, 77].
In order to fix such a problem, we can consider the following perturbation of the bilinear
form which does satisfy the Korn’s inequality and was proposed by Hansbo and Larson [78],

ah (v, w) = a(v, w) + aJ (v, w), where aJ (v, w) = 2μγ1 h −1
e [[v]]e [[w]]e ,
e∈E e
where E denotes the set of faces (interfaces) in the triangulation Th , and constant γ1 > 0
is a fixed real number away from 0 (i.e. γ1 = 21 is an acceptable choice). As shown in
Hansbo and Larson [78] the bilinear form ah (·, ·) is positive definite and the corresponding
error is of optimal (first) order in the corresponding energy norm. Moreover, the resulting
method is locking free and therefore we use such ah (·, ·) in our nonconforming scheme (see
more details in [19]. Regarding the bilinear form corresponding to the term (K −1 μ f wm h , rh)
in (36), the standard choice is just taking the usual L 2 () inner product, leading to a mass
matrix in the Raviart–Thomas–Nédélec element. Here, however, we use mass lumping in
the Raviart–Thomas space, since it actually gives an oscillation-free approximation while
maintains optimal error estimates (see [19]). Notice that such lumped mass approximation
results in a block diagonal matrix and, therefore, it allows to eliminate the Darcy velocity and
reduce the three-field formulation to a two-field formulation involving only displacement and
pressure, reducing the size of the linear system to solve at each time step and consequently
saving computational cost.
123
Considering all these modifications with respect to the three-field discrete formulation
introduced in the previous sections, we can define the corresponding composite bilinear
form L(x, y), as done in (40), and⎛ prove theT well-posedness
⎞ of the linear system A x = b,
A u α Bu 0
with x = (u, p, w)T and A = ⎝−α Bu 0 −τ Bw ⎠ , with respect to the weighted norm
0 τ BwT τ Iw
introduced in (39). Notice that we denote as Iw to the lumped mass matrix for Darcy’s
velocity, and that the first term in the flow equation is neglected by assuming that the fluid
is incompressible, as we did in the two-field formulation section, so that we consider this
challenging case also in the three-field formulation framework.
Theorem 7 If the triple (V h , W h , Q h ) is Stokes–Biot stable, then, it is fulfilled that L(·, ·) is
continuous with respect to |||·|||, and the following inf-sup condition holds,
L(x, y)
inf sup ≥ γ, (54)
0= x∈X 0 = y∈X |||x|||||| y|||
where the constants involved are independent of mesh size h, time step size τ , and the physical
parameters.
Proof The proof is as the well-posedness proof in [19], but considering the weighted norm
(39). We do not include the proof here since it uses the same techniques than those applied
in the previous proofs included in this work.

After proving the well-posedness of the system, we can consider again as norm-equivalent
preconditioner the Riesz operator with respect to the inner product associated with the
weighted norm introduced in (39), that is,
⎛ ⎞−1
Au 0 0
⎜ α2 ⎟
BD = ⎝ 0 ζ2
Mp 0 ⎠ . (55)
τ Iw + τ 2 αζ 2 Aw
2
0 0
This results in a parameter-robust block diagonal preconditioner for the solution of the con-
sidered Biot’s three-field formulation.
By considering spectrally-equivalent SPD approximations of the diagonal blocks, we can
construct the following inexact block diagonal preconditioner,
⎛ ⎞
Su 0 0
D = ⎝ 0 S p 0 ⎠ ,
B (56)
0 0 Sw
where Su , S p and Sw satisfy that
c1,u (Su u, u) ≤ (A−1 u u, u) ≤ c2,u (Su u, u), (57)
2
ζ −1
c1, p (S p p, p) ≤ M p, p ≤ c2, p (S p p, p), (58)
α2 p

2 −1
2ζ
c1,w (Sw w, w) ≤ τ Iw + τ 2Aw w, w ≤ c2,w (Sw w, w), (59)
α
where the involved constants are independent of the discretization and physical parameters.
As previously commented, again, Su can be defined by multigrid methods, S p is applied as
123
a Jacobi iteration since the original block results to be diagonal, and Sw can be defined by an
auxiliary space method as the HX-preconditioner.
Next, we consider block triangular preconditioners, arising from the FOV-equivalent frame-
work introduced in Sect. 2.2. In particular, the following block lower triangular preconditioner
can be defined,
⎛ ⎞−1
Au 0 0
⎜−α B α 2 ⎟
BL = ⎜
⎝
u
λ+ 2μd
Mp

0

⎟ ,
⎠ (60)
τ2 2μ
0 τ Bw T τ Iw + α 2 λ + d A w
as well as its inexact counterpart given by

⎛ −1 ⎞−1
Su 0 0
BL = ⎝−α Bu S −1
p 0 ⎠ , (61)
0 τ BwT Sw−1 .
where Su , S p and Sw satisfy the spectral equivalence relations (57)–(59). These operators B L
L satisfy that they are FOV-equivalent with A, since we can prove analogous theorems
and B
to Theorems 5 and 6, and therefore they result in parameter-robust preconditioners for the
GMRES method.
Analogously, we can also consider uniform right preconditioners given by the following
block upper triangular operators,
⎛ ⎞−1
A u α BuT 0 ⎛ −1 ⎞−1
⎜0 α M 2 ⎟ Su α BuT 0
BU = ⎜ −τ Bw ⎟ , B U = ⎝ 0 S −1 ⎠ .
p −τ Bw
p
⎝ λ+ 2μ
d ⎠
−1
τ Iw + τ λ + 2μ Aw
2
0 0 0 0 Sw
α2 d
(62)
Numerical experiment: poroelastic problem on a square domain with a uniform load

We consider a Terzaghi’s type problem which models a column of poroelastic material
= (0, 1) × (0, 1) which drains on the north (top) edge of the boundary, where also a
uniform unit load is applied, that is,
p = 0, σ · n = (0, −1)t , on 1 = [0, 1] × 1.
The medium is assumed to be rigid and impermeable on the rest of the boundary, namely,
the rest of the boundary conditions are:
∇ p · n = 0, u = 0, on 2 = ∂ \ 1 .
For clarity, the computational domain together with the prescribed boundary conditions are
shown in Fig. 3.
Next, we test the performance of the proposed exact, B D , B L and BU , and inexact, B D ,
L and B
B U , preconditioners with respect to different values of both the discretization and
physical parameters. In order to do that, we consider the preconditioned flexible GMRES
with a relative residual stopping criteria of 10−6 . Regarding the material properties, in this
numerical experiment, we fix the Young’s modulus E = 3 × 104 N/m2 , and we will consider
a diagonal permeability tensor K = K I, with K constant.
123
Fig. 3 Computational domain

and boundary conditions
corresponding to the poroelastic
test problem on a square domain
with a uniform load
Table 5 Three-field CR−RT−P0 BD BL BU

scheme: Iteration counts for the
τ /h 1 1 1 1 1 1 1 1 1 1 1 1
block preconditioners for 8 16 32 64 8 16 32 64 8 16 32 64
different values of the
discretization parameters 0.1 42 46 48 52 22 23 24 25 23 24 26 27
0.01 41 43 47 47 22 22 23 25 20 21 23 24
0.001 41 43 43 47 22 22 23 24 18 18 21 22
0.0001 41 43 43 46 22 22 23 24 18 18 19 19
D
B L
B U
B
τ /h 1 1 1 1 1 1 1 1 1 1 1 1
8 16 32 64 8 16 32 64 8 16 32 64
0.1 42 48 54 57 24 26 26 31 26 30 29 32
0.01 42 44 49 50 24 25 27 28 24 26 26 28
0.001 42 44 49 49 24 25 27 28 25 24 24 26
0.0001 42 46 49 50 24 25 27 28 25 29 30 30
We start by testing the robustness of the proposed preconditioners with respect to the
discretization parameters. With this purpose, we fix the values of the permeability as K =
10−6 and the Poisson’s ratio as ν = 0.2, and vary the spatial discretization parameter h
from 1/8 to 1/64 and the time step size τ from 0.1 to 0.0001. For these values, in Table 5,
we show the number of preconditioned GMRES iterations that are necessary to reduce the
residual under the prescribed tolerance. From the iteration counts shown in the table, we can
observe a robust behavior or the proposed preconditioners with respect to the discretization
parameters. Also, we can derive similar conclusions to those mentioned in the numerical
examples for the other considered numerical schemes, as for example, the good performance
provided by the inexact counterparts of the proposed block preconditioners.
Finally, we want to show the robustness of the proposed preconditioned GMRES method
with respect to the material parameters. With this aim, in Table 6, we provide iteration counts
for the proposed block preconditioners for varying physical parameters. In particular, on the
left part of Table 6, we assume a fixed value of the Poisson’s ratio ν = 0.2 and we choose the
permeability K ranging from 1 to 10−10 , whereas on the right part of such a table, we fix the
value of the permeability as K = 10−6 and vary the value of ν from 0.1 to 0.499, approaching
the incompressible material case. In both cases, we observe again a robust behavior of the
proposed preconditioners.
123
Table 6 Three-field CR−RT−P0 scheme: iteration counts for the block preconditioners with varying physical
parameters K and ν

1 10−2 10−4 10−6 10−8 10−10 0.1 0.2 0.4 0.45 0.49 0.499
BD 80 61 47 47 47 47 51 47 30 24 17 13
BL 54 41 23 23 23 23 26 23 15 12 9 7
BU 52 41 23 23 23 23 26 23 16 13 9 6
D
B 98 63 49 49 49 49 53 49 35 29 24 21
L
B 76 45 27 27 27 27 30 27 19 16 14 13
U
B 68 44 30 30 30 30 33 30 22 17 18 20
4 Conclusions
In this work, we have presented an overview about the design of parameter-robust precondi-
tioners for poromechanics problems. In particular we have considered different conforming
and non-conforming discretizations proposed by the authors for the two- and three-field
formulations of the Biot’s model and we have proposed block diagonal and triangular precon-
ditioners with a uniform performance with respect to discretization and physical parameters.
We have explained how the stability and well-posedness of the discrete problem provides a
foundation for designing such robust preconditioners which yield uniform convergence rates
for Krylov subspace methods. We present theoretical results confirming the robustness of the
proposed preconditioners with respect to mesh size, time step, and the physical parameters of
the model, and finally we demonstrate this behavior with numerical experiments in which we
consider well-known benchmark problems that are used in the literature for testing numerical
methods for poroelastic problems.
Acknowledgements The work of Carmen Rodrigo is supported in part by the Spanish project PGC2018-
099536-A-I00 (MCIU/AEI/FEDER, UE). Francisco J. Gaspar has received funding from the Spanish project
PID2019-105574GB-I00 (MCIU/AEI/FEDER, UE). Carmen Rodrigo and Francisco J. Gaspar acknowledge
support from the Diputación General de Aragón, Spain (Grupo de referencia APEDIF, ref. E24_17R). The
work of Adler, Hu and Ohm is partially supported by the National Science Foundation (NSF) under Grant
DMS-2208267. The research of Zikatanov is supported in part by the U. S.-Norway Fulbright Foundation and
the U. S. National Science Foundation Grant DMS-2208249.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Data availability Not applicable.
Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence,
and indicate if changes were made. The images or other third party material in this article are included in the
article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is
not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
123
References
1. Terzaghi, K.: Theoretical Soil Mechanics. Wiley, New York (1943). https://doi.org/10.1002/
9780470172766
2. Biot, M.A.: General theory of three-dimensional consolidation. J. Appl. Phys. 12(2), 155–164 (1941).
https://doi.org/10.1063/1.1712886
3. Biot, M.A.: Theory of elasticity and consolidation for a porous anisotropic solid. J. Appl. Phys. 26(2),
182–185 (1955). https://doi.org/10.1063/1.1721956
4. Lee, J.J., Mardal, K.-A., Winther, R.: Parameter-robust discretization and preconditioning of Biot’s con-
solidation model. SIAM J. Sci. Comput. 39(1), 1–24 (2017). https://doi.org/10.1137/15M1029473
5. Wang, H.F.: Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology.
Princeton University Press, Princeton (2001). https://doi.org/10.1515/9781400885688
6. Ben-Hatira, F., Saidane, K., Mrabet, A.: A finite element modeling of the human lumbar unit including
the spinal cord. J. Biomed. Sci. Eng. 5, 146–152 (2012). https://doi.org/10.4236/jbise.2012.53019
7. Smith, J.H., Humphrey, J.A.: Interstitial transport and transvascular fluid exchange during infusion into
brain and tumor tissue. Microvasc. Res. 73(1), 58–73 (2007). https://doi.org/10.1016/j.mvr.2006.07.001
8. Støverud, K.H., Alnæs, M., Langtangen, H.P., Haughton, V., Mardal, K.-A.: Poro-elastic modeling of
syringomyelia—a systematic study of the effects of pia mater, central canal, median fissure, white and gray
matter on pressure wave propagation and fluid movement within the cervical spinal cord. Comput. Methods
Biomech. Biomed. Eng. 19(6), 686–698 (2016). https://doi.org/10.1080/10255842.2015.1058927
9. Axelsson, O., Blaheta, R., Byczanski, P.: Stable discretization of poroelasticity problems and efficient
preconditioners for arising saddle point type matrices. Comput. Visual Sci. 15, 191–207 (2012). https://
doi.org/10.1007/s00791-013-0209-0
10. Gaspar, F.J., Lisbona, F.J., Vabishchevich, P.N.: A finite difference analysis of Biot’s consolidation model.
Appl. Numer. Math. 44(4), 487–506 (2003). https://doi.org/10.1016/S0168-9274(02)00190-3
11. Gaspar, F.J., Lisbona, F.J., Vabishchevich, P.N.: Staggered grid discretizations for the quasi-static Biot’s
consolidation problem. Appl. Numer. Math. 56(6), 888–898 (2006). https://doi.org/10.1016/j.apnum.
2005.07.002
12. Nordbotten, J.M.: Stable cell-centered finite volume discretization for Biot equations. SIAM J. Numer.
Anal. 54(2), 942–968 (2016). https://doi.org/10.1137/15M1014280
13. Murad, M.A., Loula, A.F.D.: Improved accuracy in finite element analysis of Biot’s consolidation
problem. Comput. Methods Appl. Mech. Eng. 95(3), 359–382 (1992). https://doi.org/10.1016/0045-
7825(92)90193-N
14. Murad, M.A., Loula, A.F.D.: On stability and convergence of finite element approximations of Biot’s
consolidation problem. Int. J. Numer. Methods Eng. 37(4), 645–667 (1994). https://doi.org/10.1002/
nme.1620370407
15. Murad, M.A., Thomée, V., Loula, A.F.D.: Asymptotic behavior of semidiscrete finite-element approxi-
mations of Biot’s consolidation problem. SIAM J. Numer. Anal. 33(3), 1065–1083 (1996). https://doi.
org/10.1137/0733052
16. Aguilar, G., Gaspar, F., Lisbona, F., Rodrigo, C.: Numerical stabilization of Biot’s consolidation model
by a perturbation on the flow equation. Int. J. Numer. Methods Eng. 75(11), 1282–1300 (2008). https://
doi.org/10.1002/nme.2295
17. Rodrigo, C., Gaspar, F., Hu, X., Zikatanov, L.: Stability and monotonicity for some discretizations of the
Biot’s consolidation model. Comput. Methods Appl. Mech. Eng. 298, 183–204 (2016). https://doi.org/
10.1016/j.cma.2015.09.019
18. Pé de la Riva, A., Gaspar, F., Rodrigo, C., Hu, X., Zikatanov, L.: Oscillation-free numerical schemes for
the Biot’s model with automatic stabilization. Submitted (2023)
19. Hu, X., Rodrigo, C., Gaspar, F.J., Zikatanov, L.T.: A nonconforming finite element method for the Biot’s
consolidation model in poroelasticity. J. Comput. Appl. Math. 310, 143–154 (2017). https://doi.org/10.
1016/j.cam.2016.06.003
20. Rodrigo, C., Hu, X., Ohm, P., Adler, J.H., Gaspar, F.J., Zikatanov, L.T.: New stabilized discretizations
for poroelasticity and the Stokes’ equations. Comput. Methods Appl. Mech. Eng. 341, 467–484 (2018).
https://doi.org/10.1016/j.cma.2018.07.003
21. Hong, Q., Kraus, J.: Parameter-robust stability of classical three-field formulation of Biot’s consolidation
model. Electron. Tran. Numer. Anal. 48, 202–226 (2018). https://doi.org/10.1553/etna_vol48s202
22. Oyarzúa, R., Ruiz-Baier, R.: Locking-free finite element methods for poroelasticity. SIAM J. Numer.
Anal. 54(5), 2951–2973 (2016). https://doi.org/10.1137/15M1050082
23. Lee, J.J.: Robust error analysis of coupled mixed methods for Biot’s consolidation model. J. Sci. Comput.
69(2), 610–632 (2016). https://doi.org/10.1007/s10915-016-0210-0
123
24. Boon, W.M., Kuchta, M., Mardal, K.-A., Ruiz-Baier, R.: Robust preconditioners for perturbed saddle-
point problems and conservative discretizations of biot’s equations utilizing total pressure. SIAM J. Sci.
Comput. 43(4), 961–983 (2021). https://doi.org/10.1137/20M1379708
25. Kim, J., Tchelepi, H.A.A., Juanes, R.: Stability, accuracy, and efficiency of sequential methods for coupled
flow and geomechanics. SPE J. 16(02), 249–262 (2011). https://doi.org/10.2118/119084-PA
26. Kim, J.: Sequential methods for coupled geomechanics and multiphase flow. PhD thesis, Stanford Uni-
versity (2010)
27. Mikelić, A., Wheeler, M.F.: Convergence of iterative coupling for coupled flow and geomechanics. Com-
put. Geosci. 17(3), 455–461 (2013). https://doi.org/10.1007/s10596-012-9318-y
28. Both, J.W., Borregales, M., Nordbotten, J.M., Kumar, K., Radu, F.A.: Robust fixed stress splitting for
Biot’s equations in heterogeneous media. Appl. Math. Lett. 68, 101–108 (2017). https://doi.org/10.1016/
j.aml.2016.12.019
29. Almani, T., Kumar, K., Dogru, A., Singh, G., Wheeler, M.F.: Convergence analysis of multirate fixed-
stress split iterative schemes for coupling flow with geomechanics. Comput. Methods Appl. Mech. Eng.
311, 180–207 (2016). https://doi.org/10.1016/j.cma.2016.07.036
30. Bause, M., Radu, F.A., Köcher, U.: Space-time finite element approximation of the Biot poroelasticity
system with iterative coupling. Comput. Methods Appl. Mech. Eng. 320, 745–768 (2017). https://doi.
org/10.1016/j.cma.2017.03.017
31. Borregales, M., Kumar, K., Radu, F.A., Rodrigo, C., Gaspar, F.J.: A partially parallel-in-time fixed-stress
splitting method for Biot’s consolidation model. Comput. Math. Appl. 77(6), 1466–1478 (2019). https://
doi.org/10.1016/j.camwa.2018.09.005. (7th International Conference on Advanced Computational
Methods in Engineering (ACOMEN 2017))
32. Castelletto, N., White, J.A., Tchelepi, H.A.: Accuracy and convergence properties of the fixed-stress
iterative solution of two-way coupled poromechanics. Int. J. Numer. Anal. Meth. Geomech. 39(14),
1593–1618 (2015). https://doi.org/10.1002/nag.2400
33. White, J.A., Castelletto, N., Tchelepi, H.A.: Block-partitioned solvers for coupled poromechanics: a
unified framework. Comput. Methods Appl. Mech. Eng. 303, 55–74 (2016). https://doi.org/10.1016/j.
cma.2016.01.008
34. Gaspar, F.J., Gracia, J.L., Lisbona, F.J., Oosterlee, C.W.: Distributive smoothers in multigrid for problems
with dominating grad-div operators. Numer. Linear Algebra Appl. 15(8), 661–683 (2008). https://doi.
org/10.1002/nla.587
35. Gaspar, F.J., Lisbona, F.J., Oosterlee, C.W., Wienands, R.: A systematic comparison of coupled and
distributive smoothing in multigrid for the poroelasticity system. Numer. Linear Algebra Appl. 11(2–3),
93–113 (2004). https://doi.org/10.1002/nla.372
36. Luo, P., Rodrigo, C., Gaspar, F.J., Oosterlee, C.W.: On an Uzawa smoother in multigrid for poroelasticity
equations. Numer. Linear Algebra Appl. 24(1), 2074 (2017). https://doi.org/10.1002/nla.2074.e2074nla.
2074
37. Gaspar, F.J., Rodrigo, C.: On the fixed-stress split scheme as smoother in multigrid methods for coupling
flow and geomechanics. Comput. Methods Appl. Mech. Eng. 326, 526–540 (2017). https://doi.org/10.
1016/j.cma.2017.08.025
38. Adler, J.H., He, Y., Hu, X., MacLachlan, S., Ohm, P.: Monolithic multigrid for a reduced-quadrature
discretization of poroelasticity. SIAM J. Sci. Comput. 20, 54–81 (2022). https://doi.org/10.1137/
21M1429072
39. Phoon, K.K., Toh, K.C., Chan, S.H., Lee, F.H.: An efficient diagonal preconditioner for finite element
solution of Biot’s consolidation equations. Int. J. Numer. Methods Eng. 55(4), 377–400 (2002). https://
doi.org/10.1002/nme.500
40. Bergamaschi, L., Ferronato, M., Gambolati, G.: Novel preconditioners for the iterative solution to FE-
discretized coupled consolidation equations. Comput. Methods Appl. Mech. Eng. 196(25), 2647–2656
(2007). https://doi.org/10.1016/j.cma.2007.01.013
41. Ferronato, M., Bergamaschi, L., Gambolati, G.: Performance and robustness of block constraint precon-
ditioners in finite element coupled consolidation problems. Int. J. Numer. Methods Eng. 81(3), 381–402
(2010). https://doi.org/10.1002/nme.2702
42. Haga, J.B., Osnes, H., Langtangen, H.P.: Efficient block preconditioners for the coupled equations of
pressure and deformation in highly discontinuous media. Int. J. Numer. Anal. Methods Geomech. 35(13),
1466–1482 (2011). https://doi.org/10.1002/nag.973
43. Castelletto, N., White, J.A., Ferronato, M.: Scalable algorithms for three-field mixed finite element coupled
poromechanics. J. Comput. Phys. 327, 894–918 (2016). https://doi.org/10.1016/j.jcp.2016.09.063
44. Bærland, T., Lee, J.J., Mardal, K.-A., Winther, R.: Weakly imposed symmetry and robust preconditioners
for Biot’s consolidation model. Comput. Methods Appl. Math. 17(3), 377–396 (2017). https://doi.org/
10.1515/cmam-2017-0016
123
45. Adler, J.H., Gaspar, F.J., Hu, X., Rodrigo, C., Zikatanov, L.T.: Robust block preconditioners for Biot’s
model. In: Bjørstad, P.E., Brenner, S.C., Halpern, L., Kim, H.H., Kornhuber, R., Rahman, T., Widlund,
O.B. (eds.) Domain Decomposition Methods in Science and Engineering XXIV, pp. 3–16. Springer, Cham
(2018). https://doi.org/10.1007/978-3-319-93873-8_1
46. Hong, Q., Kraus, J., Lymbery, M., Philo, F.: Conservative discretizations and parameter-robust precon-
ditioners for Biot and multiple-network flux-based poroelasticity models. Numer. Linear Algebra Appl.
26(4), 2242 (2019). https://doi.org/10.1002/nla.2242
47. Adler, J.H., Gaspar, F.J., Hu, X., Ohm, P., Rodrigo, C., Zikatanov, L.T.: Robust preconditioners for a
new stabilized discretization of the poroelastic equations. SIAM J. Sci. Comput. 42(3), 761–791 (2020).
https://doi.org/10.1137/19M1261250
48. Chen, S., Hong, Q., Xu, J., Yang, K.: Robust block preconditioners for poroelasticity. Comput. Methods
Appl. Mech. Eng. 369, 113229 (2020). https://doi.org/10.1016/j.cma.2020.113229
49. Elman, H., Howle, V.E., Shadid, J., Shuttleworth, R., Tuminaro, R.: Block preconditioners based on
approximate commutators. SIAM J. Sci. Comput. 27(5), 1651–1668 (2006). https://doi.org/10.1137/
040608817
50. Elman, H., Howle, V.E., Shadid, J., Shuttleworth, R., Tuminaro, R.: A taxonomy and comparison of
parallel block multi-level preconditioners for the incompressible Navier–Stokes equations. J. Comput.
Phys. 227(3), 1790–1808 (2008). https://doi.org/10.1016/j.jcp.2007.09.026
51. Brandt, A., McCormick, S.F., Ruge, J.W.: Algebraic multigrid (AMG) for sparse matrix equations. In:
Evans, D.J. (ed.) Sparsity and Its Applications. Cambridge University Press, Cambridge (1984)
52. Arnold, D.N., Falk, R.S., Winther, R.: Multigrid in H (div) and H (curl). Numer. Math. 85(2), 197–217
(2000). https://doi.org/10.1007/PL00005386
53. Hiptmair, R., Xu, J.: Nodal auxiliary space preconditioning in H( curl) and H(div) spaces. SIAM J. Numer.
Anal. 45(6), 2483–2509 (2007). https://doi.org/10.1137/060660588
54. Loghin, D., Wathen, A.J.: Analysis of preconditioners for saddle-point problems. SIAM J. Sci. Comput.
25(6), 2029–2049 (2004). https://doi.org/10.1137/S1064827502418203
55. Mardal, K.-A., Winther, R.: Preconditioning discretizations of systems of partial differential equations.
Numer. Linear Algebra Appl. 18(1), 1–40 (2011). https://doi.org/10.1002/nla.716
56. Klawonn, A., Starke, G.: Block triangular preconditioners for nonsymmetric saddle point problems:
field-of-values analysis. Numer. Math. 81(4), 577–594 (1999). https://doi.org/10.1007/s002110050405
57. Starke, G.: Field-Of-Values analysis of preconditioned iterative methods for nonsymmetric elliptic prob-
lems. Numer. Math. 78, 103–117 (1997). https://doi.org/10.1007/s002110050306
58. Hu, X., Adler, J.H., Zikatanov, L.T.: HAZmath: a simple finite element, graph, and solver library. https://
hazmathteam.github.io/hazmath/
59. Greenbaum, A.: Iterative Methods for Solving Linear Systems. Society for Industrial and Applied Math-
ematics, Philadelphia (1997). https://doi.org/10.1137/1.9781611970937
60. Elman, H.C.: Iterative methods for large, sparse, nonsymmetric systems of linear equations. PhD thesis,
Yale University New Haven, Conn (1982)
61. Eisenstat, S.C., Elman, H.C., Schultz, M.H.: Variational iterative methods for nonsymmetric systems of
linear equations. SIAM J. Numer. Anal. 20(2), 345–357 (1983). https://doi.org/10.1137/0720023
62. Arnold, D.N., Brezzi, F., Fortin, M.: A stable finite element for the Stokes equations. Calcolo 21(4),
337–344 (1985). https://doi.org/10.1007/BF02576171
63. Stenberg, R.: A technique for analysing finite element methods for viscous incompressible flow. Int. J.
Numer. Methods Fluids 11(6), 935–948 (1990). https://doi.org/10.1002/fld.1650110615. (The Seventh
International Conference on Finite Elements in Flow Problems (Huntsville, AL, 1989))
64. Ma, Y., Hu, K., Hu, X., Xu, J.: Robust preconditioners for incompressible MHD models. J. Comput. Phys.
316, 721–746 (2016). https://doi.org/10.1016/j.jcp.2016.04.019
65. Girault, V., Raviart, P.-A.: Finite Element Methods for Navier–Stokes Equations. Theory and Algorithms.
Springer Series in Computational Mathematics, vol. 5, p. 374. Springer, Berlin (1986). https://doi.org/
10.1007/978-3-642-61623-5
66. Kolev, T.V., Vassilevski, P.S.: Parallel auxillary space AMG solver for H(div) problems. SIAM J. Sci.
Comput. 34(6), 3079–3098 (2012). https://doi.org/10.1137/110859361
67. Abousleiman, Y., Cheng, A.H.-D., Cui, L., Detournay, E., Roegiers, J.-C.: Mandel’s problem revisited.
Géotechnique 46(2), 187–195 (1996). https://doi.org/10.1680/geot.1996.46.2.187
68. Davis, T.A.: Algorithm 832: UMFPACK, an unsymmetric-pattern multifrontal method. ACM Trans. Math.
Softw. 30(2), 196–199 (2004). https://doi.org/10.1145/992200.992206
69. Davis, T.A.: A column pre-ordering strategy for the unsymmetric-pattern multifrontal method. ACM
Trans. Math. Softw. 30(2), 165–195 (2004). https://doi.org/10.1145/992200.992205
70. Davis, T.A., Duff, I.S.: An unsymmetric-pattern multifrontal method for sparse LU factorization. SIAM
J. Matrix Anal. Appl. 18(1), 140–158 (1997). https://doi.org/10.1137/S089547989424690
123
71. Davis, T.A., Duff, I.S.: A combined unifrontal/multifrontal method for unsymmetric sparse matrices.
ACM Trans. Math. Softw. 25(1), 1–19 (1999). https://doi.org/10.1145/305658.287640
72. Raviart, P.-A., Thomas, J.M.: A mixed finite element method for 2nd order elliptic problems. In: Mathe-
matical Aspects of Finite Element Methods (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome,
1975). Springer, Berlin, pp. 292–315606 (1977). https://doi.org/10.1007/BFb0064470
73. Nédélec, J.-C.: A new family of mixed finite elements in R3 . Numer. Math. 50(1), 57–81 (1986). https://
doi.org/10.1007/BF01389668
74. Nédélec, J.-C.: Mixed finite elements in R3 . Numer. Math. 35(3), 315–341 (1980). https://doi.org/10.
1007/BF01396415
75. Crouzeix, M., Raviart, P.-A.: Conforming and nonconforming finite element methods for solving the
stationary Stokes equations I. Rev. Française Autom. Inf. Rech. Opér. Sér. Rouge 7(R–3), 33–75 (1973).
https://doi.org/10.1051/m2an/197307R300331
76. Falk, R.S.: Nonconforming finite element methods for the equations of linear elasticity. Math. Comput.
57(196), 529–550 (1991). https://doi.org/10.2307/2938702
77. Falk, R.S., Morley, M.E.: Equivalence of finite element methods for problems in elasticity. SIAM J.
Numer. Anal. 27(6), 1486–1505 (1990). https://doi.org/10.1137/0727086
78. Hansbo, P., Larson, M.G.: Discontinuous Galerkin and the Crouzeix–Raviart element: application to
elasticity. Math. Model. Numer. Anal. 37(1), 63–72 (2003). https://doi.org/10.1051/m2an:2003020
Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and
institutional affiliations.
123

Parameter-Robust Preconditioners For Biot's Model

Uploaded by

Copyright:

Available Formats

Parameter-Robust Preconditioners For Biot's Model

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Parameter-Robust Preconditioners For Biot's Model

Uploaded by

Copyright:

Available Formats

SeMA Journal (2024) 81:51–80

Parameter-robust preconditioners for Biot’s model

Carmen Rodrigo1 · Francisco J. Gaspar1 · James Adler2 · Xiaozhe Hu2 ·

Keywords Biot’s model · Poroelasticity · Well-posedness · Preconditioning ·

Mathematics Subject Classification 65F08 · 65F10 · 65N12 · 65N22 · 65N30

2.1 Norm-equivalent preconditioner

Let M : X  → X be a symmetric positive definite (SPD) operator, which induces an

c1 x2X ≤ x2M−1 ≤ c2 x2X , (6)

2.2 FOV-equivalent preconditioner

In this section, we study a more general class of preconditioners as the field-of-values-

In general, M : X → X can be any SPD operator, but we will consider M to be a SPD

3 Parameter-robust preconditioners for Biot’s model

3.1 Two-field formulation

V b = span{ϕb,T e1 , . . . , ϕb,T ed }T ∈Th , ϕb,T = αT λ1,T . . . λd+1,T ,

a(unh , v h ) − α(div v h , phn ) = (g(tn ), v h ), ∀v h ∈ V h , (13)

where ∂¯t unh := (unh − un−1 ¯ n n−1

where a(u, v) → A u , −(div u, q) → B, a p ( p, q) → A p , and (∇ p, ∇q) → L p represent

With this purpose, the key is to define a proper weighted norm on X = V h × Q h as

L(x, y) = (A u u, v) + α(Bv, p) + α(Bu, q) − τ (μ−1

L(x, y) = (A u u, u+θ w)+α(B(u + θ w), p)−α(Bu, p)+τ (K ∇ p, ∇ p)+ηh 2 (∇ p, ∇ p)

where M is the mass matrix of the pressure block.

c1,u (Su u, u) ≤ (A−1

being the involved constants independent of discretization and physical parameters.

First, we consider the following block lower triangular preconditioner:

and we need to show that the FOV-equivalence in (8) is satisfied.

Theorem 2 There exist constants and ϒ, independent of discretization or physical param-

provided that η = δ αζ 2 with δ > 0.

The inexact counterpart of (24) is also considered,

Theorem 3 If I −Su A u  A u ≤ ρ with 0 ≤ ρ < 1, there exist constants and ϒ, independent

provided that η = δ αζ 2 with δ > 0.

Table 1 Two-field stabilized BD BL BU

ν = 0.2 and varying K K = 10−6 and varying ν

3.2 Three-field formulation

a(u, v) − α( p, div v) = (ρ g, v), ∀ v ∈ V , (30)

approximation to (u(·, tm ), p(·, tm ), w(·, tm )) , at time tm = mτ, m = 1, 2, . . . . Note that

3.2.1 Stabilized P1 −RT−P0

More specifically, to discretize displacements we choose a piecewise linear continuous

Theorem 4 If the triple (V h , W h , Q h ) is Stokes–Biot stable, that is,

where Aw := BwT M −1p Bw .

c1,u (Su u, u) ≤ (A−1 u u, u) ≤ c2,u (Su u, u), (44)

(B L A x, x)(B D )−1 B L A x(B D )−1

Analogously to the two-field formulation case, we also consider an inexact approach to

Fig. 2 2D physical and

and the corresponding inexact preconditioner,

where p0 = B(1 + νu )F, B = 1 is the Skempton’s coefficient, νu = 3ν+B(1−2ν)

ν = 0.0 and varying K K = 10−6 and varying ν

3.2.2 Non conforming CR−RT−P0

as well as its inexact counterpart given by

Numerical experiment: poroelastic problem on a square domain with a uniform load

p = 0, σ · n = (0, −1)t , on 1 = [0, 1] × 1.

Fig. 3 Computational domain

Table 5 Three-field CR−RT−P0 BD BL BU

ν = 0.2 and varying K K = 10−6 and varying ν

Data availability Not applicable.

You might also like

Let M : X → X be a symmetric positive definite (SPD) operator, which induces an

c1 x2X ≤ x2M−1 ≤ c2 x2X , (6)

Theorem 3 If I −Su A u A u ≤ ρ with 0 ≤ ρ < 1, there exist constants and ϒ, independent

(B L A x, x)(B D )−1 B L A x(B D )−1