Monolithic two-level Schwarz preconditioner for Biot's consolidation model in two space dimensions
Stefan Meggendorfer Guido Kanschat Johannes Kraus

Reusing Deep Learning Models: Challenges and Directions in Software Engineering
    
James C. Davis1, Purvish Jajal1, Wenxin Jiang1,

Taylor R. Schorlemmer1, Nicholas Synovic2, and George K. Thiruvathukal2 

1Department of Electrical & Computer Engineering, Purdue University, West Lafayette, IN, USA

{davisjam,pjajal,jiang784,tschorle}@purdue.edu

2Department of Computer Science, Loyola University Chicago, Chicago, IL, USA

{nsynovic,gthiruvathukal}@luc.edu
All authors are listed in alphabetical order as equal contributors to this manuscript. Please direct correspondence to davisjam@purdue.edu and gthiruvathukal@luc.edu.

    April 29, 2024
=====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




This paper addresses the construction and analysis of a class of domain decomposition methods
for the iterative solution of the quasi-static Biot problem in three-field formulation. 
The considered discrete model arises from time discretization by the implicit Euler method and space discretization by a family of strongly mass-conserving methods exploiting -conforming approximations of the solid displacement and
fluid flux fields.
For the resulting saddle-point problem, we construct monolithic overlapping domain decomposition (DD) methods whose analysis relies on a transformation into an equivalent symmetric positive definite system and on stable decompositions of the involved finite element spaces under proper problem-dependent norms.
Numerical results on two-dimensional test problems are in accordance with the provided theoretical uniform convergence estimates for the two-level multiplicative Schwarz method.


MSC-Classes: 65M12, 65M55, 65F10


Keywords: quasi-static Biot model, poroelasticity, domain decomposition, preconditioning, mass-conservative -conforming discontinuous Galerkin






§ INTRODUCTION




Mathematical models describing the interaction of fluid flow with the mechanical deformation of porous media
go back to the pioneering works by Karl von Terzaghi and Maurice Anthony Biot in the first half of the last
century <cit.>.
Important application areas include, but are not limited to, petroleum and reservoir engineering, CO2 sequestration
and the mechanics of biological tissues. Partial saturation, the incorporation of two- or multi-phase flow, thermal
processes, and/or large deformations, and/or the consideration of multiple fluid networks lead to even more complex,
sometimes highly nonlinear mathematical models that require physically correct (structure-preserving) discretization
and robust iterative solution methods. 

The well-posedness of the basic Biot model of consolidation on the continuous and discrete levels has
been established in the works of Zenisek and Showalter using semi-group theory and Galerkin discretization
methods <cit.>.

The simplest linear model uses a two-field formulation of Biot's equations for the solid displacement   and pore
pressure p. Popular time discretizations of the quasi-static differential algebraic system are the backward Euler
or high-order strongly stable time integration methods, see, e.g., <cit.>.
A stability and convergence analysis of inf-sup stable finite element approximations of this two-field model has
first been presented in <cit.> in the setting of implicit Euler time stepping.

Discretizations based on three-field-formulation have originally been proposed
in <cit.> where standard continuous Galerkin
approximations of the displacement variable  are coupled to a mixed finite-element method for the pore pressure p
based on introducing the flux variable  to the system, which obeys the Darcy law. This approach has also been extended
to discontinuous Galerkin approximations of  in <cit.> and other nonconforming approximations,
e.g., using modified rotated bilinear elements <cit.>, or Crouzeix-Raviart elements in <cit.>.
Lately, the conservation of mass was identified essential for parameter robustness.
In <cit.>, families of strongly mass conserving discretizations based on the
-conforming discontinuous Galerkin discretization of  have been suggested and proven stable independently
of model and discretization parameters (Lamé parameters, permeability and Biot-Willis parameter, storage coefficient, time step
and mesh size), see, e.g. <cit.>, for time-dependent error estimates, see <cit.>.
More recently, mass-conservation was achieved by finite elements based on enrichment in <cit.>
and by mixed virtual finite elements in <cit.>.

A discretization of a four-field formulation of Biot's model has been proposed in <cit.>. The method uses the
(effective) stress tensor σ, the displacement , the fluid flux , and the pore pressure p as unknowns and
couples two standard mixed finite element methods, one for the mechanics and one for the flow subproblem, the former of which
is based on the Hellinger-Reissner formulation. Optimal a priori error estimates were shown for both semidiscrete and fully discrete
problems when the flux unknown is approximated in the Raviart-Thomas space and the stress tensor in the Arnold-Winther space.
The error analysis of this coupled mixed method has been complemented in <cit.> by estimates in L^∞-norm
in time and L^2-norm in space that are robust with respect to the Lamé parameters and do not require strict positivity of the
constrained storage coefficient .

When it comes to the solution of the algebraic problems arising from discretization of Biot's problem, preconditioning techniques
have been developed for various formulations and discretizations, in particular, in the framework of norm- or field-of-values-equivalent
operator preconditioners, see,
e.g. <cit.> for two-field,
<cit.> for three-field, and <cit.>
for four-field formulations. 
Here, <cit.> focus on multiple network extensions of Biot's model and <cit.> on a generalization to Biot-Brinkman equations. 
The key task in the development of such methods is the
stability analysis of the underlying saddle-point problems, which can be performed in the framework presented in <cit.>.

By a decomposition into physical subsystems, the coupled static problem can be solved implicitly, using a loose or explicit coupling,
or an iterative coupling. The latter, often provides an attractive alternative in terms of achieving high accuracy at reasonable computational
cost. The most popular procedures in this category are the undrained split, the fixed-stress split, the drained split and the fixed-strain split
iterative methods, which also generate preconditioners for Krylov methods such as the GMRES method. As shown in <cit.>,
in contrast to the drained split and the fixed-strain split methods, the undrained split and fixed-stress split methods are unconditionally stable. 
The first convergence analysis of the latter methods has been presented in <cit.> for the quasi-static Biot system.
More recent works on fixed-stress type split methods have focused on the parameter-robustness of the scheme and the optimization of the stabilization/acceleration parameter <cit.>, extensions to heterogeneous poroelastic media <cit.>, and generalizations to nonlinear poroelasticity models <cit.>. 
For a parameter-robust convergence analysis of splitting methods that extends also to multiple network poroelastic theory (MPET) equations, see <cit.>.

However, no matter whether using classical block diagonal or block triangular preconditioners or splitting methods, the efficient solution of the subsystems for the individual physical fields remains an important issue. Domain decomposition and/or (algebraic) multigrid methods provide a feasible approach to solve the arising H^1 and  subproblems although their parameter-robust performance, e.g, for nearly incompressible materials and/or low hydraulic conductivities can be quite challenging.
In any case, this approach requires to combine at least two (nearly) optimal solvers for subproblems, such as elasticity and Darcy flow, 
in order to obtain an overall (nearly) optimal method.

