# A scalable preconditioning framework for stabilized contact mechanics with
hydraulically active fractures

Andrea Franceschini<EMAIL_ADDRESS>Laura Gazzola
<EMAIL_ADDRESS>Massimiliano Ferronato
<EMAIL_ADDRESS>Department of Civil, Environmental and
Architectural Engineering, University of Padova, Padova, Italy

###### Abstract

A preconditioning framework for the coupled problem of frictional contact
mechanics and fluid flow in the fracture network is presented. The porous
medium is discretized using low-order continuous finite elements, with cell-
centered Lagrange multipliers and pressure unknowns used to impose the
constraints and solve the fluid flow in the fractures, respectively. This
formulation does not require any interpolation between different fields, but
is not uniformly inf-sup stable and requires a stabilization. For the
resulting $3\times 3$ block Jacobian matrix, we design scalable
preconditioning strategies, based on the physically-informed block
partitioning of the unknowns and state-of-the-art multigrid preconditioners.
The key idea is to restrict the system to a single-physics problem,
approximately solve it by an inner algebraic multigrid approach, and finally
prolong it back to the fully-coupled problem. Two different techniques are
presented, analyzed and compared by changing the ordering of the restrictions.
Numerical results illustrate the algorithmic scalability, the impact of the
relative number of fracture-based unknowns, and the performance on a real-
world problem.

###### keywords:

Scalable preconditioners , Contact mechanics , Darcy fracture flow

###### MSC:

[2010] 65F08 , 65N22 , 65N30 , 65N55

††journal: J. Comput. Phys.

##  1 Introduction

In recent years, attention has grown around novel technologies and
applications in the subsurface, like geothermal energy production [1, 2, 3],
hydraulic fracturing [4, 5, 6], CO2 sequestration [7, 8, 9] and underground
gas storage [10, 11, 12]. In these contexts, one of the key components is the
simultaneous simulation of frictional contact mechanics and fluid flow in
faults and fractures, which represent tightly coupled physical processes. In
fact, the aperture and slippage between the contact surfaces drive the fluid
flow in the fractures, while the pressure variation perturbs the stress state
in the surrounding medium and influences the contact mechanics itself. To
achieve the desired accuracy, large domains are usually required, with high
resolution representations of geological structures and their heterogeneous
properties [13, 14], and, specifically, of faults and fracture networks [15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32]. It is,
therefore, natural to have a growing demand towards the development of
sophisticated models of increasing size, which are computationally intensive
and require better and better performances. A key factor in this sense is the
linear solver, which is usually by far the most time-consuming component in a
real-world simulation [33, 34].

In this work, we analyze the simulation of frictional contact mechanics
coupled with the fluid flow in a fracture network and present a scalable and
efficient preconditioning framework for the linear system arising from the
discretization and linearization of the coupled problem. As to the
discretization approach, we elect to use the Discrete Fracture Model (DFM)
[18], i.e., an explicit representation of the fracture surfaces, while the
constraints are imposed with the aid of Lagrange multipliers [35, 36, 37, 38,
39]. As it is common in geological and reservoir simulations, we rely on low-
order finite elements for the mechanics and a cell-centered finite volume
scheme for the fluid flow. Lagrange multipliers are the contact forces acting
on the fracture surfaces as a cell-centered variable, thus sharing the same
representation as the fluid pressure field with no interpolation needed. The
details of this discretization scheme are described in [40]. This approach is
unstable in the Ladyzhenskaya-Babuška-Brezzi (LBB) sense, i.e., it does not
uniformly satisfy the inf-sup condition [41, Section 3.1], and requires a
stabilization. In this work, we use the global algebraic approach introduced
in the reference work [40]. The Jacobian matrix arising from the described
problem is non-symmetric with a $3\times 3$ block structure, which has to be
properly preconditioned to allow for a robust, scalable and efficient solution
with the aid of Krylov subspace solvers.

It is well known that iterative methods based on
projections/orthogonalizations onto Krylov subspaces [42] are in practice
mandatory to solve large and sparse linear systems deriving from the
discretization of PDEs, because they allow for a lower complexity, smaller
memory requirement, and better degree of algorithmic parallelism than direct
methods [43]. However, robustness, scalability and computational efficiency of
this class of methods is tightly connected with the choice of a proper
preconditioning technique [42]. Roughly speaking, preconditioners are
approximate applications of the system matrix inverse, and, from the algebraic
viewpoint, can be classified into three main categories: (i) incomplete
factorizations [44, 45, 46], (ii) approximate inverses [47, 48, 49, 50, 51,
52], and (iii) multilevel methods, i.e., domain decomposition [53, 54, 55, 56,
57] and multigrid-like techniques [58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68,
69, 70, 71]. A key feature for a modern preconditioning framework is the
algorithmic scalability, i.e., the ability to solve an increasingly refined
problem with an approximately constant number of iterations of the Krylov
solver. This property is particularly important in view of the development of
problems of increasing size by exploiting the availability of massively
parallel computational platforms. Incomplete factorizations and approximate
inverses can exhibit amazing performances, but do not have a linear complexity
with the system size. By distinction, multilevel methods can have a lower
performance on a single system, but are designed to be optimal with respect to
the scalability issue. Algebraic multigrid (AMG, [72]) is one of the most
effective multilevel approaches and consists of the complementary use of: (i)
a smoother that reduces high frequency errors, (ii) a coarse grid correction
that reduces low frequency errors, and (iii) restriction and interpolation
operators, to move from one grid to another. Starting from the original works,
e.g., [73], a wide range of multigrid approaches has appeared in the
literature, extending the applicability of this method, originally designed
for elliptic PDEs, to both non-symmetric [74, 75] and block matrices [76, 77,
78, 79, 80, 81]. Nonetheless, robustness and efficiency is still an open issue
for AMG whenever used as a black-box tool in problems with these algebraic
properties. The Jacobian matrix arising from the model considered herein is a
non-symmetric $3\times 3$ block matrix and, despite the available studies for
similar problems, none of them can be straightforwardly and effectively
applied to our case. In the context of geomechanical simulations, only a few
studies on $2\times 2$ block Jacobian systems [82, 83, 80] are found by the
authors.

The purpose of this work is to design a scalable preconditioning framework for
the $3\times 3$ block matrix arising from the coupled simulation of frictional
contact mechanics and fluid flow in the fracture network. The idea is to
exploit the inherent physics-based block subdivision and the scalability of
AMG techniques available from the literature. The full system is first
restricted to a single-physics problem, then approximately solved by AMG, and
finally prolonged back to the original size. According to the selected
restriction ordering, different approaches can be derived. In this work, we
consider two different options and investigate advantages and drawbacks in
order to find the most appropriate algorithm for real-world simulations. The
paper is organized as follows. Section 2 introduces the physical problem in
both the strong and weak forms, in order to understand the meaning and
features of each block of the Jacobian system. In Section 3, the
preconditioning framework is presented, with a detailed analysis of two
selected options. Finally, Section 4 presents a set of numerical results with
the aim of comparing the proposed approaches and investigating the algorithmic
scalability in both theoretical and real-world benchmarks. A few concluding
remarks close the paper.

##  2 Problem statement

We model the deformation of an open elastic domain
$\Omega\subset\mathbb{R}^{3}$, assuming quasi-static conditions and
infinitesimal strains within the open time interval
$\mathcal{T}=\left(0,t_{\max}\right]$. We denote by $\partial\Omega$ its
boundary, with $\overline{\Omega}=\Omega\cup\partial\Omega$, and
$\boldsymbol{n}_{\Omega}$ the outer normal vector to $\partial\Omega$, while a
set of internal boundaries $\Gamma=\cup_{i=1}^{n_{f}}\Gamma_{i}$ represents a
fracture network consisting of $n_{f}$ surfaces. The external boundary is
subdivided into two non-overlapping subsets, $\partial\Omega_{u}$ and
$\partial\Omega_{\sigma}$, where Dirichlet and Neumann boundary conditions
apply, respectively. Each fracture $\Gamma_{i}$ consists of two overlapping
surfaces, $\Gamma_{i}^{-}$ and $\Gamma_{i}^{+}$, with the orientation defined
by a unitary vector $\boldsymbol{n}_{i}$ orthogonal to the fracture plane. By
convention, we choose
$\boldsymbol{n}_{i}=\boldsymbol{n}_{i}^{-}=-\boldsymbol{n}_{i}^{+}$. The
pressure field is defined on the union $\Gamma$ of the two-dimensional (2D)
domains $\Gamma_{i}$, with $\partial\Gamma_{i}$ a one-dimensional (1D) curve
defining the boundary of each fracture and
$\overline{\Gamma}_{i}=\Gamma_{i}\cup\partial\Gamma_{i}$. The curve
$\partial\Gamma_{i}$ is subdivided into two non-overlapping subsets,
$\partial\Gamma_{i,p}$ and $\partial\Gamma_{i,q}$, where Dirichlet and Neumann
boundary conditions for the pressure field are imposed. The vector
$\boldsymbol{m}_{i}$ denotes the outer normal direction to
$\partial\Gamma_{i}$. The fluid is assumed to be incompressible, and body
forces and buoyancy effects are neglected. The projection of the stress tensor
$\boldsymbol{\sigma}$ along $\boldsymbol{n}_{i}$,
$\boldsymbol{t}=\boldsymbol{\sigma}\cdot\boldsymbol{n}_{i}^{-}=-\boldsymbol{\sigma}\cdot\boldsymbol{n}_{i}^{+}=(t_{N}\boldsymbol{n}_{i}+\boldsymbol{t}_{T})$,
is the traction vector over $\Gamma_{i}$, with $t_{N}$ and
$\boldsymbol{t}_{T}$ its normal and tangential component, respectively, with
respect to the fracture-local reference frame. The traction on $\Gamma_{i}$
controls the possible slipping and aperture of the fracture according to the
Coulomb frictional law. A schematic representation of the considered
conceptual framework is shown in Figure 1(a).

(a)

(b)

Figure 1: (a) Conceptual scheme of the elastic domain and embedded fracture
network. (b) Example of low-order discretization.

The strong form of the initial boundary value problem (IBVP) can be stated as
follows [84, 85, 86, 40]: given the fluid discharge
$q_{s}:\cup_{i=1}^{n_{f}}\Gamma_{i}\times\mathcal{T}\rightarrow\mathbb{R}$,
the prescribed boundary displacement
$\overline{\boldsymbol{u}}:\partial\Omega_{u}\times\mathcal{T}\rightarrow\mathbb{R}^{3}$
and traction
$\overline{\boldsymbol{t}}:\partial\Omega_{\sigma}\times\mathcal{T}\rightarrow\mathbb{R}^{3}$,
the prescribed fracture boundary pressure
$\overline{p}:\cup_{i=1}^{n_{f}}\partial\Gamma_{i,p}\times\mathcal{T}\rightarrow\mathbb{R}$
and flux
$\overline{q}:\cup_{i=1}^{n_{f}}\partial\Gamma_{i,q}\times\mathcal{T}\rightarrow\mathbb{R}$,
the initial displacement
$\boldsymbol{u}_{0}:\overline{\Omega}\rightarrow\mathbb{R}^{3}$ and pressure
$p_{0}:\cup_{i=1}^{n_{f}}\overline{\Gamma}_{i}\rightarrow\mathbb{R}$, find the
displacement
$\boldsymbol{u}:\overline{\Omega}\times[0,t_{\max}]\rightarrow\mathbb{R}^{3}$,
the traction
$\boldsymbol{t}:\cup_{i=1}^{n_{f}}\Gamma_{i}\times[0,t_{\max}]\rightarrow\mathbb{R}^{3}$,
and pressure
$p:\cup_{i=1}^{n_{f}}\overline{\Gamma}_{i}\times[0,t_{\max}]\rightarrow\mathbb{R}$
such that:

$\displaystyle-\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})$
$\displaystyle=\boldsymbol{0}$ in $\displaystyle\Omega\times\mathcal{T}$
$\displaystyle\mbox{(linear momentum balance)},$ (1a)
$\displaystyle\dot{g}_{N}(\boldsymbol{u})+\nabla\cdot\boldsymbol{q}(\boldsymbol{u},p)$
$\displaystyle=q_{s}$ in
$\displaystyle\cup_{i=1}^{n_{f}}\Gamma_{i}\times\mathcal{T}$
$\displaystyle\mbox{(mass balance)},$ (1b)
$\displaystyle\boldsymbol{t}-p\boldsymbol{n}_{i}$ $\displaystyle=0$ on
$\displaystyle\cup_{i=1}^{n_{f}}\Gamma_{i}\times\mathcal{T}$
$\displaystyle\mbox{(traction balance)},$ (1c) respecting the boundary
conditions $\displaystyle\boldsymbol{u}$ $\displaystyle=\bar{\boldsymbol{u}}$
on $\displaystyle\partial\Omega_{u}\times\mathcal{T}$
$\displaystyle\mbox{(prescribed boundary displacement)},$ (1d)
$\displaystyle\boldsymbol{\sigma}(\boldsymbol{u})\cdot\boldsymbol{n}_{\Omega}$
$\displaystyle=\bar{\boldsymbol{t}}$ on
$\displaystyle\partial\Omega_{\sigma}\times\mathcal{T}$
$\displaystyle\mbox{(prescribed boundary traction)},$ (1e) $\displaystyle p$
$\displaystyle=\bar{p}$ on
$\displaystyle\cup_{i=1}^{n_{f}}\partial\Gamma_{i,p}\times\mathcal{T}$
$\displaystyle\mbox{(prescribed boundary pressure)},$ (1f)
$\displaystyle\boldsymbol{q}(\boldsymbol{u},p)\cdot\boldsymbol{m}_{i}$
$\displaystyle=\bar{q}$ on
$\displaystyle\cup_{i=1}^{n_{f}}\partial\Gamma_{i,q}\times\mathcal{T}$
$\displaystyle\mbox{(prescribed boundary flux)},$ (1g) and initial conditions
$\displaystyle\boldsymbol{u}|_{t=0}$ $\displaystyle=\boldsymbol{u}_{0}$ in
$\displaystyle\overline{\Omega}$ $\displaystyle\mbox{(initial displacement)},$
(1h) $\displaystyle p|_{t=0}$ $\displaystyle=p_{0}$ in
$\displaystyle\cup_{i=1}^{n_{f}}\overline{\Gamma}_{i}$
$\displaystyle\mbox{(initial pressure)},$ (1i) subject to the constraints over
each $\Gamma_{i}$ and for every time in $\mathcal{T}$ $\displaystyle
t_{N}=\boldsymbol{t}\cdot\boldsymbol{n}_{i}$ $\displaystyle\leq 0,$
$\displaystyle
g_{N}=\llbracket\boldsymbol{u}\rrbracket\cdot\boldsymbol{n}_{i}\geq 0,\quad
t_{N}g_{N}$ $\displaystyle=0$ $\displaystyle\mbox{(normal contact
conditions)},$ (1j)
$\displaystyle\left\|\boldsymbol{t}_{T}\right\|_{2}-\tau_{\max}(t_{N})$
$\displaystyle\leq 0,$
$\displaystyle\dot{\boldsymbol{g}}_{T}\cdot\boldsymbol{t}_{T}-\tau_{\max}(t_{N})\left\|\dot{\boldsymbol{g}}_{T}\right\|_{2}$
$\displaystyle=0$ $\displaystyle\mbox{(Coulomb frictional law)}.$ (1k)

In the problem statement,
$\boldsymbol{\sigma}(\boldsymbol{u})=\textbf{C}:\nabla^{s}\boldsymbol{u}$ is
the Cauchy stress tensor, with C the fourth-order elasticity tensor;
$\boldsymbol{q}(\boldsymbol{u},p)=-(C_{f}(\boldsymbol{u})/\mu)\nabla p$ is the
fluid volumetric flux in the fracture domain according to Darcy’s law
[87]—assuming laminar flow—with $\nabla p$ the fluid pressure gradient, $\mu$
the fluid viscosity (constant), and $C_{f}$ the isotropic fracture hydraulic
conductivity modeled as in [18]:

$C_{f}=C_{f,0}+\frac{g_{N}^{3}}{12},$ (2)

with $C_{f,0}$ the conductivity related to two irregular surfaces that are in
contact [88];
$\llbracket\boldsymbol{u}\rrbracket=(\left.\boldsymbol{u}\right|_{\Gamma_{i}^{+}}-\left.\boldsymbol{u}\right|_{\Gamma_{i}^{-}})=(g_{N}\boldsymbol{n}_{i}+\boldsymbol{g}_{T})$
denotes the relative displacement across $\Gamma_{i}$, where $g_{N}$ and
$\boldsymbol{g}_{T}$ are the normal and tangential components, respectively,
and $\left.\boldsymbol{u}\right|_{\Gamma_{i}^{+}}$ and
$\left.\boldsymbol{u}\right|_{\Gamma_{i}^{-}}$ are the restrictions of
$\boldsymbol{u}$ on $\Gamma_{i}^{+}$ and $\Gamma_{i}^{-}$;
$\tau_{\text{max}}=c-t_{N}\tan(\theta)$ is the limit value provided by the
static Coulomb criterion, with $c$ and $\theta$ the cohesion and friction
angle, respectively. Since we employ a static Coulomb criterion, the
tangential velocity $\dot{\boldsymbol{g}}_{T}$ in (1k) is replaced with the
tangential displacement increment $\Delta{\boldsymbol{g}}_{T}$ [41] with
respect to the previously converged time-step.

In our framework, we assume $\Gamma$ to be fixed with no propagation. The
domain $\Gamma$ is partitioned into three portions, where the following
contact conditions occur:

  * 1.

stick on $\Gamma^{\text{stick}}$: the fracture is closed
($\llbracket\boldsymbol{u}\rrbracket=\boldsymbol{0}$) and the traction vector
$\boldsymbol{t}$ is unknown;

  * 2.

slip on $\Gamma^{\text{slip}}$: the fracture is closed in the normal direction
($g_{N}=0$ and $t_{N}$ is unknown), but a slip displacement
$\boldsymbol{g}_{T}$ between $\Gamma_{i}^{+}$ and $\Gamma_{i}^{-}$ is allowed
for, with
$\boldsymbol{t}_{T}|_{\Gamma^{\text{slip}}}=\boldsymbol{t}_{T}^{\ast}=\tau_{\text{max}}(t_{N})\frac{\boldsymbol{\Delta
g}_{T}}{||\boldsymbol{\Delta g}_{T}||_{2}}$;

  * 3.

open on $\Gamma^{\text{open}}$: the fracture is fully open and a free relative
displacement $\llbracket\boldsymbol{u}\rrbracket$ is allowed for, with
$\boldsymbol{t}=\boldsymbol{0}$.

For additional details regarding the governing formulation, we refer the
reader to [84, 85, 86, 40].

###  2.1 Discrete weak form

In the solution to the model problem (1), the traction $\boldsymbol{t}$ used
as a primary variable plays the role of Lagrange multipliers. Denoting with
$(\cdot,\cdot)_{D}$ the appropriate $L^{2}$-inner product of scalar, vector or
tensor functions in the spatial domain $D$, we introduce the finite-
dimensional subspaces $\boldsymbol{\mathcal{V}}^{h}$,
$\boldsymbol{\mathcal{M}}^{h}(t^{h}_{N})$ and $\mathcal{P}^{h}$:

$\displaystyle\boldsymbol{\mathcal{V}}^{h}\subset\boldsymbol{\mathcal{V}}$
$\displaystyle=\\{\boldsymbol{\eta}\in[H^{1}(\Omega)]^{3}:\boldsymbol{\eta}=\overline{\boldsymbol{u}}\mbox{
on }\partial\Omega_{u}\\},$ (3a)
$\displaystyle\boldsymbol{\mathcal{M}}^{h}(t^{h}_{N})\subset\boldsymbol{\mathcal{M}}(t_{N})$
$\displaystyle=\left\\{\boldsymbol{\mu}\in[L^{2}(\Gamma)]^{3}:\mu_{N}\leq
0,\left(\boldsymbol{\mu},\boldsymbol{v}\right)_{\Gamma}\leq\left(\tau_{\max}(t_{N}),||\boldsymbol{v}_{T}||_{2}\right)_{\Gamma},\boldsymbol{v}\in[H^{1/2}(\Gamma)]^{3}]\text{
with }v_{N}\geq 0\right\\},$ (3b)
$\displaystyle\mathcal{P}^{h}\subset\mathcal{P}$ $\displaystyle=\\{\chi\in
L^{2}(\Gamma)\\},$ (3c)

and the discrete approximations
$\\{\boldsymbol{u}^{h},\boldsymbol{t}^{h},p^{h}\\}$ of
$\\{\boldsymbol{u},\boldsymbol{t},p\\}$:

$\boldsymbol{u}^{h}=\sum_{i=1}^{n_{u}}\boldsymbol{\eta}_{i}u_{i}\in\boldsymbol{\mathcal{V}}^{h},\qquad\boldsymbol{t}^{h}=\sum_{j=1}^{n_{t}}\boldsymbol{\mu}_{j}t_{j}\in\boldsymbol{\mathcal{M}}^{h}(t^{h}_{N}),\qquad
p^{h}=\sum_{k=1}^{n_{p}}\chi_{k}p_{k}\in\mathcal{P}^{h},$ (4)

where, as before, the pedices $N$ and $T$ denote the components of a vector
function along the normal and tangential direction with respect to a fracture-
local reference frame on every $\Gamma_{i}$. In (4), $n_{u}$, $n_{t}$, and
$n_{p}$ denote the number of discrete displacement, traction and pressure
unknowns. The weak form of (1) reads [40]: find
$\\{\boldsymbol{u}^{h},\boldsymbol{t}^{h},p^{h}\\}\in\boldsymbol{\mathcal{V}}^{h}\times\boldsymbol{\mathcal{M}}^{h}(t^{h}_{N})\times\mathcal{P}^{h}$
such that

$\displaystyle\mathcal{R}_{u}$
$\displaystyle=(\nabla^{s}\boldsymbol{\eta},\boldsymbol{\sigma})_{\Omega}+\sum_{i=1}^{n_{f}}(\llbracket\boldsymbol{\eta}\rrbracket,\boldsymbol{t}^{h}-p^{h}\boldsymbol{n}_{i})_{\Gamma_{i}}-(\boldsymbol{\eta},\bar{\boldsymbol{t}})_{\partial\Omega_{\sigma}}=0$
$\displaystyle\forall\boldsymbol{\eta}\in\boldsymbol{\mathcal{V}}^{h}_{0},$
(5a) $\displaystyle\mathcal{R}_{t,i}$
$\displaystyle=(t_{N}^{h}-\mu_{N},g_{N})_{\Gamma_{i}}+(\boldsymbol{t}_{T}^{h}-\boldsymbol{\mu}_{T},\Delta{\boldsymbol{g}}_{T})_{\Gamma_{i}}\geq
0$
$\displaystyle\forall\boldsymbol{\mu}\in\boldsymbol{\mathcal{M}}^{h}(t^{h}_{N}),$
$\displaystyle i=1,\ldots,n_{f},$ (5b) $\displaystyle\mathcal{R}_{p,i}$
$\displaystyle=\left(\chi,\frac{\Delta g_{N}}{\Delta
t}\right)_{\Gamma_{i}}+[\chi,p^{h}]_{\mathcal{F}_{i}}-F_{\mathcal{F}_{i}}(\chi)+G_{\mathcal{F}_{i}}(\chi)-(\chi,q_{s})_{\Gamma_{i}}=0$
$\displaystyle\forall\chi\in\mathcal{P}^{h},$ $\displaystyle
i=1,\ldots,n_{f},$ (5c)

where $\boldsymbol{\mathcal{V}}^{h}_{0}$ is $\boldsymbol{\mathcal{V}}^{h}$
with homogeneous conditions along $\partial\Omega_{u}$, $\Delta t$ is the time
step size, and $[\chi,p^{h}]_{\mathcal{F}}$ is a weighted inner product
representing the classical two-point flux approximation (TPFA) scheme. This is
introduced to allow a unified presentation of the coupled finite
element/finite volume model [89, 90, 91]. In particular, we have:

$\displaystyle[\chi,p^{h}]_{\mathcal{F}}=\sum_{\mathcal{E}_{\text{int}}}(\chi|_{\varphi_{L}}-\chi|_{\varphi_{K}})\Upsilon_{KL}(p^{h}|_{\varphi_{L}}-p^{h}|_{\varphi_{K}})+\sum_{\mathcal{E}_{p}}\chi|_{\varphi_{K}}\Upsilon_{K}p^{h}|_{\varphi_{K}},$
(6)

where $\mathcal{E}_{\text{int}}$ and $\mathcal{E}_{p}$ represent the set of
edges included in $\Gamma$ and $\partial\Gamma_{p}$, respectively;
$\varphi_{K}$ and $\varphi_{L}$ are the two adjacent cells $K$ and $L$; and
$\Upsilon_{KL}$ is the harmonic average of one-sided transmissibility
$\Upsilon_{K}$ and $\Upsilon_{L}$ associated to $\varphi_{K}$ and
$\varphi_{L}$. Finally, $F_{\mathcal{F}}$ and $G_{\mathcal{F}}$ collect the
boundary conditions to be prescribed on $\partial\Gamma_{p}$ and
$\partial\Gamma_{q}$, respectively. For further details, we refer the reader
to [40].

To solve the problem (5), we transform the variational inequality (5b) into a
variational equality. For this purpose, we apply an active-set algorithm, as
described in [92, 93, 40], which allows to identify the subdivision into
stick/slip/open regions for every $\Gamma_{i}$. At a given step of the active-
set algorithm, the stick/slip/open regions of each fracture $\Gamma_{i}$ are
fixed and the inequality (5b) becomes:

$\mathcal{R}_{t,i}=\left(\boldsymbol{\mu},\boldsymbol{g}\right)_{\Gamma_{i}^{\text{stick}}}+\left(\mu_{N},g_{N}\right)_{\Gamma_{i}^{\text{slip}}}+\frac{1}{k}\left(\boldsymbol{\mu}_{T},\boldsymbol{t}_{T}^{h}-\boldsymbol{t}_{T}^{*}\right)_{\Gamma_{i}^{\text{slip}}}+\frac{1}{k}\left(\boldsymbol{\mu},\boldsymbol{t}^{h}\right)_{\Gamma_{i}^{\text{open}}}=0,\qquad
i=1,\ldots,n_{f},$ (7)

with $k$ a coefficient needed to ensure the dimensional consistency of the
equation. Introducing in (5a), (7) and (5c) the finite-dimensional bases of
$\boldsymbol{\mathcal{V}}^{h}_{0}$, $\boldsymbol{\mathcal{M}}^{h}(t_{N}^{h})$
and $\mathcal{P}^{h}$ yields the following system of nonlinear discrete
residual equations:

$\mathbf{r}(\mathbf{u}^{\ell},\mathbf{t}^{\ell},\mathbf{p}^{\ell})=\begin{bmatrix}\mathbf{r}_{u}(\mathbf{u}^{\ell},\mathbf{t}^{\ell},\mathbf{p}^{\ell})\\\
\mathbf{r}_{t}(\mathbf{u}^{\ell},\mathbf{t}^{\ell},\mathbf{p}^{\ell})\\\
\mathbf{r}_{p}(\mathbf{u}^{\ell},\mathbf{t}^{\ell},\mathbf{p}^{\ell})\end{bmatrix}=\mathbf{0},$
(8)

which is solved by a Newton-Krylov method. In (8), the algebraic vectors
$\mathbf{u}^{\ell}\in\mathbb{R}^{n_{u}}$,
$\mathbf{t}^{\ell}\in\mathbb{R}^{n_{t}}$ and
$\mathbf{p}^{\ell}\in\mathbb{R}^{n_{p}}$ collect the coefficients $u_{i}$,
$t_{j}$ and $p_{k}$ of the discrete displacement, traction and pressure fields
in (4) and $\ell$ is the active-set counter. After convergence of the Newton-
Krylov method at the $\ell$-th step of the active-set algorithm, a consistency
check is carried out in order to verify whether the assumed stick/slip/open
region subdivision meets the Coulomb frictional conditions. If not, the region
subdivision is updated and a new step is performed. The algorithm stops when
the consistency check does not require to modify the stick/slip/open region
subdivision. At this point, convergence is achieved and the solution is sought
at the following time step.

The finite element/finite volume spaces used in this work are the same as in
[40], i.e., first-order continuous finite elements for displacements and face-
centered piecewise-constant elements for tractions and pressures, as
schematically represented in Figure 1(b). To model the fractures, we use a DFM
approach with a conforming mesh [18], hence any $\Gamma_{i}$ is represented by
a set of finite element faces. Thus, displacement unknowns are located on mesh
vertices, while traction and pressure unknowns are on fracture faces (Figure
1(b)), with $n_{u}$ equal to three times the number of 3D finite element
nodes, $n_{p}$ equal to the number of 2D faces discretizing the fracture
network, and $n_{t}$ equal to $3\cdot n_{p}$. Displacements are represented in
the global reference system and tractions are represented in a face-based
local reference frame. This approach is intrinsically unstable, as it does not
fulfill the inf-sup condition [41]. In this work, we use the global algebraic
stabilization proposed in [40], which relaxes the zero jump and the
impenetrability conditions between the two fracture surfaces in the traction
balance equation, and the fluid incompressibility constraint in the mass
balance equation. Only stick and slip portions are involved in the traction
balance, being the tractions in the open part known. With the introduction of
the stabilization, equations (7) and (5c) become:

$\displaystyle\mathcal{R}_{t,i}$
$\displaystyle=\left(\boldsymbol{\mu},\boldsymbol{g}\right)_{\Gamma_{i}^{\text{stick}}}+\left(\mu_{N},g_{N}\right)_{\Gamma_{i}^{\text{slip}}}+\frac{1}{k}\left(\boldsymbol{\mu}_{T},\boldsymbol{t}_{T}^{h}-\boldsymbol{t}_{T}^{*}\right)_{\Gamma_{i}^{\text{slip}}}+\frac{1}{k}\left(\boldsymbol{\mu},\boldsymbol{t}^{h}\right)_{\Gamma_{i}^{\text{open}}}-j_{t}\left(\boldsymbol{\mu},\boldsymbol{t}^{h}\right)_{\Gamma_{i}^{\text{stick}}\cup\Gamma_{i}^{\text{slip}}}=0,$
$\displaystyle i=1,\ldots,n_{f},$ (9a) $\displaystyle\mathcal{R}_{p,i}$
$\displaystyle=\left(\chi,\frac{\Delta g_{N}}{\Delta
t}\right)_{\Gamma_{i}}+[\chi,p^{h}]_{\mathcal{F}_{i}}-F_{\mathcal{F}_{i}}(\chi)+G_{\mathcal{F}_{i}}(\chi)-(\chi,q_{s})_{\Gamma_{i}}+\frac{1}{\Delta
t}j_{p}\left(\chi,p^{h}\right)_{\Gamma_{i}}=0,$ $\displaystyle
i=1,\ldots,n_{f},$ (9b)

where $j_{t}(\boldsymbol{\mu},\boldsymbol{t}^{h})$ and $j_{p}(\chi,p^{h})$ are
the stabilizing bilinear forms for the traction and pressure field,
respectively. In particular, we have

$j_{t}\left(\boldsymbol{\mu},\boldsymbol{t}^{h}\right)=\sum_{\epsilon\in\mathcal{E}_{\text{int}}}\frac{1}{|\epsilon|}\int_{\epsilon}\llbracket\boldsymbol{\mu}\rrbracket_{\epsilon}\cdot\boldsymbol{\Upsilon}^{(\epsilon)}\cdot\llbracket\boldsymbol{t}^{h}\rrbracket_{\epsilon}\;\mathrm{d}l,$
(10)