A monolithic multigrid approach for a reduced-quadrature discretization has recently been developed in <cit.> where an additive Vanka relaxation method, based on overlapping vertex-patches, and an inexact Braess-Sarazin relaxation method is considered.
The methods rely on optimized relaxation factors depending on the Poisson ratio.
Their numerical tests confirm robustness of the method in the permeability coefficient and for moderately incompressible materials.
In comparison, our method is robust with respect to the Poisson ratio.

In this work we design a monolithic two-level Schwarz method based on overlapping vertex-patches for the quasi-static Biot problem in three-field formulation similar to <cit.>. However, instead of a reduced-quadrature discretization we apply the preconditioner to a family of strongly mass-conserving mixed finite element methods using an -conforming discontinuous Galerkin ansatz for the displacement field, as proposed in <cit.>.
This allows us to prove convergence of the method by transforming the saddle-point problem into an equivalent singularly perturbed symmetric
positive definite problem, cf. <cit.>.
In particular, we obtain a robust method with standard choices of the relaxation parameters under the assumption of two-sided bounds on the permeability tensor.

The remainder of the article is organized as follows:
we introduce Biot's consolidation model and its conservative discretization in <ref>. The monolithic two-level algorithm is presented
in <ref> and analyzed in <ref>. In <ref> we provide numerical evidence of our theoretical findings.





§ MODEL AND DISCRETIZATION




Biot's consolidation model in three-field formulation couples the displacement field  of the solid component, the seepage velocity  of the fluid, and the fluid pressure p in a fully saturated porous medium at constant temperature.
It is posed on a computational domain Ω×(0,T) consisting of a spatial domain Ω⊂^d, d∈{1,2,3}, where in this work we focus on the case d=2, and a time interval (0,T), T∈. In strong form, Biot's quasi-static model of consolidation reads

    1pt
        
            -÷(2μ
            +
            λ÷ ()  I )
                     +   ∇ p
               =     Ω× (0,T),
            
             +   ∇ p
               =    0   Ω× (0,T),
            
    
            -÷∂_t    -   ÷   -   ∂_t p
               =    g   Ω× (0,T),

with initial conditions (0)=_0, and p(0)=p_0 at time t=0.
The set of equations describes the momentum balance, Darcy's law, and the mass balance, respectively.
Here, the physical parameters λ>0 and μ>0 denote the Lamé coefficients of linear elasticity, >0 denotes the Biot-Willis constant,  is the symmetric and positive definite permeability coefficient of the porous medium, and ≥0 is the specific storage coefficient.
Further, the strain tensor is defined by := (∇+∇^T)/2.
For an introduction to the system of equations, as well as a discussion of the physical parameters we refer the reader to <cit.>.

We apply a semi-discretization in time by the backward Euler method with time step size τ>0, resulting in a
semi-discrete static problem of the form

    1pt
        
            - 2 μ÷ 
            -
            λ∇÷         +   ∇ p
               =     Ω,
            
             +   ∇ p
               =    0   Ω,
            
    
            -/τ÷   -   ÷   -   /τp
               =   g̅  Ω,

in every time step where g̅ := g -/τ÷_ old - /τp_ old with _ old and
p_ old taken from the previous time moment and all other quantities referring to the new time moment. 

We point out that the specific time discretization is irrelevant for this article, since all implicit methods will result in the necessity
of solving systems of similar kind. Furthermore, we note that /τ has the unit of speed and hence it is meaningful to add
the quantities (÷)/τ, ÷, and p/τ in the last equation.

Following <cit.>, by proper scaling and substitution of variables, we can transform (<ref>)
into a symmetric problem.
First, from (<ref>) we have

    1pt
        
            -÷
            -
            λ/2μ∇÷         +   /2μ∇ p
               =   1/2μ     Ω,
            
          /2μ^-1   +   /2μ∇ p
               =    0      Ω ,
            
    
            - ÷   -   τ/÷   -   /p
               =   τ/g̅     Ω.

Next, by substituting

    = τ/,
          p̃ = /2μ p,
           = 1/2μ,
          g̃ = τ/g̅ - ÷_ old - /p_ old,
        
    ^-1 = ^2/2μτ^-1,
           = λ/2μ,
           = 2μ/^2,

we obtain the rescaled system

    1pt
        
            -÷
            +
            ∇(÷)
                     +   ∇p̃   =        Ω,
            
          ^-1   +   ∇p̃   =    0      Ω,
            
    
            -÷   -   ÷   -   p̃   =   g̃     Ω.

Since this is the normalized system we build our analysis on, we will omit the tilde symbol from here on, but we will keep the hat to clarify the scaling of the parameters.

Throughout this paper, we assume that the scaled permeability tensor  is bounded such that

    ^2 ≤(^-1,) ≤^2
          ∀∈ [L^2(Ω)]^2,

for two positive constants  and .
Here, ·=·_L^2(Ω) and (·,·)=(·,·)_L^2(Ω) denote the L^2-norm
and L^2-inner product on Ω, respectively.
Note that condition (<ref>) follows, e.g., from the pointwise condition

    ≤_min(^-1) ≤_max(^-1) ≤,
          a. e. in Ω,

where _min(^-1) and _max(^-1) denote (pointwise) the smallest and largest eigenvalues of ^-1, respectively.
As a matter of fact, the constants in the analysis below depend on  and , so that it neither intends to cover high contrast coefficients nor the case of strong anisotropy.
In the simplest case, we consider a permeability tensor of the form = I with a scalar , constant over Ω and constant over the time interval (0,T).

This system is closed by proper boundary conditions. As usual, essential boundary conditions enter the definitions
of the function spaces, here U⊂ [H^1(Ω)]^2 and V⊂(Ω) for displacement and seepage velocity,
respectively.
We will choose a particular combination of boundary conditions in the presentation and analysis of the two-level algorithm
in Sections <ref> and <ref>, although one might also want to consider other combinations
for which certain arguments/basic estimates have to be adapted/exchanged accordingly.
The general assumptions that are required to conduct an analysis similar to the one presented in this paper are as follows.



  The boundary conditions on the spaces U and V are such that:
  
      
  * The seepage velocity  is uniquely determined.
      
  * Korn's inequality holds for  in the continuous setting as well as for _h in the discrete setting. For sufficient conditions
      in the continuous case we refer to <cit.> and the literature cited there. For the discrete case, see <cit.>.
      
  * The pressure space Q is chosen such that the operator ÷ V→ Q is surjective. The boundary conditions on U are compatible such that also ÷ U→ Q is surjective.
  


To give one example, which applies to the setting of the analysis presented in Section <ref>, Korn's inequality

    1/c_e∇^2 ≤^2 ≤∇^2

holds for all ∈ [H^1_0(Ω)]^2.

In line with this, on the whole boundary we assume a homogeneous Dirichlet (no slip) boundary condition for  and a Neumann boundary condition for p, which translates to an essential condition for ·.
Accordingly, we choose

    U=[H^1_0(Ω)]^2,      V=_0(Ω),      Q=L^2_0(Ω),

where


    _0(Ω)    = ∈(Ω) : ·=0  on ∂Ω,

in the sense of traces of  and

    L^2_0(Ω)    = q∈ L^2(Ω):∫_Ωq dx=0.


The weak formulation of (<ref>) is then:
find (,,p)∈ U× V× Q, such that

    1pt
        [ (,)
            +
            (,)                         -                      (p,)                         =                       (,)                     ∀∈ U,;                                                 (,)                         -                      (p,)                         =                         0                     ∀∈ V,;                     -(q,)                         -                      (q,)                         -                     (p,q)                         =                     (g,q)                   ∀ q∈ Q, ]

with bilinear forms

    (,)    = (,),    (,)    = (÷,÷), 
    (p,)    = (p,÷),    (,)    = (^-1,).






 §.§ Discretization


The model problem (<ref>) is discretized by a family of mixed strongly mass-conserving methods
as proposed in <cit.>, where the displacement and seepage velocity
fields are approximated in suitable -conforming spaces U_h and V_h and the pressure space  Q_h
consists of piecewise polynomial functions discontinuous at element interfaces such that the conditions

    ÷ U_h = Q_h,     ÷ V_h = Q_h

are satisfied.
The combined finite element space for the mixed method is then defined by X_h := U_h× V_h× Q_h.
Since the space U_h is not H^1-conforming, a discrete interior penalty discontinuous Galerkin bilinear form
(·,·) is used to approximate (·,·) as detailed below.

Let  be a family of shape regular triangulations of the computational domain Ω into mesh cells 
with diameter h_=() where h:=max_∈h_ denotes the mesh size.
Let  be the set of all interior faces (edges in two dimensions) of , and  be the set
of all faces on the boundary ∂Ω.
For every ∈ there are two neighboring cells _+,_-∈ such that
=∂_+∩∂_-. Let  be the unit outward normal vector pointing from _+ to _-.
Then on every face ∈ and for any ∈ [L^2(Ω)]^2 and ∈ [L^2(Ω)]^2 × 2,
we define jump  and average
 by

    = _+ - _-,      = 1/2(_+ + _- ) ,

where _±=|__± and _±=|__±.
Further on the broken Sobolev space

    [H^1(Ω,)]^2 = {∈ [L^2(Ω)]^2  | |_∈ [H^1()]^2}

we introduce the discrete norm

    = (
        ∑_∈[]∇^2
        + ∑_∈1/h[]^2
        + ∑_∈1/h[]^2
        )^1/2.

On U_h× U_h we define the discrete bilinear form

    (_h,)
        =   ∑_∈[](_h,)
        + ∑_∈η/h[](_h,) 
       - ∑_∈[](_h,)
        - ∑_∈[](_h,) 
       + ∑_∈η/h[](_h,)
        - ∑_∈[](_h,)
        - ∑_∈[](_h,).

Here, the penalty parameter η>0 is chosen large enough to ensure coercivity of (·,·), i.e.,
there is a positive constant c such that

    (,) ≥ c ^2   ∀∈ U_h.

In addition we have continuity of (·,·) in the norm ·, i.e.,

    (_h,) ≤ c _h  ∀_h,∈ U_h.

We note that since U_h⊂(Ω) 
the jumps in face terms in (<ref>)
are zero in normal direction and hence can equivalently be defined using tangential components only. 

The mass-conserving mixed method based on the finite element space X_h can then be represented as

    (_h_hp_h,q)
        = (q),

where the discrete bilinear form (·,·) X_h× X_h→ is defined by

    (_h_h,q)
        =   (_h,) + (_h,) + (_h,) 
        - (p_h,+)
        - (q,_h+_h) - (p_h,q)

and the right hand side is given by

    (q) = (,) + (g_h,q),

with  chosen as the L^2-projection g.
System (<ref>) is consistent and uniformly well-posed, which has been proven in <cit.>
and follows essentially by the special choice of -conforming discretization spaces, coercivity (<ref>),
continuity (<ref>) and the discrete inf-sup conditions,
see, e.g. <cit.>)

    inf_q∈ Q_hsup_∈ U_h(÷,q)/q≥ > 0

and

    inf_q∈ Q_hsup_∈ V_h(÷,q)/q≥ > 0,

where

    = ( ^2 + ÷^2 )^1/2.






§ TWO-LEVEL ALGORITHM





In this section, we define a monolithic two-level algorithm for the mixed problem (<ref>).
The method differs from partitioned solvers which rely on and exploit the block structure of the system.

Let {Ω_j}_j=1^J be an overlapping covering of the domain Ω into subdomains Ω_j such that the following two standard assumptions are met, see <cit.>.

[Finite covering]

The partition {Ω_j}_j=1^J can be colored using at most N^c colors, in such a way that subregions with the same color are disjoint.


[Sufficient overlap]

There is a constant δ>0 that measures the size of overlaps between the subdomains Ω_j, for j=1,...,J, such that

    c_1 h≤δ≤ c_2 H

with constants c_1>0 and c_2>0.


In our experiments, we focus on the vertex patch smoother, where the subdomains Ω_j can consist of all mesh cells sharing the vertex j, i.e., the vertex patch associated with vertex j.
In the special case of a uniform partition of the domain into rectangles this would be a union of 4 mesh cells for interior vertices.

With the subdomains Ω_j we associate subspaces X_j⊂ X_h.
To this end, we define

    U_j := U_h∩_0(Ω_j),  
            V_j := V_h∩_0(Ω_j),  
            Q_j := Q_h∩ L^2_0(Ω_j),

continued by zero on Ω∖Ω_j, and set

    X_j := U_j× V_j× Q_j.

Moreover, for H>h let  be a triangulation of the coarse domain Ω_0=Ω such
that the associated subspace X_0:=X_H ⊂ X_h. Then, X_0 plays the role of a global coarse space and
we can decompose X_h into an overlapping sum by

    X_h=X_0 + ∑_j=1^J X_j.

For j=1,…,J we introduce local bilinear forms (·,·) X_j× X_j→ as restrictions
of (·,·) to X_j and define projections _j X_h→ X_j by

    ( _j_h_h,_j_jq_j) = (_h_h,_j_jq_j)
            ∀_j_jq_j∈ X_j.

Additionally, on the coarse space X_0=X_H, we introduce a coarse bilinear form [0](·,·) X_0× X_0→ as [0](·,·)=(·,·) and define _0 X_h→ X_0 by

    [0]( _0_h_h,_H_Hq_H) = (_h_h,_H_Hq_H)
            ∀_H_Hq_H∈ X_0.

Note, that the bilinear form (·,·) is not inherited from (·,·) although X_H⊂ X_h, because it differs in the face and boundary terms.
Thus, _0 is not a projection, since consequently [0](·,·) differs from (·,·).

The bilinear form (·,·) X_h× X_h→ given in (<ref>) defines
an operator := _h: X_h → X'_h via the relation

    ⟨_h, _h ⟩_X'_h × X_h = (_h, _h)     ∀_h, _h ∈ X_h,

where X'_h denotes the dual space of X_h and ⟨·, ·⟩_X'_h × X_h the corresponding duality pairing.
With this notation we define the two-level additive Schwarz preconditioner : X'_h → X_h by

    = ∑_j=0^J _j ^-1,