with $\llbracket\cdot\rrbracket_{\epsilon}$ denoting the jump of a quantity
across the generic internal edge $\epsilon$ and
$\boldsymbol{\Upsilon}^{(\epsilon)}$ is a positive definite second-order
tensor providing the appropriate scaling. The discrete formulation of
$\boldsymbol{\Upsilon}^{(\epsilon)}$ is fully provided in [40]. The
contribution $j_{p}({\chi},{p}^{h})$ is computed as the normal projection of
$j_{t}$ with respect to the surface $\Gamma_{i}$.

At a given active-set iteration $\ell$, the Newton linearization of (8), which
now includes also the stabilization terms, generates a sequence of linear
systems and vector updates. To advance by one Newton iteration $k$, we have
to:

$\displaystyle\text{solve}\;\mathcal{J}^{\ell,(k)}\delta\mathbf{x}=-\mathbf{r}^{\ell,(k)}$
$\displaystyle\Rightarrow\qquad\begin{bmatrix}A&C_{1}&Q_{1}\\\ C_{2}&-H&0\\\
Q_{2}&0&T\\\ \end{bmatrix}^{\ell,(k)}\begin{bmatrix}\delta\mathbf{u}\\\
\delta\mathbf{t}\\\
\delta\mathbf{p}\end{bmatrix}=-\begin{bmatrix}\mathbf{r}_{u}\\\
\mathbf{r}_{t}\\\ \mathbf{r}_{p}\end{bmatrix}^{\ell,(k)},$ (11)
$\displaystyle\text{update}\;\mathbf{x}^{\ell,(k+1)}=\mathbf{x}^{\ell,(k)}+\delta\mathbf{x}$
$\displaystyle\Rightarrow\qquad\begin{bmatrix}\mathbf{u}\\\ \mathbf{t}\\\
\mathbf{p}\end{bmatrix}^{\ell,(k+1)}=\begin{bmatrix}\mathbf{u}\\\
\mathbf{t}\\\
\mathbf{p}\end{bmatrix}^{\ell,(k)}+\begin{bmatrix}\delta\mathbf{u}\\\
\delta\mathbf{t}\\\ \delta\mathbf{p}\end{bmatrix}.$

The submatrices in the $3\times 3$ block Jacobian read:

$\displaystyle[A]_{ij}$
$\displaystyle=\left(\nabla^{s}\boldsymbol{\eta}_{i},\textbf{C}:\nabla^{s}\boldsymbol{\eta}_{j}\right)_{\Omega},$
$\displaystyle i=1,n_{u},$ $\displaystyle j=1,n_{u},$ (12a)
$\displaystyle[C_{1}]_{ij}$
$\displaystyle=\left(\llbracket\boldsymbol{\eta}_{i}\rrbracket,\boldsymbol{\mu}_{j}\right)_{\Gamma},$
$\displaystyle i=1,n_{u},$ $\displaystyle j=1,n_{t},$ (12b)
$\displaystyle[Q_{1}]_{ij}$
$\displaystyle=-\left(\llbracket\eta_{i,N}\rrbracket,\chi_{j}\right)_{\Gamma},$
$\displaystyle i=1,n_{u},$ $\displaystyle j=1,n_{p},$ (12c)
$\displaystyle[C_{2}]_{ij}$
$\displaystyle=\left(\boldsymbol{\mu}_{i},\llbracket\boldsymbol{\eta}_{j}\rrbracket\right)_{\Gamma^{\ell,\text{stick}}}+\left({\mu}_{i,N},\llbracket{\eta}_{j,N}\rrbracket\right)_{\Gamma^{\ell,\text{slip}}}-\frac{1}{k}\left(\boldsymbol{\mu}_{i,T},\left(\frac{\partial\boldsymbol{t}_{T}^{*}}{\partial\Delta\boldsymbol{g}_{T}}\right)\biggr{|}^{\ell,(k)}\cdot\llbracket\boldsymbol{\eta}_{j,T}\rrbracket\right)_{\Gamma^{\ell,\text{slip}}},$
$\displaystyle i=1,n_{t},$ $\displaystyle j=1,n_{u},$ (12d)
$\displaystyle[H]_{ij}$
$\displaystyle=-\frac{1}{k}\left(\boldsymbol{\mu}_{i,T},\boldsymbol{\mu}_{j,T}\right)_{\Gamma^{\ell,\text{slip}}}+\frac{1}{k}\left(\boldsymbol{\mu}_{i,T},\left(\frac{\partial\boldsymbol{t}_{T}^{*}}{\partial
t_{N}}\right)\biggr{|}^{\ell,(k)}\mu_{j,N}\right)_{\Gamma^{\ell,\text{slip}}}-\frac{1}{k}\left(\boldsymbol{\mu}_{i},\boldsymbol{\mu}_{j}\right)_{\Gamma^{\ell,\text{open}}}+j_{t}\left(\boldsymbol{\mu}_{i},\boldsymbol{\mu}_{j}\right)_{\Gamma^{\ell,\text{stick}}\cup\Gamma^{\ell,\text{slip}}},$
$\displaystyle i=1,n_{t},$ $\displaystyle j=1,n_{t},$ (12e)
$\displaystyle[Q_{2}]_{ij}$ $\displaystyle=\frac{1}{\Delta
t}\left(\chi_{i},\llbracket\eta_{j,N}\rrbracket\right)_{\Gamma^{\ell,\text{open}}}+\left.\frac{\partial\left([\chi_{i},p^{h}]_{\mathcal{F}}\right)}{\partial
u_{j}}\right|^{\ell,(k)}-\left.\frac{\partial\left(F_{\mathcal{F}}(\chi_{i})\right)}{\partial
u_{j}}\right|^{\ell,(k)},$ $\displaystyle i=1,n_{p},$ $\displaystyle
j=1,n_{u},$ (12f) $\displaystyle[T]_{ij}$
$\displaystyle=\left.\frac{\partial\left([\chi_{i},p^{h}]_{\mathcal{F}}\right)}{\partial
p_{j}}\right|^{\ell,(k)}+\frac{1}{\Delta
t}j_{p}\left({\chi}_{i},{\chi}_{j}\right)_{\Gamma^{\ell}},$ $\displaystyle
i=1,n_{p},$ $\displaystyle j=1,n_{p}.$ (12g)

The partial derivatives appearing in (12) are reported in [40, Appendix A].

###  2.2 Linear system

We focus our attention on the linear system solution and the design of robust,
scalable and efficient preconditioners for the $3\times 3$ block matrix of
equation (11). The global matrix $\mathcal{J}$ is large, sparse, and non-
symmetric, with properties that change with the evolution of the
stick/slip/open regions in the fracture network. A representative evolution of
the non-zero pattern of $\mathcal{J}$ during a full simulation is shown in
Figure 2. The features that follow are worth summarizing.

  1. 1.

The first block row of $\mathcal{J}$ includes the contributions arising from
the linear momentum balance of the 3D domain $\Omega$. All the submatrices do
not depend on the fracture state and can be assembled once at the beginning of
the whole simulation if an elastic constitutive law is used. In particular,
$A$ is the classical symmetric positive definite (SPD) elastic stiffness
matrix, while $C_{1}$ and $Q_{1}$ are tall rectangular blocks collecting a
surface measure of the fracture elements and transferring tractions and
pressures to the 3D body as applied forces.

  2. 2.

In the second block row of $\mathcal{J}$, $C_{2}$ varies as the
stick/slip/open fracture regions evolve through the active-set algorithm, in
both the entry values and the non-zero pattern (Figure 2). If all the fracture
elements are in stick mode, we have that $C_{2}=C_{1}^{T}$, otherwise the
frictional law derivatives appear and $C_{2}\neq C_{1}^{T}$.

  3. 3.

When all fractures belong to the stick region, $H$ is the symmetric positive
semidefinite (SPSD) stabilization matrix. In case of sliding, non-symmetric
diagonal $2\times 2$ blocks arise, one for each traction component along the
local tangential direction to the fracture surface. In the open regions the
rows of $H$ have a single non-zero entry in the main diagonal, with no
contribution from the stabilization term (Figure 2). In any case, $H$ is
singular and cannot be regularly inverted.

  4. 4.

The third block row of $\mathcal{J}$ includes the contributions arising from
the fluid mass balance on the fracture network. The coupling between fluid
flow and fracture mechanics is controlled by $Q_{2}$. In particular, when all
fracture elements are in stick mode, $Q_{2}=0$ and $\mathcal{J}$ is reducible
with a $2\times 2$ symmetric saddle-point matrix as leading block. Otherwise,
contributions from the flux derivative with respect to the displacements
appear, i.e., $Q_{2}$ entries depend on the current pressure solution (Figure
2). By distinction with $C_{1}^{T}$ and $C_{2}$, there is no simple
relationship between $Q_{1}^{T}$ and $Q_{2}$, in both the entry values and the
non-zero pattern. Denoting with $\Psi$ the matrix-to-matrix operator returning
a zero row if the corresponding element index belongs to
$\Gamma^{\text{stick}}\cup\Gamma^{\text{slip}}$ and the original row if the
element index belongs to $\Gamma^{\text{open}}$, $Q_{2}$ can be written as:

$Q_{2}=-\frac{\Psi\left(Q_{1}^{T}\right)}{\Delta t}+F_{u},$ (13)

where $F_{u}$ collects the contributions from the flux derivatives with
respect to the displacements.

  5. 5.

$T$ is the sum of the standard transmissibility matrix arising from the TPFA
discretization in the 2D domain $\Gamma$ and the stabilization contribution.
As such, it is SPD with the 5-point stencil of a 2D discrete Laplacian.
Moreover, $T$ has a block diagonal structure for all non-intersecting
fractures. Observe also that traction and pressure fields are always
decoupled.

From the observations above, it appears that matrix $\mathcal{J}$ changes
nature with the evolution of the fracture conditions, moving from a reducible
matrix with a symmetric saddle-point leading block to a general non-symmetric
and indefinite matrix. The objective of our work is to define a unique
preconditioning framework ensuring robustness, scalability and computational
efficiency for any working situation.

(a)

(b)

(c)

Figure 2: Non-zero pattern of $\mathcal{J}$ with the evolution of the
stick/slip/open fracture regions.

##  3 Preconditioning framework

A preconditioner $\mathcal{M}^{-1}$ of $\mathcal{J}$ is a non-singular
operator whose application to a vector resembles as much as possible the
action of $\mathcal{J}^{-1}$. The exact application of $\mathcal{J}^{-1}$ to
some vector $\mathbf{w}\in\mathbb{R}^{n_{u}+n_{t}+n_{p}}$ provides the vector
$\mathbf{v}\in\mathbb{R}^{n_{u}+n_{t}+n_{p}}$ such that:

$\left\\{\begin{array}[]{ccccccc}A\mathbf{v}_{u}&+&C_{1}\mathbf{v}_{t}&+&Q_{1}\mathbf{v}_{p}&=&\mathbf{w}_{u}\\\
C_{2}\mathbf{v}_{u}&-&H\mathbf{v}_{t}&&&=&\mathbf{w}_{t}\\\
Q_{2}\mathbf{v}_{u}&&&+&T\mathbf{v}_{p}&=&\mathbf{w}_{p}\end{array}\right.,$
(14)

with $\mathbf{v}_{u},\mathbf{w}_{u}\in\mathbb{R}^{n_{u}}$,
$\mathbf{v}_{t},\mathbf{w}_{t}\in\mathbb{R}^{n_{t}}$, and
$\mathbf{v}_{p},\mathbf{w}_{p}\in\mathbb{R}^{n_{p}}$ natural subvectors of
$\mathbf{v},\mathbf{w}$, respectively. The objective is to approximate the
solution to the multi-physics system (14) by exploiting the physics-based
variable partitioning. The system is first reduced to a single-physics
problem, and then prolonged back to the full multi-physics space. According to
the selected sequence of reductions, different algorithms may arise.

###  3.1 Method no. 1: t-p-u approach

Traction and pressure variables live on the fractures and are mutually
decoupled independently on the stick/slip/ open region partitioning.
Therefore, it is natural to exploit this condition and perform a simultaneous
reduction of both variable sets onto the displacement space. This corresponds
to compute $\mathbf{v}_{t}$ and $\mathbf{v}_{p}$ from the second and third
equation of (14), respectively, and introduce them in the first equation, thus
eliminating both physics from the equilibrium equation. Recall, however, that
$H$ is singular, so a regular surrogate is needed. A block diagonal
approximation can be used instead, where each block is the $3\times 3$ local
stabilization matrix computed for each fracture element. Denoting with
$\widetilde{H}$ such a block-diagonal approximation, we have:

$\displaystyle\mathbf{v}_{t}$
$\displaystyle\simeq-\widetilde{H}^{-1}\left(\mathbf{w}_{t}-C_{2}\mathbf{v}_{u}\right),$
(15a) $\displaystyle\mathbf{v}_{p}$
$\displaystyle=T^{-1}\left(\mathbf{w}_{p}-Q_{2}\mathbf{v}_{u}\right).$ (15b)

With (15a) and (15b), the first equation of (14) becomes:

$\left(A+C_{1}\widetilde{H}^{-1}C_{2}-Q_{1}T^{-1}Q_{2}\right)\mathbf{v}_{u}\simeq\mathbf{w}_{u}+C_{1}\widetilde{H}^{-1}\mathbf{w}_{t}-Q_{1}T^{-1}\mathbf{w}_{p},$
(16)

which is a single-physics equilibrium equation on the 3D domain where the
elimination of fracture tractions and pressures introduces fictitious
stiffness contributions. The matrix at the left-hand side of (16) is the Schur
complement $S$:

$S=A+C_{1}\widetilde{H}^{-1}C_{2}-Q_{1}T^{-1}Q_{2}.$ (17)

Solution to (16) provides $\mathbf{v}_{u}$, which, introduced into equations
(15), yields the final vector $\mathbf{v}$. The multi-physics reduction order
performed in this case is traction-pressure-displacement (t-p-u) and is
schematically summarized in Figure 3.

Figure 3: Schematic representation of the multi-physics reduction
preconditioning framework: t-p-u (left) and t-u-p (right) approach.

The computation and inversion of $S$ in (17) cannot be performed exactly. The
Schur complement is explicitly approximated by $\widetilde{S}$:

$S\simeq\widetilde{S}=A+C_{1}\widetilde{H}^{-1}C_{2}-Q_{1}\widetilde{T}^{-1}Q_{2},$
(18)

where $\widetilde{T}$ is a diagonal surrogate for $T$. The inverse
$\widetilde{S}^{-1}$ can be applied inexactly by means of an AMG operator,
which can be efficiently used in mechanical problems preserving a linear
complexity with respect to the problem size. This is a key property to
guarantee the solver scalability. Recent examples of effective AMG
preconditioners are, for instance, taken from the References [67, 69, 70, 71].
In this work, we use an aggregation-based multigrid as the reference AMG
operator. Specifically, the application of $\widetilde{S}^{-1}$ is
approximated by GAMG [94], the state-of-the-art aggregation based multigrid
provided by the PETSc package [95].

The construction and application of the resulting preconditioning operator
$\mathcal{M}_{1}^{-1}$ with the t-p-u approach is summarized in Algorithms 1
and 2. The $\text{Bdiag}(K,n)$ operator gives a matrix with the $n\times n$
diagonal blocks of $K$, while the $\text{AMG}(K,\mathbf{y})$ operator applies
the selected AMG preconditioner of $K$ to the vector $\mathbf{y}$. Since
$n_{u}$ is generally much larger than $n_{t}$ and $n_{p}$, the cost for
applying the exact inverse of $\widetilde{H}$ and $T$ is negligible with
respect to the AMG algorithm for $\widetilde{S}$. Hence, the latter can be
roughly assumed as the cost per iteration for the $\mathcal{M}_{1}^{-1}$
application.

Algorithm 1 Preconditioner Construction: t-p-u approach [$\widetilde{H}$,
$\widetilde{S}$]=cpt_tpu($\mathcal{J}$).

1:$\widetilde{H}=\text{Bdiag}(H,3)$;

2:$\widetilde{T}=\text{diag}(T)$;

3:$\widetilde{S}=A+C_{1}\widetilde{H}^{-1}C_{2}-Q_{1}\widetilde{T}^{-1}Q_{2}$;

Algorithm 2 Preconditioner Application: t-p-u approach [$\mathbf{v}_{u}$,
$\mathbf{v}_{t}$, $\mathbf{v}_{p}$]=app_tpu($\mathcal{J}$, $\widetilde{H}$,
$\widetilde{S}$, $\mathbf{w}_{u}$, $\mathbf{w}_{t}$, $\mathbf{w}_{p}$).

1:$\mathbf{t}_{p}=T^{-1}\mathbf{w}_{p}$;

2:$\mathbf{t}_{t}=\widetilde{H}^{-1}\mathbf{w}_{t}$;

3:$\mathbf{t}_{u}=\mathbf{w}_{u}+C_{1}\mathbf{t}_{t}-Q_{1}\mathbf{t}_{p}$;

4:$\mathbf{v}_{u}=\text{AMG}(\widetilde{S},\mathbf{t}_{u})$;

5:$\mathbf{s}_{p}=Q_{2}\mathbf{v}_{u}$;

6:$\mathbf{s}_{t}=C_{2}\mathbf{v}_{u}$;

7:$\mathbf{v}_{p}=\mathbf{t}_{p}-T^{-1}\mathbf{s}_{p}$;

8:$\mathbf{v}_{t}=-\mathbf{t}_{t}+\widetilde{H}^{-1}\mathbf{s}_{t}$;

From an algebraic viewpoint, the preconditioning operator
$\mathcal{M}_{1}^{-1}$ arising from the t-p-u approach can be written as an
inexact block LDU factorization of $\mathcal{J}$. Using the permutation matrix
$\mathcal{Q}_{1}$:

$\mathcal{Q}_{1}=\begin{bmatrix}0&I_{t}&0\\\ 0&0&I_{p}\\\ I_{u}&0&0\\\
\end{bmatrix},$ (19)

where $I_{n}$ is the identity matrix in $\mathbb{R}^{n_{n}\times n_{n}}$ and
$0$ the zero matrix of proper size, the block LDU factorization reads:

$\mathcal{Q}_{1}\mathcal{J}\mathcal{Q}_{1}^{T}\simeq\mathcal{L}_{1}\mathcal{D}_{1}\mathcal{U}_{1},$
(20)

with:

$\displaystyle\mathcal{L}_{1}$ $\displaystyle=\begin{bmatrix}I_{t}&0&0\\\
0&I_{p}&0\\\ -C_{1}\widetilde{H}^{-1}&Q_{1}T^{-1}&I_{u}\\\ \end{bmatrix},$
$\displaystyle\mathcal{D}_{1}$
$\displaystyle=\begin{bmatrix}-\widetilde{H}&0&0\\\ 0&T&0\\\
0&0&\widetilde{S}\\\ \end{bmatrix}$ and $\displaystyle\mathcal{U}_{1}$
$\displaystyle=\begin{bmatrix}I_{t}&0&-\widetilde{H}^{-1}C_{2}\\\
0&I_{p}&T^{-1}Q_{2}\\\ 0&0&I_{u}\\\ \end{bmatrix}.$ (21)

Hence, the final algebraic expression of $\mathcal{M}_{1}^{-1}$ is:

$\mathcal{M}_{1}^{-1}=\begin{bmatrix}I_{t}&0&\widetilde{H}^{-1}C_{2}\\\
0&I_{p}&-T^{-1}Q_{2}\\\ 0&0&I_{u}\\\
\end{bmatrix}\begin{bmatrix}-\widetilde{H}^{-1}&0&0\\\ 0&T^{-1}&0\\\
0&0&\widetilde{S}^{-1}\\\ \end{bmatrix}\begin{bmatrix}I_{t}&0&0\\\
0&I_{p}&0\\\ C_{1}\widetilde{H}^{-1}&-Q_{1}T^{-1}&I_{u}\\\ \end{bmatrix}.$
(22)

######  Remark 1.

The multi-physics reduction approach proposed herein can be equivalently
recast in other ways as well. Since we use a twofold approximation for
$T^{-1}$, i.e., exact in $\mathcal{L}_{1}$, $\mathcal{U}_{1}$ and inexact in
$\widetilde{S}$, $\mathcal{M}_{1}^{-1}$ can be regarded as a member of the
mixed constraint preconditioner class [96, 97, 98]. Similarly, the upper and
lower block triangular factors in (22) play the role of decoupling operators
for the original multi-physics problem and are the outcome of the general-
purpose algebraic procedure defined in [99]. Finally, $\mathcal{M}_{1}^{-1}$
can be also regarded as an example of application in a $3\times 3$ block non-
symmetric context of the multigrid reduction framework [100, 101], where
fracture and body variables play the role of fine and coarse nodes,
respectively, and $\widetilde{H}$ replaces $H$ in matrix $\mathcal{J}$.

Let us introduce the matrices:

$\displaystyle E_{H}$ $\displaystyle=I_{t}-\widetilde{H}^{-1}H,$ (23a)
$\displaystyle E_{S}$ $\displaystyle=I_{u}-\widetilde{S}^{-1}S,$ (23b)

which can be regarded as a matrix measure of the quality of the approximations
$\widetilde{H}$ and $\widetilde{S}$ introduced in $\mathcal{M}_{1}^{-1}$. The
following result holds.

######  Proposition 3.1.

The eigenvalues $\lambda_{1}\in\mathbb{C}$ of the preconditioned matrix with
the t-p-u approach are either 1, with multiplicity $n_{p}$, or such that:

$\left|\lambda_{1}-1\right|\leq\left(1+\zeta\right)\max\\{\varepsilon_{H},\varepsilon_{S}\\},$
(24)

with $\varepsilon_{H}=\|E_{H}\|$, $\varepsilon_{S}=\|E_{S}\|$,
$\zeta=\|\mathcal{Z}\|$ and

$\mathcal{Z}=\left[\begin{array}[]{cc}\widetilde{H}^{-1}C_{2}\widetilde{S}^{-1}C_{1}&-\widetilde{H}^{-1}C_{2}\\\
\widetilde{S}^{-1}C_{1}&0\end{array}\right]$ (25)

for any compatible matrix norm.

###### Proof.

Recalling equations (19) and (22) and introducing the error matrices (23), the
preconditioned matrix with the t-p-u approach reads:

$\mathcal{M}_{1}^{-1}\mathcal{Q}_{1}\mathcal{J}\mathcal{Q}_{1}^{T}=\left[\begin{array}[]{ccc}I_{t}+\left(\widetilde{H}^{-1}C_{2}\widetilde{S}^{-1}C_{1}-I_{t}\right)E_{H}&0&-\widetilde{H}^{-1}C_{2}E_{S}\\\
-T^{-1}Q_{2}\widetilde{S}^{-1}C_{1}E_{H}&I_{p}&T^{-1}Q_{2}E_{S}\\\
\widetilde{S}^{-1}C_{1}E_{H}&0&I_{u}-E_{S}\end{array}\right],$ (26)

which has $n_{p}$ unitary eigenvalues. The remaining $n_{t}+n_{u}$ eigenvalues
are those of the matrix $\mathcal{H}$ obtained by dropping the second block
row and column from (26):

$\mathcal{H}=I_{t+u}+\left(\mathcal{Z}-I_{t+u}\right)\mathcal{E}_{1},$ (27)

with $I_{t+u}$ the identity matrix of order $n_{t}+n_{u}$ and
$\mathcal{E}_{1}$:

$\mathcal{E}_{1}=\left[\begin{array}[]{cc}E_{H}&0\\\
0&E_{S}\end{array}\right].$ (28)

The eigenvalues of $\mathcal{H}$ satisfy the bound (24), thus closing the
proof. ∎

######  Remark 2.

Proposition 3.1 shows that the distance of $\widetilde{H}^{-1}H$ from the
identity and the approximation quality of $\widetilde{S}^{-1}$ are key factors
for the overall performance of $\mathcal{M}_{1}^{-1}$. While $\widetilde{H}$
is fixed, notice, however, that $\widetilde{S}$ of equation (18) has algebraic
properties that change with the fracture state and the evolution of the
stick/slip/open regions throughout the active-set algorithm. In stick mode,
the contribution $C_{1}\widetilde{H}^{-1}C_{2}$ is symmetric positive
semidefinite and $Q_{2}=0$. Hence, $\widetilde{S}$ is SPD. Also in slip mode
the contribution $C_{1}\widetilde{H}^{-1}C_{2}$ is positive definite and
$Q_{2}=0$, so $\widetilde{S}$ remains positive definite, though slightly non-
symmetric. With open elements, however, $Q_{2}$ depends on the current
pressure solution and no theoretical considerations can be made in general. In
these conditions, $\widetilde{S}$ is an indefinite non-symmetric matrix.

###  3.2 Method no. 2: t-u-p approach

An alternative multi-physics reduction sequence relies on the scheme sketched
in the rightmost panel of Figure 3. Introducing the traction variables (15a)
into the first equation of system (14) yields:

$\left(A+C_{1}\widetilde{H}^{-1}C_{2}\right)\mathbf{v}_{u}+Q_{1}\mathbf{v}_{p}\simeq\mathbf{w}_{u}+C_{1}\widetilde{H}^{-1}\mathbf{w}_{t},$
(29)

where the matrix:

$S_{1}=A+C_{1}\widetilde{H}^{-1}C_{2}$ (30)

is the first-level Schur complement. From a physical viewpoint, $S_{1}$ is an
elasticity matrix with fictitious stiffness contributions arising along the
fractures from the traction elimination. Then, a second reduction is needed by
computing $\mathbf{v}_{u}$ from (29) and introducing it in the third equation
of (14):

$\left(T-Q_{2}S_{1}^{-1}Q_{1}\right)\mathbf{v}_{p}\simeq\mathbf{w}_{p}-Q_{2}S_{1}^{-1}\left(\mathbf{w}_{u}+C_{1}\widetilde{H}^{-1}\mathbf{w}_{t}\right).$
(31)

The matrix at the left-hand side of equation (31) is the second-level Schur
complement:

$S_{2}=T-Q_{2}S_{1}^{-1}Q_{1},$ (32)

which, from a physical point of view, represents a modified transmissibility
matrix including the effect of the stiffness of the 3D medium surrounding the
fractures. Hence, the multi-physics reduction order is traction-displacement-
pressure (t-u-p).

The computation of $S_{1}$ can be performed exactly, but its inverse has to be
approximated. Since its nature is the same as that of $S$ of equation (17), we
can effectively use an AMG operator, such as GAMG. We denote with
$\widetilde{S}_{1}^{-1}$ the operator that approximately applies $S_{1}^{-1}$.
By distinction, $S_{2}$ cannot be computed exactly. Recalling the physical
interpretation of $T$ and $S_{1}$, we can approximate the contribution
$Q_{2}S_{1}^{-1}Q_{1}$ by the diagonal fixed-stress matrix introduced as a
preconditioner in [102, 103]. Denoting by $S_{K}$ such a matrix, we
approximate $S_{2}$ as:

$S_{2}\simeq\widetilde{S}_{2}=T-S_{K},$ (33)

where the diagonal entries of $S_{K}$ are:

$[S_{K}]_{k}=\frac{|\Omega_{k}|}{\overline{K}_{k}},\qquad k=1,\ldots,n_{p},$
(34)

with $|\Omega_{k}|$ a measure of the volume of the cells surrounding the
$k$-th fracture element and $\overline{K}_{k}$ an estimate of the associated
bulk modulus.

######  Remark 3.

The computation of $|\Omega_{k}|$ and $\overline{K}_{k}$ in equation (34) can
be carried out using the information arising from the discretization grid and
the material properties. A more general algebraic strategy, however, can be
implemented by following the ideas sketched in [103]. Recalling the
definitions of $Q_{1}$ and $Q_{2}$, provided in (12c) and (12f), respectively,
it can be observed that the $k$-th column $Q_{1}^{(k)}$ of $Q_{1}$ and the
$k$-th row $Q_{2}^{(k),T}$ of $Q_{2}$, $k=1,\ldots,n_{p}$, are sparse vectors
with non-zero entries only in a small number of components corresponding to
the indices of the degrees of freedom associated to some nodes of the 3D cells
surrounding the $k$-th 2D fracture element. Let us denote with
$\mathcal{K}_{1}^{(k)}$ and $\mathcal{K}_{2}^{(k)}$ the subsets of
$\mathcal{K}=\\{1,\ldots,n_{u}\\}\subset\mathbb{N}$ containing the indices of
the non-zero components of $Q_{1}^{(k)}$ and $Q_{2}^{(k),T}$, respectively. In
general, $\mathcal{K}_{1}^{(k)}\neq\mathcal{K}_{2}^{(k)}$, with
$\mathcal{K}_{2}^{(k)}$ possibly being the empty set. Setting
$\mathcal{K}^{(k)}=\mathcal{K}_{1}^{(k)}\cup\mathcal{K}_{2}^{(k)}$, we define
$R^{(k)}$ as the restriction operator from $\mathbb{R}^{n_{u}}$ to
$\mathbb{R}^{|\mathcal{K}^{(k)}|}$ such that:

$S_{1}^{(k)}=R^{(k),T}S_{1}R^{(k)}$ (35)

is the sub-matrix of $S_{1}$ made by the entries lying in the rows and columns
with indices in $\mathcal{K}^{(k)}$. Since $S_{1}^{(k)}$ is a diagonal block
of an SPD matrix, it is non-singular and can be regularly inverted. The $k$-th
diagonal entry of $S_{K}$ can be therefore computed as:

$[S_{K}]_{k}=Q_{2}^{(k),T}S_{1}^{(k),-1}Q_{1}^{(k)},\qquad k=1,\ldots,n_{p},$
(36)

which provides a fully algebraic interpretation of the classical fixed-stress
contribution (34).

The construction and application of the resulting preconditioning operator
$\mathcal{M}_{2}^{-1}$ with the t-u-p approach is summarized in Algorithms 3
and 4. The FS($Q_{2},S_{1},Q_{1}$) operator gives the diagonal matrix $S_{K}$
with entries computed as in equation (36). Since $n_{u}$ is generally
significantly larger than $n_{t}$ and $n_{p}$, the cost for applying the exact
inverse of $\widetilde{H}$ and $\widetilde{S}_{2}$ is negligible with respect
to the AMG algorithm for $\widetilde{S}_{1}^{-1}$. Therefore, the cost per
iteration for the $\mathcal{M}_{2}^{-1}$ application is roughly two AMG calls,
i.e., twice that of $\mathcal{M}_{1}^{-1}$.

Algorithm 3 Preconditioner Construction: t-u-p approach [$\widetilde{H}$,
$S_{1}$, $\widetilde{S}_{2}$]=cpt_tup($\mathcal{J}$).

1:$\widetilde{H}=\text{Bdiag}(H,3)$;

2:$S_{1}=A+C_{1}\widetilde{H}^{-1}C_{2}$;

3:$\widetilde{S}_{2}=T-\text{FS}(Q_{2},S_{1},Q_{1})$;

Algorithm 4 Preconditioner Application: t-u-p approach [$\mathbf{v}_{u}$,
$\mathbf{v}_{t}$, $\mathbf{v}_{p}$]=app_tup($\mathcal{J}$, $\widetilde{H}$,
$S_{1}$, $\widetilde{S}_{2}$, $\mathbf{w}_{u}$, $\mathbf{w}_{t}$,
$\mathbf{w}_{p}$).

1:$\mathbf{t}_{t}=\widetilde{H}^{-1}\mathbf{w}_{t}$;

2:$\mathbf{t}_{u}=\mathbf{w}_{u}+C_{1}\mathbf{t}_{t}$;

3:$\mathbf{s}_{u}=\text{AMG}(S_{1},\mathbf{t}_{u})$;

4:$\mathbf{t}_{p}=\mathbf{w}_{p}-Q_{2}\mathbf{s}_{u}$;

5:$\mathbf{v}_{p}=\widetilde{S}_{2}^{-1}\mathbf{t}_{p}$;

6:$\mathbf{t}_{u}=Q_{1}\mathbf{v}_{p}$;

7:$\mathbf{v}_{u}=\text{AMG}(S_{1},\mathbf{t}_{u})$;

8:$\mathbf{v}_{u}\leftarrow\mathbf{s}_{u}-\mathbf{v}_{u}$;

9:$\mathbf{s}_{t}=C_{2}\mathbf{v}_{u}$;

10:$\mathbf{v}_{t}=-\mathbf{t}_{t}+\widetilde{H}^{-1}\mathbf{s}_{t}$;

Similarly to the t-p-u approach, the preconditioning operator
$\mathcal{M}_{2}^{-1}$ can be written as an inexact block LDU factorization of
$\mathcal{J}$. With the permutation matrix $\mathcal{Q}_{2}$:

$\mathcal{Q}_{2}=\begin{bmatrix}0&I_{t}&0\\\ I_{u}&0&0\\\ 0&0&I_{p}\\\
\end{bmatrix},$ (37)

the block LDU factorization reads:

$\mathcal{Q}_{2}\mathcal{J}\mathcal{Q}_{2}^{T}\simeq\mathcal{L}_{2}\mathcal{D}_{2}\mathcal{U}_{2},$
(38)

with:

$\displaystyle\mathcal{L}_{2}$ $\displaystyle=\begin{bmatrix}I_{t}&0&0\\\
-C_{1}\widetilde{H}^{-1}&I_{u}&0\\\ 0&Q_{2}\widetilde{S}_{1}^{-1}&I_{p}\\\
\end{bmatrix},$ $\displaystyle\mathcal{D}_{2}$
$\displaystyle=\begin{bmatrix}-\widetilde{H}&0&0\\\ 0&S_{1}&0\\\
0&0&\widetilde{S}_{2}\\\ \end{bmatrix}$ and $\displaystyle\mathcal{U}_{2}$
$\displaystyle=\begin{bmatrix}I_{t}&-\widetilde{H}^{-1}C_{2}&0\\\
0&I_{u}&\widetilde{S}_{1}^{-1}Q_{1}\\\ 0&0&I_{p}\\\ \end{bmatrix}.$ (39)

The final algebraic expression of $\mathcal{M}_{2}^{-1}$ is therefore:

$\mathcal{M}_{2}^{-1}=\left[\begin{array}[]{ccc}I_{t}&\widetilde{H}^{-1}C_{2}&-\widetilde{H}^{-1}C_{2}\widetilde{S}_{1}^{-1}Q_{1}\\\
0&I_{u}&-\widetilde{S}_{1}^{-1}Q_{1}\\\
0&0&I_{p}\end{array}\right]\left[\begin{array}[]{ccc}-\widetilde{H}^{-1}&0&0\\\
0&\widetilde{S}_{1}^{-1}&0\\\
0&0&\widetilde{S}_{2}^{-1}\end{array}\right]\left[\begin{array}[]{ccc}I_{t}&0&0\\\
C_{1}\widetilde{H}^{-1}&I_{u}&0\\\
-Q_{2}\widetilde{S}_{1}^{-1}C_{1}\widetilde{H}^{-1}&-Q_{2}\widetilde{S}_{1}^{-1}&I_{p}\end{array}\right].$
(40)

######  Remark 4.

Like the t-p-u approach, also $\mathcal{M}_{2}^{-1}$ can be equivalently
recast in other ways. For instance, it can be viewed again as a mixed
constraint preconditioner applied to the matrix:

$\hat{\mathcal{J}}=\left[\begin{array}[]{cc}\mathcal{A}&\mathcal{B}_{1}\\\
\mathcal{B}_{2}&T\end{array}\right],\qquad\mbox{ with
}\;\mathcal{A}=\left[\begin{array}[]{cc}-H&C_{2}\\\
C_{1}&A\end{array}\right],\quad\mathcal{B}_{1}=\left[\begin{array}[]{c}0\\\
Q_{1}\end{array}\right],\quad\mathcal{B}_{2}=\left[\begin{array}[]{cc}0&Q_{2}\end{array}\right],$
(41)

where the inverse of the leading block $\mathcal{A}$ is approximated by an
inner inexact constraint preconditioner [104, 105]. Alternatively,
$\mathcal{M}_{2}^{-1}$ can be regarded as a double application of the
multigrid reduction framework: (i) to matrix $\mathcal{A}$, with
$\widetilde{H}$ instead of $H$, and (ii) to matrix $\mathcal{S}$:

$\mathcal{S}=\left[\begin{array}[]{cc}S_{1}&Q_{1}\\\
Q_{2}&T\end{array}\right].$ (42)

Let us introduce the matrices:

$\displaystyle E_{S_{1}}$ $\displaystyle=I_{u}-\widetilde{S}_{1}^{-1}S_{1},$
(43a) $\displaystyle E_{S_{2}}$
$\displaystyle=I_{p}-\widetilde{S}_{2}^{-1}\hat{S}_{2},$ (43b)

where $\hat{S}_{2}=T-Q_{2}\widetilde{S}_{1}^{-1}Q_{1}$. Along with $E_{H}$
already introduced in equation (23a), $E_{S_{1}}$ and $E_{S_{2}}$ are matrix
measures of the quality of the approximations $\widetilde{S}_{1}$ and
$\widetilde{S}_{2}$ introduced in $\mathcal{M}_{2}^{-1}$. The following result
holds.

######  Proposition 3.2.

The eigenvalues $\lambda_{2}\in\mathbb{C}$ of the preconditioned matrix with
the t-u-p approach are such that:

$|\lambda_{2}-1|\leq(1+\gamma)\max\\{\varepsilon_{H},\varepsilon_{S_{1}},\varepsilon_{S_{2}}\\},$
(44)

with $\varepsilon_{H}=\|E_{H}\|$, $\varepsilon_{S_{1}}=\|E_{S_{1}}\|$,
$\varepsilon_{S_{2}}=\|E_{S_{2}}\|$, $\gamma=\|\mathcal{G}\|$ and

$\mathcal{G}=\left[\begin{array}[]{ccc}C_{2H}\hat{I}_{Q}C_{1S}&-C_{2H}\hat{I}_{Q}&C_{2H}Q_{1S}\\\
\hat{I}_{Q}C_{1S}&-Q_{1S}Q_{2S}&Q_{1S}\\\
-Q_{2S}C_{1S}&Q_{2S}&0\end{array}\right],$ (45)
$C_{1S}=\widetilde{S}_{1}^{-1}C_{1},\quad C_{2H}=\widetilde{H}^{-1}C_{2},\quad
Q_{1S}=\widetilde{S}_{1}^{-1}Q_{1},\quad
Q_{2S}=\widetilde{S}_{2}^{-1}Q_{2},\quad\hat{I}_{Q}=I_{u}+Q_{1S}Q_{2S},$ (46)

for any compatible matrix norm.

###### Proof.

Recalling equations (37) and (40) and introducing the error matrices (23a) and
(43), the preconditioned matrix with the t-u-p approach reads:

$\mathcal{M}_{2}^{-1}\mathcal{Q}_{2}\mathcal{J}\mathcal{Q}_{2}^{T}=\left[\begin{array}[]{ccc}I_{t}+\left[-I_{t}+\widetilde{H}^{-1}C_{2}\left(I_{u}+\widetilde{S}_{1}^{-1}Q_{1}\widetilde{S}_{2}^{-1}Q_{2}\right)\widetilde{S}_{1}^{-1}C_{1}\right]E_{H}&-\widetilde{H}^{-1}C_{2}\left(I_{u}+\widetilde{S}_{1}^{-1}Q_{1}\widetilde{S}_{2}^{-1}Q_{2}\right)E_{S_{1}}&\widetilde{H}^{-1}C_{2}\widetilde{S}_{1}^{-1}Q_{1}E_{S_{2}}\\\
\left(I_{u}+\widetilde{S}_{1}^{-1}Q_{1}\widetilde{S}_{2}^{-1}Q_{2}\right)\widetilde{S}_{1}^{-1}C_{1}E_{H}&I_{u}-\left(I_{u}+\widetilde{S}_{1}^{-1}Q_{1}\widetilde{S}_{2}^{-1}Q_{2}\right)E_{S_{1}}&\widetilde{S}_{1}^{-1}Q_{1}E_{S_{2}}\\\
-\widetilde{S}_{2}^{-1}Q_{2}\widetilde{S}_{1}^{-1}C_{1}E_{H}&\widetilde{S}_{2}^{-1}Q_{2}E_{S_{1}}&I_{p}-E_{S_{2}}\end{array}\right].$
(47)

By making use of the definitions (45) and (46), equation (47) can be re-
written as:

$\mathcal{M}_{2}^{-1}\mathcal{Q}_{2}\mathcal{J}\mathcal{Q}_{2}^{T}=\mathcal{I}+\left(\mathcal{G}-\mathcal{I}\right)\mathcal{E}_{2},$
(48)

with:

$\mathcal{I}=\left[\begin{array}[]{ccc}I_{t}&0&0\\\ 0&I_{u}&0\\\
0&0&I_{p}\end{array}\right],\qquad\mathcal{E}_{2}=\left[\begin{array}[]{ccc}E_{H}&0&0\\\
0&E_{S_{1}}&0\\\ 0&0&E_{S_{2}}\end{array}\right].$ (49)

The eigenvalues $\lambda_{2}$ of the matrix in equation (48) satisfy the bound
(44), thus closing the proof. ∎

######  Remark 5.

The outcome of Proposition 3.2 is very close to that of Proposition 3.1, so
the two approaches are expected to behave similarly. However, with
$\mathcal{M}_{2}^{-1}$ an additional contribution to the error arises,
$\varepsilon_{S_{2}}$, and there is no guarantee that a cluster of eigenvalues
is exactly 1. Moreover, as already observed, the cost for the
$\mathcal{M}_{2}^{-1}$ application is approximately twice that of
$\mathcal{M}_{1}^{-1}$. To decrease this cost of $\mathcal{M}_{2}^{-1}$, one
can apply an incomplete block factorization, i.e., only
$\mathcal{D}_{2}^{-1}\mathcal{L}_{2}^{-1}$. The preconditioner
$\mathcal{M}_{2}^{\star,-1}$ could behave similarly to the original block
preconditioner $\mathcal{M}_{2}^{-1}$, but almost halving the application
cost.

######  Remark 6.

As with the t-p-u approach, the Schur complements $S_{1}$ and
$\widetilde{S}_{2}$, equation (30) and (33), respectively, have algebraic
properties that change with the fracture state and the evolution of the
stick/slip/open partitioning. $S_{1}$ is symmetric in stick mode and non-
symmetric otherwise, but in any case is positive definite. By distinction,
$\widetilde{S}_{2}$ is symmetric anyway, but can be indefinite with open
fracture elements, because $Q_{2}$ depends on the current pressure solution
and no a-priori considerations can be done.

Recalling equation (13), we have that the Schur complement $S$ of the t-p-u
approach is linked to $S_{1}$ by:

$S=S_{1}+\frac{1}{\Delta
t}Q_{1}\widetilde{T}^{-1}\Psi\left(Q_{1}^{T}\right)-Q_{1}\widetilde{T}^{-1}F_{u}.$
(50)

The possible indefiniteness of $S$ is related to the contribution depending on
matrix $F_{u}$. The second Schur complement $S_{2}$ reads:

$S_{2}=T+\frac{1}{\Delta
t}\Psi\left(Q_{1}^{T}\right)\widetilde{S}_{1}^{-1}Q_{1}-F_{u}S_{1}^{-1}Q_{1}.$
(51)

Again, the sum of the first two contributions is positive definite, while the
indefiniteness arises from the term depending on $F_{u}$. There is no general
indication on the actual probability of either $S$ or $S_{2}$ to become
indefinite during a full simulation.

##  4 Numerical results

Three sets of numerical experiments are used to investigate the robustness,
scalability and computational performance of the proposed preconditioning
framework. The first set (Test 1) consists of a small size single-fracture
problem and is used to analyze and compare the robustness of the t-p-u and
t-u-p approaches. The second set (Test 2) simulates the behavior of a number
of uniformly discretized fractures with the aim at investigating the weak
scalability. Finally, we consider a large-size realistic application (Test 3),
representing a tilted well in a hydraulic fracturing stimulation process, in
order to verify the computational efficiency in a meaningful context.

Figure 4: Flow rate in time for all the test cases.

In all test cases, a linearly increasing flow discharge $q_{s}$ is first
injected in the fractures, with a maximum value of $5\cdot 10^{-4}$ m3s-1, and
then extracted, following the time history depicted in Figure 4. The overall
process covers $t_{\max}=6$ s and is discretized into 12 uniform time steps
($\Delta t=0.5$ s). The homogeneous rock material has Young’s modulus and
Poisson’s ratio equal to $3\cdot 10^{3}$ MPa and $0.25$, respectively. The
friction coefficient for the fractures is $0.577$, i.e., the friction angle
$\theta$ is $30^{\circ}$, with the cohesion $c$ set to zero. According to
[88], the initial conductivity value $C_{f,0}$ of equation (2) is equal to
$10\ \text{mD}\cdot\text{m}=9.87\cdot 10^{-15}\ \text{m}^{2}\cdot\text{m}$. On
the corners of the fracture boundary $\partial\Gamma_{i}$ embedded in the 3D
body $\Omega$, a constant zero pressure ($p_{0}=0$ MPa) is imposed.

For each time step, an outer loop (active-set strategy) is coupled with an
inner loop (Newton’s method). Starting from the last converged stick/slip/open
region partitioning, we iterate with Newton’s method until the 2-norm of the
non-linear residual (8) is reduced by a factor $10^{-6}$. At convergence, a
consistency check is carried out to verify the active/inactive status of each
fracture element. If the initial stick/slip/open region partitioning has to be
updated, another inner Newton’s loop is performed. At each Newton’s iteration,
the system (11) is solved by a right-preconditioned full GMRES algorithm
[106], initialized by the zero vector and stopped when the linear residual is
reduced by a factor $10^{-8}$. Hence, for each entire simulation we consider
the total number of iterations needed for the outer active-set strategy,
$N_{\ell}$, the inner Newton’s loop, $N_{N}$, and the preconditioned GMRES
method, $N_{G}$. The average GMRES iteration count for a single linear system
is denoted by $\overline{N}_{G}$, with $N_{G,\min}$ and $N_{G,\max}$ the
minimum and maximum value, respectively, required during the entire
simulation.

###  4.1 Test 1: Robustness

The purpose of this test case is to verify the robustness of the proposed
t-p-u and t-u-p approaches and highlight the main differences. A relatively
small-size problem is set up, with a single 2-m long vertical fracture
completely passing through the 3D domain, which is a box with sizes $10\times
10\times 0.15$ m (Figure 5). The top and bottom surfaces are fixed along the
$z$-direction, while the vertical sides approximately parallel to the fracture
are compressed by a uniform normal load $\sigma_{0}=100\ \text{MPa}$. The out-
of-plane displacement of the two remaining external surfaces of the 3D body is
prevented. The computational grid consists of 2,944 nodes, 2,046 elements and
120 fracture elements, with the following number of unknowns: $n_{u}=8,832$,
$n_{t}=360$, and $n_{p}=120$.

Figure 5: Test 1: Domain configuration, computational grid and injection location. The vertical exaggeration factor is 2. $t$ [s] | 0.5 | 1.0 | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0 | 5.5 | 6.0 | $N_{\ell}$ | $N_{N}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
t-p-u | 2.0 | 5.0 | 4.0 | 4.6 | 7.8 | 9.3 | 5.9 | 5.4 | 4.9 | 4.6 | 4.5 | 2.0 | 58 | 253
t-u-p | 2.0 | 5.0 | 4.0 | 4.0 | 6.5 | 7.4 | 7.2 | 5.0 | 2.0 | 2.0 | 2.0 | 2.0 | 33 | 153
t-p-u | 80.0 | 72.0 | 72.0 | 62.2 | 57.3 | 43.0 | – | – | – | – | – | – |  22∗ |  79∗
t-u-p | 80.0 | 72.0 | 72.0 | 61.8 | 56.1 | 30.9 | 16.0 | 27.7 | 46.2 | 67.6 | 81.0 | 81.0 | 43 | 199
Table 1: Test 1: Average number of linear iterations $\overline{N}_{G}$ for
the t-p-u and t-u-p approach with nested direct solvers (upper rows) and AMG
algorithms (lower rows). The total number of active-set and Newton iterations,
$N_{\ell}$ and $N_{N}$, are also given. * means that $N_{\ell}$ and $N_{N}$
refer to the first six steps only.

We start the analysis by using nested direct methods to apply
$\widetilde{S}^{-1}$ in the t-p-u approach and $\widetilde{S}_{1}^{-1}$ in the
t-u-p approach. The performance obtained with this configuration can be viewed
as the best potential outcome with regard to the average number of linear
iterations $\overline{N}_{G}$. The result is reported in the upper rows of
Table 1. The two approaches behave quite similarly as far as
$\overline{N}_{G}$ is concerned. The very small average number of GMRES
iterations to converge confirms that the approximations introduced in the
construction of $\widetilde{H}$, $\widetilde{S}$ and $\widetilde{S}_{2}$ are
pretty much acceptable. Recall, however, that the t-u-p application cost per
iteration is roughly twice that of t-p-u, so that the first approach appears
to guarantee an overall better computational performance.

This result can significantly change when AMG solvers are introduced instead
of inner direct methods, as it is in practice mandatory for large-size real-
world simulations. The values of $\overline{N}_{G}$ for each time step,
$N_{\ell}$ and $N_{N}$ are provided in the lower rows of Table 1. It can be
noticed that again the reported values are almost the same for the two
approaches in the first six simulation steps. Notice that the significant
increase of $\overline{N}_{G}$ with respect to the use of nested direct
solvers is mainly due to the bad elemental aspect ratio, which is as small as
$\sim\\!10^{-2}$, used to create this test case. Such a distorted grid
negatively affects the conditioning of $\widetilde{S}$ and $S_{1}$, which the
AMG inner preconditioner is not able to fully address. The main difference is
met at the seventh time step, where GMRES convergence cannot be achieved with
the t-p-u approach. Although the t-p-u reduction method is generally more
efficient than t-u-p, it appears to be less robust in some configurations of
the stick/slip/open fracture regions. Similar behaviors have also been
observed with other numerical experiments as well.

The reason for such an outcome stems from the algebraic properties of the
(approximate) Schur complement matrices, $\widetilde{S}$ in the t-p-u approach
and $S_{1}$-$\widetilde{S}_{2}$ in the t-u-p approach, arising at the seventh
time step of the simulation. To this aim, we analyze the eigenvalue
distributions of such matrices (Figure 6). With the current stick/slip/open
region configuration, $\widetilde{S}$ turns out to be non-symmetric and one
eigenvalue with a negative real part arises. Notice also the large ratio
between the maximum and minimum eigenvalue modulus (Table 2). In this
condition, the AMG method used to approximate the application of
$\widetilde{S}^{-1}$ loses its theoretical properties and is no longer
effective. Replacing AMG with another indefinite inner preconditioner for
$\widetilde{S}$, such as an incomplete LU factorization with partial fill-in,
can fix this issue, at the cost of losing the method scalability. By
distinction, with the t-u-p approach we have the theoretical guarantee that
the first-level Schur complement $S_{1}$ is positive definite independently on
the fracture condition. As it can be seen from Figure 6 and Table 2, here AMG
appears to work quite effectively despite the conditioning of $S_{1}$ is the
same as $\widetilde{S}$. The second-level Schur complement $\widetilde{S}_{2}$
is symmetric by construction, but can be indefinite. In this case,
$\widetilde{S}_{2}$ turns out to be still positive definite, but in any case
its inversion by a nested direct solver would ensure an effective
preconditioner behavior.

Figure 6: Test 1: Eigenvalue distribution in the complex plane for $\widetilde{S}$, $S_{1}$ and $\widetilde{S}_{2}$ at the seventh time step. For $\widetilde{S}$ a zoom around zero in arithmetic scale is also reported. | $|\lambda_{\max}|$ | $|\lambda_{\min}|$ | $|\lambda_{\max}|/|\lambda_{\min}|$
---|---|---|---
$\widetilde{S}$ | $2.52\cdot 10^{5}$ | $2.56\cdot 10^{-1}$ | $9.88\cdot 10^{5}$
$S_{1}$ | $2.52\cdot 10^{5}$ | $2.56\cdot 10^{-1}$ | $9.88\cdot 10^{5}$
$\widetilde{S}_{2}$ | $7.74\cdot 10^{-1}$ | $1.01\cdot 10^{-4}$ | $7.66\cdot 10^{3}$
Table 2: Test 1: Maximum and minimum eigenvalue modulus for $\widetilde{S}$,
$S_{1}$ and $\widetilde{S}_{2}$ at the seventh time step.

On summary, we can conclude that the t-p-u approach can be computationally
more efficient than t-u-p. However, it may lack of robustness in large
problems whenever the stick/slip/open region partitioning in a full simulation
gives rise to an indefinite matrix $\widetilde{S}$. The t-u-p approach, though
more expensive, is also much more robust, because it concentrates the source
of the possible numerical issues into $\widetilde{S}_{2}$, which is a
symmetric matrix with a 2D graph connection that can be effectively addressed
by a nested direct solver. For this reason, the t-u-p approach is to be
preferred in a full simulation, where unpredictable fracture configurations
may arise. Therefore, in the next numerical experiments we will focus on t-u-p
approach alone.

###  4.2 Test 2: Weak scalability

To investigate the weak scalability of the proposed algorithm, we first
consider the test case shown in Figure 7 (Test 2a), consisting of a unitary
cube with four vertical fractures. Fluid injection and extraction is
prescribed at the center of each fracture, simulating the action of a
horizontal well. The external faces parallel to the fractures are subjected to
a compressive constant load ($\sigma_{0}=10$ MPa), while the displacement on
the other boundary faces is prevented. The pressure solution on the fracture
surfaces at time $t=3$ s, i.e., at the maximum fluid injection rate, is also
shown in Figure 7.

Figure 7: Test 2a: Domain configuration (left) and pressure solution at time $t=3$ s (right). | coarse | fine
---|---|---
$n_{u}$ | 55,080 | 408,045
$n_{t}$ | 2,700 | 10,800
$n_{p}$ | 900 | 3,600
total | 58,680 | 422,445
Table 3: Test 2a: Problem size for the coarse and fine mesh.

Two uniform grid refinements are considered, with the linear elemental size
varying from $h=0.04\ \text{m}$ to $h=0.02\ \text{m}$. The number of unknowns
for the coarse and fine discretization is listed in Table 3. The overall
performance of the non-linear solution algorithm is reported in Table 4, where
the number of: (i) steps of the active-set algorithm, and (ii) iterations of
the inner Newton’s loop, are reported for every simulation time for both the
coarse and fine discretization. The top panel of Figure 8 graphically
summarizes the same pieces of information.

|  | time [s]
---|---|---
| $\ell$ | 0.5 | 1.0 | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0 | 5.5 | 6.0
| 1 | 2 | 2 | 2 | 2 | 5 | 6 | 6 | 6 | 4 | 2 | 2 | 2
coarse | 2 |  |  |  | 5 | 4 | 5 | 5 | 5 | 2 |  |  |
| 3 |  |  |  |  | 3 |  | 5 |  |  |  |  |
| 1 | 2 | 2 | 2 | 6 | 6 | 5 | 7* | 6* | 6 | 3 | 2 | 2
| 2 |  |  | 5 | 5 | 5 | 6 | 7 | 6 | 5 | 2 |  |
fine | 3 |  |  | 4 | 4 | 4 | 5 | 5 | 5 | 5 |  |  |
| 4 |  |  | 4 | 3 |  | 4 | 5 | 4 | 4 |  |  |
| 5 |  |  |  |  |  |  | 4 | 4 |  |  |  |
Table 4: Test 2a: Newton’s iterations $N_{N}$ for each simulation time and
active-set step $\ell$ for the coarse (top rows) and the fine mesh (bottom
rows). * means that Newton’s scheme does not converge and starting from a
completely closed configuration is required. Figure 8: Test 2a: Non-linear
convergence profiles (a) and number of linear iterations for each Jacobian
system solution.

The simulation is built so that all operating modes are experimented along the
fracture surfaces. At the beginning, all elements are in stick conditions,
then the slip and open regions progressively increase until $t=3$ s, where
almost all fractures are open. Starting from $t=3.5$ s, the elements start to
close again until they return to the initial condition at $t=6$ s. When a
large number of fracture elements belongs to either the slip or the open
region, the non-linear problem is more difficult and requires more active-set
steps and inner Newton’s iterations to converge. Quite intuitively, a larger
number of fracture elements should require more non-linear iterations, as it
can be appreciated in Table 4. The coarse mesh totals $N_{\ell}=20$ and
$N_{N}=75$, while the fine mesh requires $N_{\ell}=35$ and $N_{N}=154$. By
distinction, the linear solver performance is practically unaffected by the
mesh refinement. The bottom panel of Figure 8 and Table 5 show the number of
GMRES iterations $N_{G}$ at every Jacobian system solution and the average
value $\overline{N}_{G}$ for each time step, respectively. Notice that the
linear iteration count tends to increase when the non-linear problem is more
difficult, i.e., around $t=3$ s. This is due to the different space and time
refining indeed, while the space discretization is halved, the time
discretization remains the same between the coarse and fine grids, and the
coupling among the different physical processes involved in the simulation
changes. However, the proposed solution method appears to be fully scalable
with respect to the grid size.

$t$ [s] | 0.5 | 1.0 | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0 | 5.5 | 6.0
---|---|---|---|---|---|---|---|---|---|---|---|---
coarse | 17.0 | 17.0 | 17.0 | 16.8 | 21.2 | 21.3 | 28.8 | 18.3 | 18.0 | 17.0 | 17.0 | 17.0
fine | 18.0 | 17.0 | 15.5 | 16.9 | 18.2 | 25.1 | 32.5 | 20.4 | 15.4 | 18.0 | 17.0 | 17.0
Table 5: Test 2a: Average number of linear iterations $\overline{N}_{G}$ for
the coarse and fine mesh.

The investigated problem couples embedded 2D structures, i.e., the fractures,
with the variables living in 3D domain. Therefore, a progressive grid
refinement changes also the relative size of the blocks appearing in the
Jacobian matrix $\mathcal{J}$, thus potentially modifying the overall problem
conditioning. To analyze the behavior of the proposed algorithm with different
2D-to-3D ratios, i.e., the value of $(n_{t}+n_{p})$ with respect to $n_{u}$, a
second test case is introduced (Test 2b) consisting of a unitary cube with 7
vertical fractures (Figure 9). The model has been regularly refined six times,
with the mesh size varying from $h=0.1\ \text{m}$ to $h=0.0167\ \text{m}$. The
number of unknowns for each refinement level, along with the percentage of 3D
and 2D variables with respect to the total, is listed in Table 6. As with Test
2a, fluid injection and extraction is prescribed at the center of each
fracture simulating the action of a horizontal well. The external faces
parallel to the fractures are subjected to a compressive constant load
($\sigma_{0}=10$ MPa), while the displacement on the other boundary faces is
prevented. Figure 9 also shows the pressure solution at $t=3$ s.

level | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
cells | $10\times 10\times 10$ | $20\times 20\times 20$ | $30\times 30\times 30$ | $40\times 40\times 40$ | $50\times 50\times 50$ | $60\times 60\times 60$
$n_{u}$ | 4,668 | 31,050 | 97,176 | 221,046 | 420,600 | 714,018
$n_{t}$ | 972 | 3,888 | 8,749 | 15,552 | 24,300 | 34,992
$n_{p}$ | 324 | 1,296 | 2,916 | 5,184 | 8,100 | 11,664
total | 5,964 | 36,234 | 108,840 | 241,782 | 453,060 | 760,674
3D | 78.3% | 85.7% | 89.3% | 91.4% | 92.8% | 93.9%
2D | 21.7% | 14.3% | 10.7% | 8.6% | 7.2% | 6.1%
Table 6: Test 2b: Mesh size and percentage of 3D ($n_{u}$) and 2D
($n_{t}+n_{p}$) variables with respect to the total for different refinement
levels.

(a)

(b)

Figure 9: Test 2b: Domain configuration (a) and pressure solution at time
$t=3$ s (b).

Observe that after six refinements the size of the 2D blocks decreases from
21.7% to 6.1% only of the overall size of $\mathcal{J}$. Moreover, even with
the last refinement level, which totals more than 750,000 unknowns, the size
of $\widetilde{S}_{2}$ is around 11,500, thus fully justifying the use of a
nested direct solver. As already observed in Test 2a, it is expected that the
overall number of non-linear iterations, i.e., $N_{\ell}$ and $N_{N}$,
increases as the grid is progressively refined. This is observed in Figure
10(a), which provides the relative variation of $N_{\ell}$ and $N_{N}$ with
respect to the outcome obtained with the coarsest grid. Figure 10(b) provides
the average, maximum and minimum number of GMRES iterations required by the
linear solver. The average value is practically constant around 17 iterations,
with the oscillations between the maximum and minimum iteration count
comprised between 22 and 13. Hence, a very stable behavior of the proposed
algorithm is obtained also changing the relative size of the matrix blocks in
$\mathcal{J}$.

(a)

(b)

Figure 10: Test 2b: active-set and Newton’s iterations relative to those for
the coarsest grid (a) and average, maximum and minimum GMRES iterations (b)
for each refinement level.

###  4.3 Test 3: Computational efficiency

Finally, the performance of the proposed preconditioning framework is verified
in a realistic application. We consider a test case simulating a tilted well
that intersects several fractures. The problem reproduces the situation met in
real-world applications of hydraulic fracturing stimulation. The well
inclination is -15∘ with respect to the horizontal plane. The model analyzes
the dynamics of 9 fractures located along the well in a $5.0\times 1.6\times
2.3$ m3 box. The fractures have the same size, but different relative
positions with respect to the well, i.e., the well does not intersect all of
them at the same location. The domain undergoes a compressive load parallel to
the fractures ($\sigma_{0}=10$ MPa), while displacements are prevented on the
other boundary faces. The pressure at the four corners of every fractures is
set to 0. The model totals 342,642 nodes and 5,184 fracture elements,
corresponding to $n_{u}=1,027,926$, $n_{t}=15,552$, $n_{p}=5,184$, and an
overall system size of 1,048,662 unknowns. The problem turns out to be
particularly challenging because of the grid distortion and the different
stick/slip/open region partitioning simultaneously obtained in each fracture.
An example of the pressure solution obtained at $t=3.5$ s is reported in
Figure 11(a).

(a)

(b)

Figure 11: Test 3: Pressure solution at $t=3.5$ s (a) and GMRES convergence
profiles for all the Newton iterations at the same time (b).

The overall non-linear simulation is very demanding, with several active-set
steps and restarts for the inner Newton’s loop. The performance of the non-
linear solution algorithm is summarized in Table 7, which provides the active-
set steps and Newton’s loop iterations at every simulation time. The irregular
and non-uniform behavior of the different fractures causes convergence
difficulties at almost every time-step, with an overall number of active-set
steps and Newton’s iterations equal to $N_{\ell}=49$ and $N_{N}=238$. The full
simulation, run on a single node of a small cluster equipped with 16 Intel(R)
Xeon(R) Gold 6130 @ 2.10GHz CPU and 6 TB of RAM, requires a total CPU time of
$6472.5$ s.

| time [s]
---|---
$\ell$ | 0.5 | 1.0 | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0 | 5.5 | 6.0
1 | 2 | 2 | 6 | 8 | 8 | 16 | 2 | 5 | 5 | 4 | 3 | 2
2 |  | 5 | 6 | 6 | 6 | 6 | 6 | 6 | 6 | 4 | 2 |
3 |  | 4 | 5 | 5 | 5 | 6 | 5 | 5 | 4 | 2 |  |
4 |  | 5 | 5 | 5 | 5 | 5 | 5 | 5 | 4 | 3 |  |
5 |  |  | 4 | 5 | 5 | 4 | 4 | 4 | 5 |  |  |
6 |  |  |  |  |  |  |  | 4 | 4 |  |  |
$\overline{N}_{G}$ | 25.0 | 22.3 | 23.3 | 25.9 | 29.2 | 26.8 | 22.9 | 21.2 | 20.5 | 22.1 | 21.3 | 21.0
Table 7: Test 3: Newton’s iterations for each simulation time and active-set
step $\ell$. The average number of GMRES iterations $\overline{N}_{G}$ per
time step is also reported.

Despite the challenges posed by the overall simulation, the performance of the
linear solver proves very stable and efficient. The average number of GMRES
iterations is reported in Table 7 for each simulation time. Over the entire
simulation, we have $\overline{N}_{G}=24.2$, with $N_{G,\min}=16$ and
$N_{G,\max}=185$. The latter corresponds to a Newton step where a large number
of elements move from one region to another. As an example, Figure 11(b) shows
all the convergence profiles obtained by a right-preconditioned GMRES
accelerated by the t-u-p approach for the Newton loop at $t=3.5s$. It can be
noticed the great stability of the solver behavior, even if the fracture state
changes significantly. The average CPU time required in this simulation for a
single system solution is $33.4$ s.

##  5 Conclusions

The simulation of frictional contact mechanics with fluid flow in the fracture
network is an important problem in several engineering applications. The
mathematical model can be numerically solved with the aid of a blended finite
element/finite volume formulation, giving rise to a strongly non-linear
problem addressed by an active-set strategy coupled with an inner Newton
iteration. At each Newton step, a linear system with a $3\times 3$ block
Jacobian matrix has to be solved. Since standard global approaches cannot be
effectively used with the resulting non-symmetric and indefinite matrix, this
work focused on the development of a robust, scalable and efficient
preconditioning framework for the solution of the inner linear problem.

The algebraic properties of the $3\times 3$ block Jacobian matrix change
during a full simulation with the evolution of the stick/slip/open region
partitioning of the fractures. In particular, different couplings may arise
and disappear, with the Schur complements possibly being symmetric positive
definite, only positive definite or non-symmetric and indefinite. The proposed
preconditioning framework exploits the physics-based variable partitioning and
the use of multigrid techniques for the sake of algorithmic scalability. The
basic idea relies on restricting the system to a single-physics problem,
approximately solve it by an inner AMG, and prolong the solution back to the
full multi-physics problem. In particular, two multi-physics reduction
sequences are developed, denoted as t-p-u (traction-pressure-displacement) and
t-u-p (traction-displacement-pressure) approaches, and compared in a set of
numerical examples. The results that follow are worth summarizing.

  * 1.

Theoretical analyses show that the proposed approaches are expected to have a
similar convergence rate, with a slight advantage for the t-p-u approach
because of a more clustered eigenvalue distribution for the preconditioned
matrix and a smaller application cost. Indeed, this is confirmed by the
numerical experiments if nested direct solvers are used in the preconditioner
application, but this approach may soon lose robustness when AMG methods are
introduced. The reason stems from the possible indefiniteness of the arising
Schur complement, which is avoided in the t-u-p approach.

  * 2.

The proposed approach proves to be algorithmically scalable with respect to
the computational grid size and the relative size of the discrete fracture
network to the full 3D domain. Although the non-linear problem can become
harder to solve, the iteration count for the inner linear solver is
independent on the discretization size.

  * 3.

The application in a realistic configuration, simulating a hydraulic
fracturing stimulation through a tilted well, shows the computational
efficiency of the proposed approach. Despite the difficulty met by the non-
linear algorithm, due to the combination of variable stick/slip/open operating
modes for the different fractures, the linear solver exhibits a very stable
behavior throughout the full simulation and a remarkable efficiency also in a
sequential implementation.

Future developments regard the implementation of the proposed preconditioning
framework in high performance computing infrastructures, in order to fully
exploit the algorithmic scalability and test the actual parallel efficiency,
and the extension to (multi-phase) fluid flow in the porous matrix as well.

## Acknowledgements

Partial funding was provided by TotalEnergies through the FC-MAELSTROM
project. Portions of this work were performed within the 2020 INdAM-GNCS
project “Optimization and advanced linear algebra for PDE-governed problems”.
Computational resources were provided by University of Padova Strategic
Research Infrastructure Grant 2017: “CAPRI: Calcolo ad Alte Prestazioni per la
Ricerca e l’Innovazione”.

## References

  * [1] S.-Y. Pan, M. Gao, K. J. Shah, J. Zheng, S.-L. Pei, P.-C. Chiang, Establishment of enhanced geothermal energy utilization plans: barriers and strategies, Renew. Energ. 132 (2019) 19–32.  doi:10.1016/j.renene.2018.07.126.
  * [2] X. Wei, Z. J. Feng, Y. S. Zhao, Numerical simulation of thermo-hydro-mechanical coupling effect in mining fault-mode hot dry rock geothermal energy, Renew. Energ. 139 (2019) 120–135.  doi:10.1016/j.renene.2019.02.070.
  * [3] P. Asai, P. Panja, J. McLennan, J. Moore, Efficient workflow for simulation of multifractured enhanced geothermal systems (EGS), Renew. Energ. 131 (2019) 763–777.  doi:10.1016/j.renene.2018.07.074.
  * [4] L. Williams, B. K. Sovacool, The discursive politics of ‘fracking’: Frames, storylines, and the anticipatory contestation of shale gas development in the United Kingdom, Global Environ. Change 58 (2019) 101935.  doi:10.1016/j.gloenvcha.2019.101935.
  * [5] H. Tan, J. Xu, G. Wong-Parodi, The politics of Asian fracking: Public risk perceptions towards shale gas development in China, Energy Res. Social Sci. 54 (2019) 46–55.  doi:10.1016/j.erss.2019.03.007.
  * [6] M. Krzaczek, M. Nitka, J. Kozicki, J. Tejchman, Simulations of hydro-fracking in rock mass at meso-scale using fully coupled DEM/CFD approach, Acta Geotech. 15 (2) (2020) 297–324.  doi:10.1007/s11440-019-00799-6.
  * [7] C. Fan, D. Elsworth, S. Li, L. Zhou, Z. Yang, Y. Song, Thermo-hydro-mechanical-chemical couplings controlling CH4 production and CO2 sequestration in enhanced coalbed methane recovery, Energy 173 (2019) 1054–1077.  doi:10.1016/j.energy.2019.02.126.
  * [8] Q. Li, W. Fei, J. Ma, M. Jing, X. Wei, Coupled CO2 sequestration simulation using ABAQUS and ECLIPSE, Environ. Geotech. 40 (2019) 1–10.  doi:10.1680/jenge.18.00036.
  * [9] P. Liu, T. Zhang, S. Sun, A tutorial review of reactive transport modeling and risk assessment for geologic CO2 sequestration, Comput. Geosci. 127 (2019) 1–11.  doi:10.1016/j.cageo.2019.02.007.
  * [10] P. Zhou, H. Yang, B. Wang, J. Zhuang, Seismological investigations of induced earthquakes near the Hutubi underground gas storage facility, J. Geophys. Res.: Solid Earth 124 (8) (2019) 8753–8770.  doi:10.1029/2019JB017360.
  * [11] V. I. Karev, Geomechanical Approach to Improving the Efficiency of the Operation of Underground Gas Storages, in: Physical and Mathematical Modeling of Earth and Environment Processes (2018), Springer, 2019, pp. 150–158.  doi:10.1007/978-3-030-11533-3_16.
  * [12] P. A. L. P. Firme, D. Roehl, C. Romanel, Salt caverns history and geomechanics towards future natural gas strategic storage in Brazil, J. Nat. Gas Sci. Eng. 72 (2019) 103006.  doi:10.1016/j.jngse.2019.103006.
  * [13] M. Ferronato, G. Gambolati, C. Janna, P. Teatini, Geomechanical issues of anthropogenic CO2 sequestration in exploited gas fields, Energy Conv. Manag. 51 (10) (2010) 1918–1928.  doi:10.1016/j.enconman.2010.02.024.
  * [14] N. Castelletto, G. Gambolati, P. Teatini, Geological CO2 sequestration in multi-compartment reservoirs: Geomechanical challenges, J. Geophys. Res.: Solid Earth 118 (5) (2013) 2417–2428.  doi:10.1002/jgrb.50180.
  * [15] M. D. Zoback, Reservoir geomechanics, Cambridge University Press, 2010.  doi:10.1017/cbo9780511586477.
  * [16] R. E. Goodman, R. L. Taylor, T. L. Brekke, A model for the mechanics of jointed rock, J. Soil Mech. Found. Div. 94 (3) (1968) 637–659.  doi:10.1061/JSFEAQ.0001133.
  * [17] M. Ferronato, G. Gambolati, C. Janna, P. Teatini, Numerical modelling of regional faults in land subsidence prediction above gas/oil reservoirs, Int. J. Numer. Anal. Methods Geomech. 32 (6) (2008) 633–657.  doi:10.1002/nag.640.
  * [18] T. A. Garipov, M. Karimi-Fard, H. A. Tchelepi, Discrete fracture model for coupled flow and geomechanics, Comput. Geosci. 20 (1) (2016) 149–160.  doi:10.1007/s10596-015-9554-z.
  * [19] R. R. Settgast, P. Fu, S. D. C. Walsh, J. A. White, C. Annavarapu, F. J. Ryerson, A fully coupled method for massively parallel simulation of hydraulically driven fractures in 3-dimensions, Int. J. Numer. Anal. Methods Geomech. 41 (5) (2017) 627–653.  doi:10.1002/nag.2557.
  * [20] M. Shakiba, K. Sepehrnoori, Using embedded discrete fracture model (EDFM) and microseismic monitoring data to characterize the complex hydraulic fracture networks, in: SPE annual technical conference and exhibition, Society of Petroleum Engineers, 2015, pp. 1–23.  doi:10.2118/175142-MS.
  * [21] G. Ren, J. Jiang, R. M. Younis, A fully coupled XFEM-EDFM model for multiphase flow and geomechanics in fractured tight gas reservoirs, Procedia Comput. Sci. 80 (2016) 1404–1415.  doi:10.1016/j.procs.2016.05.449.
  * [22] D. L. Y. Wong, F. Doster, S. Geiger, E. Francot, F. Gouth, Investigation of Water Coning Phenomena in a Fractured Reservoir Using the Embedded Discrete Fracture Model (EDFM), in: 81st EAGE Conference and Exhibition 2019, Society of Petroleum Engineers, 2019, pp. 1–5.  doi:10.3997/2214-4609.201901303.
  * [23] K. Wu, W. Yu, J. Miao, Integrating complex fracture modeling and EDFM to optimize well spacing in shale oil reservoirs, in: 53rd US Rock Mechanics/Geomechanics Symposium, 2019.
URL https://onepetro.org/ARMAUSRMS/proceedings-pdf/ARMA19/All-
ARMA19/ARMA-2019-0251/1127420/arma-2019-0251.pdf

  * [24] D. Deb, K. Das, Extended finite element method (XFEM) for analysis of cohesive rock joint, Geotech. Geol. Eng. 28 (5) (2010) 643–659.  doi:10.1007/s10706-010-9323-7.
  * [25] Y. L. Zhang, X. T. Feng, Extended finite element simulation of crack propagation in fractured rock masses, Mater. Res. Innovations 15 (sup1) (2011) s594–s596.  doi:10.1179/143307511X12858957677037.
  * [26] S. Mohammadi, XFEM fracture analysis of composites, John Wiley & Sons, 2012.  doi:10.1002/9781118443378.
  * [27] B. Flemisch, A. Fumagalli, A. Scotti, A review of the XFEM-based approximation of flow in fractured porous media, in: Advances in Discretization Methods, Springer, 2016, pp. 47–76.  doi:10.1007/978-3-319-41246-7_3.
  * [28] S. Berrone, S. Pieraccini, S. Scialò, Flow simulations in porous media with immersed intersecting fractures, J. Comput. Phys. 345 (2017) 768–791.  doi:10.1016/j.jcp.2017.05.049.
  * [29] M. Vahab, N. Khalili, Numerical investigation of the flow regimes through hydraulic fractures using the X-FEM technique, Eng. Fract. Mech. 169 (2017) 146–162.  doi:10.1016/j.engfracmech.2016.11.017.
  * [30] A. R. Khoei, M. Vahab, M. Hirmand, An enriched–FEM technique for numerical simulation of interacting discontinuities in naturally fractured porous media, Comput. Meth. Appl. Mech. Eng. 331 (2018) 197–231.  doi:10.1016/j.cma.2017.11.016.
  * [31] S. Berrone, S. Scialo, F. Vicini, Parallel meshing, discretization, and computation of flow in massive discrete fracture networks, SIAM J. Sci. Comput. 41 (4) (2019) C317–C338.  doi:10.1137/18M1228736.
  * [32] S. Berrone, D. Grappein, S. Pieraccini, S. Scialo, A three-field based optimization formulation for flow simulations in networks of fractures on nonconforming meshes, SIAM J. Sci. Comput. 43 (2) (2021) B381–B404.  doi:10.1137/20M1319188.
  * [33] S. Koric, A. Gupta, Sparse matrix factorization in the implicit finite element method on petascale architecture, Comput. Meth. Appl. Mech. Eng. 302 (2016) 281–292.  doi:10.1016/j.cma.2016.01.011.
  * [34] A. Franceschini, V. A. Paludetto Magri, G. Mazzucco, N. Spiezia, C. Janna, A robust adaptive algebraic multigrid linear solver for structural mechanics, Comput. Meth. Appl. Mech. Eng. 352 (2019) 389–416.  doi:10.1016/j.cma.2019.04.034.
  * [35] P. Hild, Y. Renard, A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics, NUMMATH. 115 (1) (2010) 101–129.  doi:10.1007/s00211-009-0273-z.
  * [36] B. Jha, R. Juanes, Coupled multiphase flow and poromechanics: A computational model of pore pressure effects on fault slip and earthquake triggering, Water Resour. Res. 5 (2014) 3776–3808.  doi:10.1002/2013WR015175.
  * [37] A. Franceschini, M. Ferronato, C. Janna, P. Teatini, A novel Lagrangian approach for the stable numerical simulation of fault and fracture mechanics, J. Comput. Phys. 314 (2016) 503–521.  doi:10.1016/j.jcp.2016.03.032.
  * [38] R. L. Berge, I. Berre, E. Keilegavlen, J. M. Nordbotten, B. Wohlmuth, Finite volume discretization for poroelastic media with fractures modeled by contact mechanics, Int. J. Numer. Meth. Eng. 121 (4) (2020) 644–663.  doi:10.1002/nme.6238.
  * [39] M. Köppel, V. Martin, J. E. Roberts, A stabilized Lagrange multiplier finite-element method for flow in porous media with fractures, GEM 10 (1) (2019) 7.  doi:10.1007/s13137-019-0117-7.
  * [40] A. Franceschini, N. Castelletto, J. A. White, H. A. Tchelepi, Algebraically stabilized Lagrange multiplier method for frictional contact mechanics with hydraulically active fractures, Comput. Meth. in Appl. Mech. Eng. 368 (2020) 113161\.  doi:10.1016/j.cma.2020.113161.
  * [41] B. Wohlmuth, Variationally consistent discretization schemes and numerical algorithms for contact problems, Acta Numer. 20 (2011) 569–734.  doi:10.1017/s0962492911000079.
  * [42] Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003.  doi:10.1137/1.9780898718003.
  * [43] T. A. Davis, Direct methods for sparse linear systems, SIAM, 2006.  doi:10.1137/1.9780898718881.
  * [44] Y. Saad, ILUT: A dual threshold incomplete LU factorization, Numer. Linear Algebra Appl. 1 (4) (1994) 387–402.  doi:10.1002/nla.1680010405.
  * [45] C. J. Lin, J. J. Moré, Incomplete Cholesky factorizations with limited memory, SIAM J. Sci. Comput. 21 (1) (1999) 24–45.  doi:10.1137/S1064827597327334.
  * [46] M. Benzi, Preconditioning techniques for large linear systems: A survey, J. Comput. Phys. 182 (2) (2002) 418–477.  doi:10.1006/jcph.2002.7176.
  * [47] M. Benzi, C. D. Meyer, M. Tůma, A sparse approximate inverse preconditioner for the conjugate gradient method, SIAM J. Sci. Comput. 17 (5) (1996) 1135–1149.  doi:10.1137/S1064827594271421.
  * [48] W.-P. Tang, Toward an effective sparse approximate inverse preconditioner, SIAM J. Matrix Anal. Appl. 20 (4) (1999) 970–986.  doi:10.1137/S0895479897320071.
  * [49] T. Huckle, Factorized sparse approximate inverses for preconditioning, J. Supercomput. 25 (2) (2003) 109–117.  doi:10.1023/A:1023988426844.
  * [50] C. Janna, M. Ferronato, G. Gambolati, A block FSAI-ILU parallel preconditioner for symmetric positive definite linear systems, SIAM J. Sci. Comput. 32 (5) (2010) 2468–2484.  doi:10.1137/090779760.
  * [51] C. Janna, M. Ferronato, Adaptive pattern research for block FSAI preconditioning, SIAM J. Sci. Comput. 33 (6) (2011) 3357–3380.  doi:10.1137/100810368.
  * [52] C. Janna, M. Ferronato, F. Sartoretto, G. Gambolati, FSAIPACK: A software package for high-performance factored sparse approximate inverse preconditioning, ACM Trans. Math. Software (TOMS) 41 (2) (2015) 1–26.  doi:10.1145/2629475.
  * [53] C. Janna, M. Ferronato, G. Gambolati, Enhanced block FSAI preconditioning using domain decomposition techniques, SIAM J. Sci. Comput. 35 (5) (2013) S229–S249.  doi:10.1137/120880860.
  * [54] V. Dolean, P. Jolivet, F. Nataf, An introduction to domain decomposition methods: Algorithms, Theory, and Parallel Implementation, SIAM, 2015.  doi:10.1137/1.9781611974065.
  * [55] S. Zampini, PCBDDC: A class of robust dual-primal methods in PETSc, SIAM J. Sci. Comput. 38 (5) (2016) S282–S306.  doi:10.1137/15M1025785.
  * [56] S. Badia, A. F. Martín, J. Principe, Multilevel balancing domain decomposition at extreme scales, SIAM J. Sci. Comput. 38 (1) (2016) C22–C52.  doi:10.1137/15M1013511.
  * [57] R. Li, Y. Saad, Low-rank correction methods for algebraic domain decomposition preconditioners, SIAM J. Matrix Anal. Appl. 38 (3) (2017) 807–828.  doi:10.1137/16M110486X.
  * [58] S. F. McCormick, J. W. Ruge, Multigrid methods for variational problems, SIAM J. Numer. Anal. 19 (5) (1982) 924–929.  doi:10.1137/0719067.
  * [59] K. Stüben, Algebraic multigrid (AMG): experiences and comparisons, Appl. Math. Comput. 13 (3-4) (1983) 419–451.  doi:10.1016/0096-3003(83)90023-1.
  * [60] A. Brandt, Algebraic multigrid theory: The symmetric case, Appl. Math. Comput. 19 (1-4) (1986) 23–56.  doi:10.1016/0096-3003(86)90095-0.
  * [61] K. Stüben, A review of algebraic multigrid, in: Numerical Analysis: Historical Developments in the 20th Century, Elsevier, 2001, pp. 331–359.  doi:10.1016/B978-0-444-50617-7.50015-X.
  * [62] Y. Notay, Aggregation-based algebraic multigrid for convection-diffusion equations, SIAM J. Sci. Comput. 34 (4) (2012) A2288–A2316.  doi:10.1137/110835347.
  * [63] M. Brezina, R. Falgout, S. MacLachlan, T. Manteuffel, S. McCormick, J. Ruge, Adaptive smoothed aggregation ($\alpha$ SA) multigrid, SIAM Rev. 47 (2) (2005) 317–346.  doi:10.1137/050626272.
  * [64] P. Vaněk, J. Mandel, M. Brezina, Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems, Computing 56 (3) (1996) 179–196.  doi:10.1007/BF02238511.
  * [65] M. Brezina, R. Falgout, S. MacLachlan, T. Manteuffel, S. McCormick, J. Ruge, Adaptive algebraic multigrid, SIAM J. Sci. Comput. 27 (4) (2006) 1261–1286.  doi:10.1137/040614402.
  * [66] A. Brandt, J. Brannick, K. Kahl, I. Livshits, Bootstrap AMG, SIAM J. Sci. Comput. 33 (2) (2011) 612–632.  doi:10.1137/090752973.
  * [67] A. Brandt, J. Brannick, K. Kahl, I. Livshits, Bootstrap Algebraic Multigrid: Status Report, Open Problems, and Outlook, Numerical Mathematics: Theory, Methods and Applications 8 (1) (2015) 112–135.  doi:10.4208/nmtma.2015.w06si.
  * [68] D. Pasetto, M. Ferronato, M. Putti, A reduced order model-based preconditioner for the efficient solution of transient diffusion equations, Int. J. Numer. Meth. Eng. 109 (8) (2017) 1159–1179.  doi:10.1002/nme.5320.
  * [69] P. D’Ambra, S. Filippone, P. S. Vassilevski, BootCMatch: A software package for bootstrap AMG based on graph weighted matching, ACM Trans. Math. Software (TOMS) 44 (4) (2018) 1–25.  doi:10.1145/3190647.
  * [70] P. D’Ambra, P. S. Vassilevski, Improving solve time of aggregation-based adaptive AMG, Numer. Linear Algebra Appl. 26 (6) (2019) e2269.  doi:10.1002/nla.2269.
  * [71] V. A. Paludetto Magri, A. Franceschini, C. Janna, A novel algebraic multigrid approach based on adaptive smoothing and prolongation for ill-conditioned systems, SIAM J. Sci. Comput. 41 (1) (2019) A190–A219.  doi:10.1137/17M1161178.
  * [72] J. Xu, L. Zikatanov, Algebraic multigrid methods, Acta Numer. 26 (2017) 591–721.  doi:10.1017/S0962492917000083.
  * [73] J. W. Ruge, K. Stüben, Algebraic multigrid, in: Multigrid methods, SIAM, 1987, pp. 73–130.  doi:10.1137/1.9781611971057.ch4.
  * [74] T. A. Manteuffel, J. Ruge, B. S. Southworth, Nonsymmetric algebraic multigrid based on local approximate ideal restriction ($\ell$ AIR), SIAM J. Sci. Comput. 40 (6) (2018) A4105–A4130.  doi:10.1137/17M1144350.
  * [75] T. A. Manteuffel, S. Münzenmaier, J. Ruge, B. Southworth, Nonsymmetric reduction-based algebraic multigrid, SIAM J. Sci. Comput. 41 (5) (2019) S242–S268.  doi:10.1137/18M1193761.
  * [76] R. Webster, Stabilisation of AMG solvers for saddle-point stokes problems, Int. J. Numer. Methods Fluids 81 (10) (2016) 640–653.  doi:10.1002/fld.4199.
  * [77] S. C. Brenner, H. Li, L. Y. Sung, Multigrid methods for saddle point problems: Stokes and Lamé systems, NUMMATH. 128 (2) (2014) 193–216.  doi:10.1007/s00211-014-0607-3.
  * [78] L. Chen, Multigrid methods for saddle point systems using constrained smoothers, Comput. Math. Appl. 70 (12) (2015) 2854–2866.  doi:10.1016/j.camwa.2015.09.020.
  * [79] S. C. Brenner, D. S. Oh, L. Y. Sung, Multigrid methods for saddle point problems: Darcy systems, NUMMATH. 138 (2) (2018) 437–471.  doi:10.1007/s00211-017-0911-9.
  * [80] T. A. Wiesner, M. Mayr, A. Popp, M. W. Gee, W. A. Wall, Algebraic multigrid methods for saddle point systems arising from mortar contact formulations, Int. J. Numer. Methods Eng. 122 (15) (2021) 3749–3779.  doi:10.1002/nme.6680.
  * [81] S. C. Brenner, S. Liu, L.-Y. Sung, Multigrid methods for saddle point problems: Optimality systems, J. Comput. Appl. Math. 372 (2020) 112733.  doi:10.1016/j.cam.2020.112733.
  * [82] B. T. Aagaard, M. G. Knepley, C. A. Williams, A domain decomposition approach to implementing fault slip in finite-element models of quasi-static and dynamic crustal deformation, J. Geophys. Res. Solid Earth 118 (6) (2013) 3059–3079.  doi:10.1002/jgrb.50217.
  * [83] A. Franceschini, N. Castelletto, M. Ferronato, Block preconditioning for fault/fracture mechanics saddle-point problems, Comput. Meth. Appl. Mech. Eng. 344 (2019) 376–401.  doi:10.1016/j.cma.2018.09.039.
  * [84] N. Kikuchi, J. T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM, Philadelphia, PA, USA, 1988.  doi:10.1137/1.9781611970845.
  * [85] T. A. Laursen, Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis, Springer-Verlag Berlin Heidelberg, 2003.  doi:10.1007/978-3-662-04864-1.
  * [86] P. Wriggers, Computational Contact Mechanics, 2nd Edition, Springer-Verlag Berlin Heidelberg, 2006.  doi:10.1007/978-3-540-32609-0.
  * [87] P. A. Witherspoon, J. S. Y. Wang, K. Iwai, J. E. Gale, Validity of cubic law for fluid flow in a deformable rock fracture, Water Resour. Res. 16 (6) (1980) 1016–1024.  doi:10.1029/WR016i006p01016.
  * [88] A. Kamenov, D. Zhu, A. D. Hill, J. Zhang, Laboratory measurement of hydraulic fracture conductivities in the Barnett shale, in: SPE Hydraulic Fracturing Technology Conference, Society of Petroleum Engineers, 2013, pp. 216–227.  doi:10.2118/163839-ms.
  * [89] R. Eymard, T. Gallouët, R. Herbin, Finite Volume Methods, in: P. G. Ciarlet, J. L. Lions (Eds.), Handbook of Numerical Analysis, Vol. 7, Elsevier, 2000, pp. 713–1018.  doi:10.1016/S1570-8659(00)07005-8.
  * [90] R. Eymard, T. Gallouët, R. Herbin, A new finite volume scheme for anisotropic diffusion problems on general grids: convergence analysis, C. R. Math. 344 (6) (2007) 403–406.  doi:10.1016/j.crma.2007.01.024.
  * [91] L. Agélas, D. A. Di Pietro, R. Eymard, R. Masson, An abstract analysis framework for nonconforming approximations of diffusion problems on general meshes, Int. J. Finite Vol. 7 (1) (2010) 1–29.
URL https://hal.archives-ouvertes.fr/hal-00318390

  * [92] J. Nocedal, S. Wright, Numerical optimization, Springer Science & Business Media, 2006.  doi:10.1007/b98874.
  * [93] H. Antil, D. P. Kouri, M. D. Lacasse, D. Ridzal, Frontiers in PDE-Constrained Optimization, Springer-Verlag New York, 2018.  doi:10.1007/978-1-4939-8636-1.
  * [94] D. A. May, P. Sanan, K. Rupp, M. G. Knepley, B. F. Smith, Extreme-scale multigrid components within PETSc, in: Proceedings of the Platform for Advanced Scientific Computing Conference, 2016, pp. 1–12.  doi:10.1145/2929908.2929913.
  * [95] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Karpeyev, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, PETSc users manual, Tech. Rep. ANL-95/11 - Revision 3.13, Argonne National Laboratory (2020).
URL https://www.mcs.anl.gov/petsc

  * [96] L. Bergamaschi, M. Ferronato, G. Gambolati, Mixed Constraint Preconditioners for the iterative solution of FE coupled consolidation equations, J. Comput. Phys. 227 (23) (2008) 9885–9897.  doi:10.1016/j.jcp.2008.08.002.
  * [97] M. Ferronato, C. Janna, G. Gambolati, Mixed constraint preconditioning in computational contact mechanics, Comput. Meth. Appl. Mech. Eng. 197 (45) (2008) 3922–3931.  doi:10.1016/j.cma.2008.03.008.
  * [98] M. Ferronato, L. Bergamaschi, G. Gambolati, Performance and robustness of block constraint preconditioners in finite element coupled consolidation problems, Int. J. Numer. Methods Eng. 81 (3) (2010) 381–402.  doi:10.1002/nme.2702.
  * [99] M. Ferronato, A. Franceschini, C. Janna, N. Castelletto, H. A. Tchelepi, A general preconditioning framework for coupled multiphysics problems with application to contact- and poro-mechanics, J. Comput. Phys. 398 (2019) 108887\.  doi:10.1016/j.jcp.2019.108887.
  * [100] Q. M. Bui, L. Wang, D. Osei-Kuffuor, Algebraic multigrid preconditioners for two-phase flow in porous media with phase transitions, Adv. Water Resour. 114 (2018) 19–28.  doi:10.1016/j.advwatres.2018.01.027.
  * [101] Q. M. Bui, D. Osei-Kuffuor, N. Castelletto, J. A. White, A scalable multigrid reduction framework for multiphase poromechanics of heterogeneous media, SIAM J. Sci. Comput. 42 (2) (2020) B379–B396.  doi:10.1137/19M1256117.
  * [102] J. A. White, N. Castelletto, H. A. Tchelepi, Block-partitioned solvers for coupled poromechanics: A unified framework, Comput. Methods Appl. Mech. Eng. 303 (2016) 55–74.  doi:10.1016/j.cma.2016.01.008.
  * [103] N. Castelletto, J. A. White, M. Ferronato, Scalable algorithms for three-field mixed finite element coupled poromechanics, J. Comput. Phys. 327 (2016) 894–918.  doi:10.1016/j.jcp.2016.09.063.
  * [104] L. Bergamaschi, M. Ferronato, G. Gambolati, Novel preconditioners for the iterative solution to FE-discretized coupled consolidation equations, Comp. Meth. Appl. Mech. Eng. 196 (25-28) (2007) 2647–2656.  doi:10.1016/j.cma.2007.01.013.
  * [105] C. Janna, M. Ferronato, G. Gambolati, Parallel inexact constraint preconditioning for ill-conditioned consolidation problems, Comput. Geosci. 16 (3) (2012) 661–675.  doi:10.1007/s10596-012-9276-4.
  * [106] Y. Saad, M. H. Schultz, GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems, SIAM J. Sci. Stat. Comput. 7 (3) (1986) 856–869.  doi:10.1137/0907058.