where  is a relaxation factor depending on the overlap that is also used for tuning the method.
Further, the multiplicative two-level Schwarz preconditioner  is defined by

    = (I - ) ^-1,
           = (I - _J)⋯(I - _1)(I - _0).

Note, that in practice only local problems are solved and the global inverse ^-1 is actually never computed in the application of the preconditioners.






§ TWO-LEVEL CONVERGENCE ANALYSIS





We begin this section by stating the main results:



Assume that H ≤ c h for some positive constant c and that the bounds (<ref>) on  hold.
Then, the additive Schwarz method (<ref>) defines a uniform preconditioner for the system (<ref>).
Consequently, the GMRES method converges uniformly, i.e.,
it achieves a prescribed accuracy within a number of iterations
that is bounded independently of the mesh size h and any other model parameters than , when applied to solve a linear system with .




   Under the assumptions of <ref> the stationary iterative method based on the multiplicative Schwarz  
   preconditioner (<ref>) converges uniformly, i.e., with a residual contraction factor strictly less than 1, independent of the mesh size h and any other model parameters than .



The dependence of the presented two-level convergence estimates on the parameter  is caused by the same dependence of the constant in the stable decomposition of the divergence-free subspace V_h^0 of V_h under the norm induced by the bilinear form (·,·), see (<ref>) in the proof of Lemma <ref>.
We will use a very crude upper bound at this point, which is given by the global contrast / for reasons of simplicity because this aspect is not the focus of the presented analysis.
Improved bounds for Darcy have been developed based on weighted Poincaré inequalities under the assumption of certain quasi-monotonicity properties of the coefficient function ^-1 in  (·,·), see <cit.>. 


Due to the lack of an abstract convergence theory for subspace correction methods for indefinite problems, 
we follow <cit.> and transform the saddle point problem into an equivalent, singularly perturbed 
symmetric positive definite (SPD) system, which is possible for any fixed storage coefficient > 0, as 
detailed in Section <ref>.

The analysis of domain decomposition methods for this singularly perturbed problem requires a stable decomposition derived in subsections <ref> and <ref>.
Thus, the convergence with respect to the singularly perturbed problem follows by the theory for symmetric positive definite problems, which is proved in <ref>.
Finally, the proofs of <ref> can be found in subsection <ref>. They rely on the fact that the estimates for the singularly perturbed problem are uniform for → 0.




 §.§ An equivalent, singularly perturbed problem


As shown in <cit.>, strong mass conservation can be recovered from system (<ref>)
by using -conforming discretizations for the displacement field and seepage velocity in combination with piecewise
polynomial approximations of the pressure if the polynomial degrees match. To be more precise, the discrete mass balance
equation

    -÷_h - ÷_h - =

is fulfilled point-wise if the condition (<ref>) is met.
This allows us to substitute

    = - ^-1(÷_h + ÷_h + )

and obtain the equivalent, singularly perturbed system

    (_h,) + (_h,) + ^-1(_h+_h,)    = (,) - ^-1(,÷)
    (_h,) +  ^-1(_h+_h,)    = - ^-1(,÷)

for all (,)∈ = U_h× V_h.
We rewrite this system as

    (_h_h,)
        = ()

with

    (_h_h,)
        = (_h_h,)
        + (_h_h,),

where the bilinear forms (·,·) and (·,·)
are defined on × by

    (_h_h,)
           = (_h,) + (_h,), 
        
    (_h_h,)
           = (_h,) + ^-1(_h+_h,+)

and the right hand side is given by

    () = (,) - ^-1(,÷( + )).



  
  The form (·,·) is symmetric and coercive with respect to the weighted norm
  
    _h = (_h^2 + ÷_h^2 + ^-12_h^2 + ^-1÷ (_h+_h)^2 )^1/2

  for all _h=(_h,_h)∈, i.e.,
  
    (_h_h,_h_h)
        ≥ c_h^2.



Similar to <ref> we choose local subspaces ⊂ associated with the subdomains Ω_j,
such that

    = U_j× V_j

as well as a global coarse space W_0=U_H× V_H. Hence, every ∈ admits a decomposition of the form

    = ∑_j=0^J _j,
             w_j∈.


For j=1,…,J we introduce the local bilinear forms (·,·)×→ as restrictions of
(·,·) to  and a coarse bilinear form [0](·,·) W_0× W_0→
as [0](·,·)=(·,·). Analogously to (<ref>) and (<ref>) we define
the projections _j W_h→ W_j by

    ( _j_h_h,_j_j) = (_h_h,_j_j)
            ∀_j_j∈ W_j

and the projection-like operator _0 W_h→ W_0 by

    [0]( _0_h_h,_H_H) = (_h_h,_H_H)
            ∀_H_H∈ W_0.


Furthermore, the bilinear form (·,·)×→ given in (<ref>) defines
an operator :=_h: →['] via the relation

    ⟨_h, _h ⟩_['] × = (_h, _h),   ∀_h, _h ∈.


Thus, for the singularly perturbed problem the two-level additive Schwarz preconditioner : ['] → is defined by

    = ∑_j=0^J _j ^-1,

and the multiplicative two-level Schwarz preconditioner  is defined by

    = (I - ) ^-1,
           = (I-_J)⋯(I-_1)(I-_0).







 §.§ Decomposition of the spaces according to ker(D)



Our goal is to prove that for any ∈ W_h any decomposition of the form (<ref>)  
is stable under the global and local bilinear forms (·,·)  W_h × W_h→
and (·,·)×→, respectively.
As a preparatory, step we decompose the space  into the orthogonal sum

    = ⊕^⊥

of the kernel of the summed divergence operator

    W_h^0 :=  = {∈ W_h| (,) = 0   ∀∈ W_h}

and its (·,·)-orthogonal complement

    ^⊥ = {∈ W_h | (,)=0  ∀∈},

cf. <cit.>.
To this end, we introduce the discrete Helmholtz-decompo­sitions

    U_h = U_h^0 ⊕ U_h^⊥,
            
        V_h = V_h^0 ⊕ V_h^⊥,

of U_h and V_h into their divergence free subspaces U_h^0 and V_h^0
and their (·,·)-ortho­gonal and (·,·)-orthogonal
complements, respectively,

    U_h^⊥   = {∈ U_h | (,)=0  ∀∈ U_h^0},
    
        V_h^⊥   = {∈ V_h | (,)=0  ∀∈ V_h^0},

cf. <cit.>.
Accordingly, we denote the Helmholtz-decompositions of the local and coarse spaces U_j and V_j into their divergence-free subspaces and their (·,·)-orthogo­nal and (·,·)-orthogonal complements by

    U_j = U_j^0 ⊕ U_j^⊥,
            
        V_j = V_j^0 ⊕ V_j^⊥.


The following lemma characterizes W_h^0 = ⊂ W_h.


For >0 there holds

    [0] :=  = U_h^0 × V_h^0,
            [⊥] := ^⊥ = U_h^⊥× V_h^⊥.




Setting = and = in
  
    (,) = (,)
        + ^-1(+,+)

shows that a necessary condition for (,)^T ∈ is ÷ = 0, i.e., ∈ U_h^0.
Further, for any ∈ U_h^0, choosing = and =- shows that in this case ∈ V_h^0
is necessary and sufficient to guarantee that (,)^T ∈[0].


Note, that for all _h^0 = (_h^0,_h^0) ∈[0] and _h = (_h,_h)∈, there holds

    (_h^0, _h) =
      (_h^0,_h) + (_h^0,_h).





 §.§ Stable decompositions



The main result of this section is <ref>. Before we can prove it, we begin by recalling known facts from the literature
and showing some auxiliary results.
We start with the existence and some standard properties of a partition of unity for the overlapping covering {Ω_j}_j=1^J in <ref>, see <cit.>.



Given <ref>, there is a piecewise bilinear
partition of unity {θ_j}_j=1^J relative to the overlapping partition {Ω_j}_j=1^J, i.e.,
supp(θ_j) ⊂Ω_j, such that for all x∈Ω there holds

    ∑_j=1^J θ_j(x) = 1,    
        0 ≤θ_j(x) ≤ 1,     j=1,2,…,J.

Moreover,

    [(Ω)]∇θ_j≤c/δ,

where δ is the constant from <ref>.


In order to represent divergence free functions in the displacement
space U_h^0, we introduce the space

    S(Ω,) = { s∈ H^1_0(Ω) |  s_|∈ H^2()},

which is equipped with the discrete norm

    s = (
        ∑_∈[]∇^2 s^2
        + ∑_∈1/h[]∂_n s^2
        + ∑_∈1/h[]∂_n s^2
        )^1/2.

Following <cit.>, the bilinear form
(·,·) on the divergence free subspace is
algebraically equivalent to the bilinear form of a C^0-interior
penalty formulation of a corresponding biharmonic problem by assigning
velocities _h =  s, where _h belongs to the discrete stream function space S(Ω,).
Combining arguments from <cit.> and <cit.>, it is easy to see that in two space dimensions there holds the elementary identity

    s
      = s,   ∀ s∈ S(Ω,),

where = (∂_2,-∂_1).

Next, we recall a result showing that stream functions in S(Ω,) have a stable decomposition with respect to the norm (<ref>).

 
Every s∈ S(Ω,) admits a decomposition s=∑_j=0^J s_j with s_j∈ S(Ω_j,), s_0∈ S(Ω_0,), such that

    [H]s_0^2 + ∑_j=1^J s_j^2 ≤ c (1+H^4/δ^4) s^2

for some constant c>0.



  Compare the proof of Lemma 4.2 in <cit.>.


The above results suffice to prove the stability of the decomposition of the divergence free subspace of U_h with respect to the elasticity energy norm.


Every function ^0∈ U_h^0 admits a decomposition of the form ^0=∑_j=0^J ^0_j with ^0_j∈ U_j^0, which satisfies the bound

    (^0_0,^0_0) + ∑_j=1^J (^0_j,^0_j)
        ≤ c (1 + H^4/δ^4) (^0,^0)

for some constant c>0.



For every ^0∈ U_h^0 there exists a unique stream function s∈ S(Ω,), such that

    ^0 =  s.

Now we decompose s = ∑_j=0^J s_j and choose

    ^0_j =  s_j,

such that ^0_j∈ U_j^0 and by linearity of the curl

    ^0 = ∑_j=0^J ^0_j.

Continuity of (·,·), the stability of the decomposition in (<ref>), the identity (<ref>), 
and Korn's inequality (<ref>) yield

    (^0_0,^0_0) + ∑_j=1^J (^0_j,^0_j)
        ≤
        c ([H]s_0^2 + ∑_j=1^J s_j^2 )
    ≤ c (1+H^4/δ^4)s^2 
        = c (1+H^4/δ^4)  s^2
    ≤ c (1+H^4/δ^4) (^0,^0).



In order to show that also U_h^⊥ has a stable decomposition with respect to the norm induced by the elasticity bilinear form, we will adapt the construction from <cit.>.

For this purpose, we recall the definition of the spaces

    _0(Ω)    := {∈(Ω): · =0 ∂Ω}, 
    _0(Ω)    := { v ∈(Ω): v · =0 ∂Ω},

and decompose them into

    _0(Ω)    = H_0^÷,0(Ω) ⊕ H_0^÷,⊥(Ω), 
    _0(Ω)    = H_0^,0(Ω) ⊕ H_0^,⊥(Ω) ,

where H_0^÷,⊥(Ω) and H_0^,⊥(Ω) denote the orthogonal complements of the
divergence-free and curl-free subspaces

    H_0^÷,0(Ω)    := {∈_0(Ω): ÷ =0 },  
    
    H_0^,0(Ω)    := { v ∈_0(Ω):  v = 0},

of _0(Ω) and _0(Ω) in the (·,·)_(Ω) and (·,·)_(Ω) inner products, respectively.
Note that the decompositions (<ref>) and (<ref>) are both L^2-orthogonal as well.

Now, we introduce a new space U^+ containing semicontinuous functions,
which can be characterized as the image of U_h under an orthogonal projection
Θ^⊥_0(Ω)→ H_0^÷,⊥(Ω), i.e., U^+:=Θ^⊥(U_h), which is defined by

    Θ^⊥ :=  -  w,

where w∈ H_0^,⊥(Ω) satisfies

    ( w,  v) = (,  v),   ∀ v∈ H_0^,⊥(Ω).

Further, we define the projection P^h_0(Ω)→ U^+ by

    (÷(P^h -), ÷) = 0,   ∀∈ U^+.

Then, for P^h there holds

    P^h ^⊥≤ c ^⊥,   ∀^⊥∈ U_h^⊥

with a constant c depending on Ω.
According to <cit.> 
there is another constant c, independent of the mesh size, but depending on shape-regularity of the mesh and the polynomial degree, such that

    P^h ^⊥≤ c ^⊥.

Moreover, from <cit.> for convex Ω, we have the approximation property

    u^⊥-P^h ^⊥≤ c   h ÷^⊥    ∀^⊥∈ U_h^⊥.

For the L^2-projection [H] L^2(Ω)^d→ U_0 onto the coarse space U_0 there holds

    - [H] ≤ c   H .

Additionally, we have H^1-stability of [H], that is, 

    [H]≤ c .

We are ready to prove the following lemma.

 
Every function ^⊥∈ U_h^⊥ admits a decomposition of the form ^⊥=∑_j=0^J ^1_j with ^1_j∈ U_j,
which satisfies the bound

    (^1_0,^1_0) + ∑_j=1^J (^1_j,^1_j)
        ≤ c (1 + H^2/δ^2) (^⊥,^⊥)

for some constant c>0.



We choose

    ^1_0    = [H] P^h ^⊥, 
    ^1_j    = (θ_j (^⊥-^1_0)),

where  is the canonical interpolation into U_h and {θ_j}_j=1^J the piecewise linear partition of unity defined
in (<ref>).
We point at the fact that ^1_j∈ U_j is not necessarily a function in U_j^⊥.

Let = ^⊥-^1_0.
Noting that on each cell ∈, the function θ_j is polynomial, such that there holds

    (θ_j )≤ c θ_j .

By continuity of the bilinear form (·,·), the bounds (<ref>) and (<ref>),
and coercivity of (·,·), we get

    (^1_0,^1_0)
         ≤ c[H]P^h ^⊥^2 
       ≤ cP^h ^⊥^2
      ≤ c^⊥^2
      ≤ c (^⊥,^⊥)

and

    ^2
        ≤ c (^⊥^2 + [H] P^h ^⊥^2 )
        ≤ c (^⊥,^⊥).

By (<ref>), (<ref>), and using h≤ H, we obtain

    ≤^⊥ - P^h ^⊥ + P^h ^⊥ - [H]P^h ^⊥
       ≤ c   h ÷^⊥ + c   H ^⊥
       ≤ c   H ^⊥.

Using the continuity of (·,·), (<ref>), the properties of θ_j,
(<ref>), (<ref>), and the coercivity of (·,·), as well as the 
finite covering <ref>, we get

    ∑_j=1^J (^1_j,^1_j)
         ≤ c ∑_j=1^J (θ_j )^2
      ≤ c ∑_j=1^J θ_j ^2
      
       ≤ c ∑_j=1^J ( ∑_∈( [()]∇θ_j^2 []^2 + []θ_j ∇^2 )
      
               
             + ∑_∈1/h[]θ_j ^2
             + ∑_∈1/h[]θ_j ^2 )
      
       ≤ c 1/δ^2^2
      + c ^2
      
       ≤ c (1+H^2/δ^2) (^⊥,^⊥).

Finally, combining (<ref>) and (<ref>) results in the desired estimate (<ref>).


The last preparatory step is to show the stability of the decomposition of
[0]=⋃_j=0^J [0] with respect to (·,·).

 
Every ^0=(^0,^0)∈[0] admits a decomposition of the form
^0 = ∑_j=0^J ^0_j with ^0_j ∈[0], which satisfies the bound

    ∑_j=0^J (^0_j,^0_j)
        ≤(1+H^4/δ^4) (^0,^0),

where = c for some constant c>0 independent of the model parameters and of J.



By <ref>, every ^0∈ U_h^0 has a decomposition ^0=∑_J=0^J ^0_j with ^0_j∈ U_j^0, such that

    (^0_0,^0_0) + ∑_j=1^J (^0_j,^0_j) ≤ c (1 + H^4/δ^4) (^0,^0),

with ^0_j =  s_j, where s_j∈ S(Ω,) is chosen as in the proof of <ref>.

Further, by a classical result in domain decomposition theory that goes back to <cit.>, every r ∈ S(Ω,)
has a decomposition r=∑_j=0^J r_j, such that

    ∑_j=0^J r_j^2 ≤ c (1+H/δ) r^2,

where r = ∇ r.
By choosing ^0_j= r_j, we get for the second term in (<ref>)

    ∑_j=0^J (^0_j,^0_j)
           ≤∑_j=0^J  r_j^2
        = ∑_j=0^J r_j^2
        
       ≤ c  (1+H/δ) r^2
        = c  (1+H/δ)  r^2
        
       ≤ c  /(1+H/δ) (^0,^0)

where we have used (<ref>). The assertion of the lemma follows then
by choosing ^0_j = (^0_j, ^0_j) = ( s_j,  r_j) ∈[0] for all j=0,…,J.


We are ready to prove the main result of this section.


Every ∈ W_h admits a decomposition of the form = ∑_j=0^J _j
with _j∈, which satisfies the bound

    ∑_j=0^J (_j,_j)
        ≤(1+H^4/δ^4) (,),

where =c max, for some constant c>0
that is independent of the model parameters in the discrete bilinear form (·,·) defined in (<ref>) and independent of the number of subdomains J as well as of the discretization parameters h and τ.



For =  0 the result is trivial. So let ∈ W_h, 0, be arbitrary.
For the sake of simplicity, we indicate dependencies of constants by subscripts.
To start with, we decompose  in the form

    :=    = ^0^0 + ^⊥^⊥,

with (^0,^0)∈ U^0× V^0 and (^⊥,^⊥)∈ U^⊥× V^⊥.
Since we need to control the term (+,+), we 
define = ^⊥ + ^⊥
and decompose it into =^0+^⊥, ^0∈ V^0, where ^0 and ^⊥ are orthogonal with respect to
(·,·).
Then, the decomposition of  that we will use in our proof reads as


    =    = ^0^0 + ^⊥-^⊥ +  0
    
        =   ^0^0 +  0^0 + ^⊥-^⊥ +  0^⊥
        =: ^0+^1+^2+^3.


We decompose each component of each summand in (<ref>) according to

    ∑_j=0^J _j^0 = ^0,
              ∑_j=0^J _j^0 = ^0,
              ∑_j=0^J _j^0 = ^0, 
       ∑_j=0^J _j^1 = ^⊥,
              ∑_j=0^J _j^1 = ^⊥,
              ∑_j=0^J _j^1 = ^⊥,

where _j^0∈ U_j^0, _j^0∈ V_j^0, _j^0∈ V_j^0, _j^1∈ U_j, _j^1∈ V_j and _j^1∈ V_j are specified below. 
The superscript 1 of _j^1, _j^1 and _j^1 indicates that these terms are not orthogonal to _j^0, _j^0 and _j^0 with respect to the inner products (·,·), (·,·) and (·,·), respectively.
Now, we define _j:=_j^0 + _j^1 + _j^2 + _j^3, where

    _j^0    = _j^0_j^0∈ W_j^0,        _j^1    =  0_j^0∈ W_j^0, 
    _j^2    = _j^1-_j^1∈ W_j,    _j^3    =  0_j^1∈ W_j,

and herewith estimate

    ∑_j=0^J (_j,_j)
           ≤ c∑_j=0^J ((_j^0,_j^0) + (_j^1,_j^1) .
       . + (_j^2,_j^2) + (_j^3,_j^3) ).

For the decomposition in the orthogonal complement W_h^⊥, we will use the stability estimates

    (^0,^0)    ≤(,) = (^⊥+^⊥,^⊥+^⊥),  
    (^⊥,^⊥)    ≤(,) = (^⊥+^⊥,^⊥+^⊥)  .

Further, we note that

    (^⊥,^⊥) = (,) = (^⊥+^⊥,^⊥+^⊥)
        .

Then, the stability of the L^2-decomposition (<ref>) and the estimate (<ref>) yield

    ∑_j=0^J (_j^1,_j^1)
        = ∑_j=0^J    (_j^0,_j^0)
        ≤ c_H,δ,(^0,^0)
        
       ≤ c_H,δ,((^⊥,^⊥) + (^⊥,^⊥) ).

Moreover, by <ref> and <cit.>, we have

    ∑_j=0^J (_j^2,_j^2)
           = ∑_j=0^J ((_j^1,_j^1) + (_j^1,_j^1) + (_j^1,_j^1))
        
       ≤ c_H,δ((^⊥,^⊥) + (^⊥,^⊥) + (^⊥,^⊥) ).

Again, by <cit.>, stability estimate (<ref>), and equality (<ref>), 
we obtain

    ∑_j=0^J (_j^3,_j^3)
           = ∑_j=0^J ( (_j^1,_j^1) + ^-1(_j^1,_j^1) ) 
       ≤ c_H,δ( (^⊥,^⊥) + ^-1(^⊥,^⊥) ) 
       ≤ c_H,δ( (^⊥,^⊥) + (^⊥,^⊥) + ^-1(^⊥+^⊥,^⊥+^⊥) ).

Due to the special choice of the decomposition (<ref>) the term (^⊥,^⊥) appears in (<ref>),
which we further estimate using the Poincaré and Korn inequalities

    (^⊥,^⊥) ≤ c_Ω(^⊥,^⊥).


Thus, by (<ref>), <ref> and collecting the estimates
(<ref>), (<ref>), (<ref>), and
(<ref>), we obtain with ^⊥=(^⊥,^⊥)^T∈ W_h^⊥

    ∑_j=0^J    (_j,_j)
        
    ≤    c_H,δ,( (^0,^0) + (^0,^0)
        
          + (^⊥,^⊥) + (^⊥,^⊥) + (^⊥,^⊥)
        + ^-1(^⊥+^⊥,^⊥+^⊥) )
        
    ≤    c_H,δ,,Ω(
        (^0,^0) + (^⊥,^⊥)
        )
        
    
        =     c_H,δ,,Ω(,)

where the constant c_H,δ,,Ω for a proper overlap δ depends only on Ω and the bounds  and .





 §.§ Results for the singularly perturbed problem



We now prove convergence of the Schwarz methods for the singularly perturbed problem by applying the theory for symmetric positive definite problems as in <cit.>.
Therefore, we further establish a local stability estimate in <ref> and strengthened Cauchy-Schwarz inequalities in (<ref>), which together with the stable decomposition in <ref> and  <ref> on the covering of the domain prove the results of this section, i.e., <ref>.

Local stability, according to <cit.>, means that the estimates

    (_j,_j) ≤(_j,_j)     ∀_j∈

hold true for all 0≤ j≤ J with a constant >0.
Since we have chosen exact bilinear forms on the local spaces, i.e. (·,·) = (·,·) restricted to ,
local stability is trivially fulfilled with =1 for all 1≤ j≤ J.
Thus, it remains to consider the case j=0.


For each _0∈ W_0 it holds

    (_0,_0) ≤H/h[0](_0,_0).




On the coarse space, we have chosen a non-inherited form for the approximation
of the elasticity bilinear form, namely (·,·) instead of (·,·), which differs in the face
and boundary terms because of the different cell size H≥ h. Due to the continuity of
the coarse displacement functions _0 = _H ∈ U_0 on every coarse cell, the jump terms []_H
vanish for all faces  that lie in the interior of a coarse cell. Therefore, we have

    (_0,_0) ≤H/h(_0,_0).

Since all other terms of the coarse space operator [0](·,·)
are chosen to be exact, the statement of <ref> follows by estimate (<ref>).



The convergence theory for multiplicative Schwarz methods in <cit.> requires that 0<<2.
In a typical application scenario, the coarse mesh size might be H=2h resulting in a factor Hh=2 which violates the above condition.
Hence, an additional relaxation factor <1 has to be introduced in the definition of the coarse bilinear form, i.e.,
[0](·,·)=(·,·),
that scales (·,·) such that Hh<2.
However, we would like to point out that this scaling is not required in our experiments and we choose =1 in these cases.


By <ref> and standard arguments, see <cit.>, there hold the strengthened Cauchy-Schwarz inequalities

    (_i,_j)≤(_i,_i)^1/2(_j,_j)^1/2

for _i∈ W_i, _j∈, 1≤ i,j≤ J, and constants 0≤≤1.
Alongside <cit.> we denote the spectral radius of ℰ = {} by .

The maximal and minimal eigenvalues of  are characterized by the Rayleigh quotients

    _max()
           = sup_∈ W_h(,)/(,),
        
    _min()
           = inf_∈ W_h(,)/(,).

We are now ready to state and prove the main results of this section, <ref>.



The two-level additive Schwarz preconditioner (<ref>) for the singularly perturbed system (<ref>) satisfies

    _min()
        ≥((1+H^4/δ^4))^-1,
            _max()
        ≤H/h( + 1 ),

where the constant  depends on the bounds  and  from (<ref>), but is independent of the mesh size h and the other model parameters , and .
In particular, we obtain a uniform preconditioner under the assumption that H ≤ c h for some constant c.



  By <ref> we are able to apply the Schwarz convergence theory for symmetric positive definite problems, cf. <cit.>, to the singularly perturbed problem.
  In particular, the stable decomposition from <ref>, the local stability estimate from <ref>, and the strengthened Cauchy-Schwarz inequalities (<ref>) constitute the assumptions of the abstract condition number estimate in <cit.>. Hence, we can apply this result which concludes the proof.




Under the same assumptions as in <ref>, the multiplicative Schwarz method (<ref>) for the singularly perturbed problem converges with a contraction number depending on the bounds  and  from equation (<ref>), but independent of the mesh size h and the other model parameters , and .




  We invoke <cit.>, where again the stable decomposition from <ref>, the local stability estimate from <ref>, and the strengthened Cauchy-Schwarz inequalities (<ref>) constitute the assumptions. For the coarse grid space, we additionally refer to <ref>.





 §.§ Proofs of the main results





  <Ref> guarantees the boundedness of the spectrum of  with respect to the singularly perturbed problem, uniform with respect to the perturbation parameter ^-1 for any >0.
  Thus, the equivalence of the three-field formulation and the singularly perturbed problem yields the convergence proof.



  <Ref> guarantees the convergence of the multiplicative Schwarz method (<ref>) with respect to the singularly perturbed problem with a contraction number independent of the perturbation parameter ^-1, and as a consequence, uniform for any >0.
  Thus, the equivalence of the three-field formulation and the singularly perturbed problem yields the convergence proof.





§ NUMERICAL EXPERIMENTS





For testing the numerical performance of the two-level Schwarz
preconditioners defined in <ref>, we follow the
setting in <cit.> and construct a test case
with solution in
[H^1_0(Ω)]^2×_0(Ω)× L^2_0(Ω).
For a more in depth investigation with different kinds of boundary
conditions, as well as a comparison to multilevel algorithms we refer the
reader to <cit.>.
The implementation of the solvers is based on the finite element library  <cit.>.

Consider (<ref>) with isotropic permeability
= I and right hand side given by

    =
      [ 900 ∂_xϕ - ∂_y^3ϕ - ∂_x^2∂_yϕ; 900 ∂_yϕ + ∂_x^3ϕ + ∂_x∂_y^2ϕ ],     
      g = 900 ϕ - (900ϕ - 1),

where ϕ is defined on the square Ω=(0,1)×(0,1) by

    ϕ = x^2(x-1)^2y^2(y-1)^2.

Homogeneous Dirichlet boundary conditions are prescribed for the solid
displacement  as well as for the normal direction of the
seepage velocity  on the whole boundary ∂Ω, i.e.,

    =  0   on ∂Ω,    · = 0   on ∂Ω.

Here,  denotes the unit outward normal vector. The solution is
defined only up to an additive constant for the pressure. Thus, we
search for a mean-value free solution in the pressure component
satisfying

    ∫_Ω p dx=0.

Then, the unique solution to this system is given by

    =
      [  ∂_y ϕ; -∂_x ϕ ],     
      p = 900 ϕ - 1,      = -900∇ p.


The mesh consists of an equidistant subdivision of Ω⊂^2 into rectangular cells.
The system is discretized with triplets RT_k× RT_k× Q_k, k≥ 0, of equal order Raviart-Thomas and discontinuous cell-wise polynomial finite element functions, such that the matching condition

    ÷ RT_k = Q_k

is fulfilled.
To solve the algebraic equations we use the Krylov subspace method GMRES preconditioned by the multiplicative two-level Schwarz method (<ref>), which combines a multiplicative subdomain solver with a multiplicative application of the coarse grid solver.
As stopping criterion we use a reduction of the unpreconditioned starting residual

measured in the L^2-norm by a factor 10^8.
The subdomains Ω_j form vertex patches, which are unions of four cells for all interior vertices with an overlap of size δ=h.
Patches consisting of two cells for vertices at the boundary ∂Ω and single cells for vertices in corners are excluded from the set of patches, since we could not observe any beneficial effects of their inclusion for calculations with homogeneous boundary conditions, cf. <cit.>.
To realize the homogeneous boundary conditions on the boundary ∂Ω_j of the vertex patches, we build the patch matrices only out of interior degrees of freedom of each vertex patch Ω_j.
The coarse space is always assembled one level below the actual level, i.e., H=2h leading to a constant factor H/δ=2.

The coarse and local problems are solved using a singular value decomposition in order to deal properly with the one-dimensional kernel of their operators caused by the nonuniqueness of the pressure solution on each subdomain.



In <ref> (left side) we observe uniform convergence of the method in a set of sample calculations with RT_2× RT_2× Q_2 finite elements, from which we conclude that no additional relaxation of the coarse grid bilinear form is needed albeit required by the theory. This observation gave reason to <ref>.
Furthermore, throughout all calculations with polynomial degree k=2, we found that the iteration counts remain constant for a mesh size below h=1/32.
This holds true for all variations of the model parameters in the system, i.e., , ^-1, and .
A comparison calculation in <ref> (right-hand side) also confirms this for more refined meshes with the multilevel Schwarz analogue  of the two-level algorithm (<ref>).
For a definition of  let there be given a hierarchy of nested meshes [h_0]⊏⋯⊏[h_L] with corresponding finite element spaces X_h_0⊂⋯⊂ X_h_L for mesh sizes h_0≤…≤ h_L=h.
Here, the symbol ⊏ means that mesh [h_ℓ+1] results by a refinement of [h_ℓ].
Then,  is defined by

    = (I - ^MG) ^-1,
         ^MG
        = ∏_j=1^J_L(I - _j^L) ⋯∏_j=1^J_1(I - _j^1) (I - _0^0),

where J_ℓ denotes the number of subdomains on each level ℓ=0,…,L, and _j^ℓ are the projection operators of each subdomain on each level.
In the multilevel setting we observe about two iterations more than in the two-level case, which is due to the relaxed coarse scale solve.
Nevertheless, the numbers show the independence of the iteration counts with respect to the levels on a mesh with h=1/32 already.





Now, we turn our attention to an examination of the robustness of the method with respect to the scaled model parameters in the system, i.e., , ^-1, and .
In this test scenario we choose again polynomial degree k=2 and show the converged iteration counts for a fixed mesh with mesh size h=1/32.
The robustness of the multiplicative two-level Schwarz method with respect to varying the scaled parameters  and ^-1 is demonstrated in <ref>.
Here we have set the scaled storage capacity coefficient =0, which actually represents the most difficult case of the saddle point problem for the solver.
The iteration counts of GMRES are at most six for the whole range of 10^-6≤≤ 10^6, and 10^-6≤^-1≤ 10^6.
For small ^-1 the iteration numbers are even smaller, which indicates the dependence
on  that we have seen in the convergence proof in <ref>.
On the other hand, we do not observe an increase in iterations for large ^-1, which corresponds to small permeabilities.
The situation that ^-1 takes small values is also realistic and occurs if the product of μ and τ and the unscaled permeability is large,
cf. the scaling of the physical parameters in (<ref>).

Moreover, the scaled parameter  is defined as the ratio of the Lamé coefficients that might be small in certain cases or become large for incompressible materials.
For variations in , the table shows nearly constant iterations for each fixed ^-1; only for ^-1=10^2 there is a slight shift in the numbers when >10^2.

Variations in the scaled storage coefficient  do not effect the performance of the method for values ≤ 10^2 as can be seen in <ref>.
Only for very large values of > 10^2 the iterations increase very slightly before they drop to 1.
This result is presented for =^-1=1 only but holds true for all possible combinations of  and ^-1.

Note that in realistic calculations the scaled storage coefficient  might take values even larger than one if the poroelastic medium contains soft inclusions modeled by a storage capacity greater than zero, together with small values of the Biot-Willis constant α and/or a large Lamé coefficient μ, cf. (<ref>).



Last, we turn our attention to the lower order case for polynomial degrees k=0 and k=1.
In this situation, we observed some sort of stiffness of the method compared to polynomial degrees k ≥ 2.
In a sample calculation in the lowest-order case k=0 in <ref> with = 100, ^-1 = 1, and = 0 we need more mesh refinements for some parameter settings (down to h=1/128) until the iteration
counts remain unchanged.
And we notice a stronger influence of the variation in the permeability on the iteration counts as they grow faster for larger values of ^-1.
Nevertheless, the method shows good performance and reasonable robustness also in the lowest-order case, i.e., employing the RT_0× RT_0× Q_0 space triplet.




§ CONCLUSION


In this article, we proposed two-level Schwarz methods to solve -conforming discretizations of Biot's quasi-static consolidation model.
We rigorously proved convergence of the monolithic algorithms by transforming the three-field formulation into an equivalent singularly perturbed symmetric positive definite system, for which a stable decomposition and local stability bounds form the basis for applying the abstract Schwarz convergence theory for noninherited forms.
Finally, independence with respect to the perturbation parameter resulted in the proof of convergence, robust with respect to the model parameters under the assumption that the permeability tensor can be bounded in the sense of (<ref>).
The numerical tests demonstrated the performance and robustness of the multiplicative two-level Schwarz method for different polynomial degrees of the mixed finite element spaces and confirmed the theoretical results, even for nearly incompressible and low-permeable materials.
Beyond the theoretical results, a multilevel Schwarz approach showed comparable performance at lower cost.

The three-dimensional case as well as the investigation of high-frequency-high-contrast problems are the subjects of future work.




§ ACKNOWLEDGMENT


The third author acknowledges the funding by the German Science Fund (Deut­sche Forschungsgemeinschaft, DFG) as part of the project “Physics-oriented solvers for multicompartmental poromechanics” under grant number 456235063.





alpha


