1 Introduction
In recent years, attention has grown around novel technologies and applications in the subsurface, like geothermal energy production pan2019establishment ; wei2019numerical ; asai2019efficient , hydraulic fracturing williams2019discursive ; tan2019politics ; krzaczek2020simulations , CO sequestration fan2019thermo ; li2019coupled ; liu2019tutorial and underground gas storage zhou2019seismological ; karev2019geomechanical ; firme2019salt . 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 fergamjantea10 ; castelletto2013geological , and, specifically, of faults and fracture networks zoback2010reservoir ; goodman1968model ; ferronato2008numerical ; GarKarTch16 ; Set_etal17 ; shakiba2015using ; ren2016fully ; wong2019investigation ; wu2019integrating ; deb2009extended ; zhang2011extended ; mohammadi2012xfem ; flemisch2016review ; Berrone2017768 ; vahab2017numerical ; khoei2018enriched ; Berrone2019C317 ; Berrone2021B381 . 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 timeconsuming component in a realworld simulation koric2016sparse ; franceschini2019robust .
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) GarKarTch16 , i.e., an explicit representation of the fracture surfaces, while the constraints are imposed with the aid of Lagrange multipliers hild2010stabilized ; JhaJua14 ; FraFerJanTea16 ; berge2020finite ; koppel2019stabilized . As it is common in geological and reservoir simulations, we rely on loworder finite elements for the mechanics and a cellcentered finite volume scheme for the fluid flow. Lagrange multipliers are the contact forces acting on the fracture surfaces as a cellcentered variable, thus sharing the same representation as the fluid pressure field with no interpolation needed. The details of this discretization scheme are described in fr2020alg . This approach is unstable in the LadyzhenskayaBabuškaBrezzi (LBB) sense, i.e., it does not uniformly satisfy the infsup condition (wohlmuth2011variationally, , Section 3.1), and requires a stabilization. In this work, we use the global algebraic approach introduced in the reference work fr2020alg . The Jacobian matrix arising from the described problem is nonsymmetric with a 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 saad2003iterative 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 davis2006direct . However, robustness, scalability and computational efficiency of this class of methods is tightly connected with the choice of a proper preconditioning technique saad2003iterative
. 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
saad1994ilut ; lin1999incomplete ; benzi2002preconditioning , (ii) approximate inverses benzi1996sparse ; tang1999toward ; huckle2003factorized ; janfergam10 ; janfer11 ; janna2015fsaipack , and (iii) multilevel methods, i.e., domain decomposition janfergam13 ; dolean2015introduction ; zampini2016pcbddc ; badia2016multilevel ; li2017low and multigridlike techniques mccormick1982multigrid ; stuben1983algebraic ; brandt1986algebraic ; stuben2001review ; notay2012aggregation ; brezina2005adaptive ; vanvek1996algebraic ; brezina2006adaptive ; brandt2011bootstrap ; brandt2014bootstrap ; Pasetto20171159 ; dambra2018bootcmatch ; dambra2019improving ; paludetto2019novel . 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, xu2017algebraic ) 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., ruge1987algebraic , a wide range of multigrid approaches has appeared in the literature, extending the applicability of this method, originally designed for elliptic PDEs, to both nonsymmetric manteuffel2018nonsymmetric ; manteuffel2019nonsymmetric and block matrices webster2016stabilisation ; brenner2014multigrid ; chen2015multigrid ; brenner2018multigrid ; wiesner2021algebraic ; brenner2020multigrid . Nonetheless, robustness and efficiency is still an open issue for AMG whenever used as a blackbox tool in problems with these algebraic properties. The Jacobian matrix arising from the model considered herein is a nonsymmetric 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 block Jacobian systems aagaard2013domain ; franceschini2019block ; wiesner2021algebraic are found by the authors.The purpose of this work is to design a scalable preconditioning framework for the 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 physicsbased block subdivision and the scalability of AMG techniques available from the literature. The full system is first restricted to a singlephysics 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 realworld 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 realworld benchmarks. A few concluding remarks close the paper.
2 Problem statement
We model the deformation of an open elastic domain , assuming quasistatic conditions and infinitesimal strains within the open time interval . We denote by its boundary, with , and
the outer normal vector to
, while a set of internal boundaries represents a fracture network consisting of surfaces. The external boundary is subdivided into two nonoverlapping subsets, and , where Dirichlet and Neumann boundary conditions apply, respectively. Each fracture consists of two overlapping surfaces, and , with the orientation defined by a unitary vector orthogonal to the fracture plane. By convention, we choose . The pressure field is defined on the union of the twodimensional (2D) domains , with a onedimensional (1D) curve defining the boundary of each fracture and . The curve is subdivided into two nonoverlapping subsets, and , where Dirichlet and Neumann boundary conditions for the pressure field are imposed. The vector denotes the outer normal direction to. The fluid is assumed to be incompressible, and body forces and buoyancy effects are neglected. The projection of the stress tensor
along , , is the traction vector over , with and its normal and tangential component, respectively, with respect to the fracturelocal reference frame. The traction on 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.The strong form of the initial boundary value problem (IBVP) can be stated as follows KikOde88 ; Lau03 ; Wri06 ; fr2020alg : given the fluid discharge , the prescribed boundary displacement and traction , the prescribed fracture boundary pressure and flux , the initial displacement and pressure , find the displacement , the traction , and pressure such that:
in  (1a)  
in  (1b)  
on  (1c)  
respecting the boundary conditions  
on  (1d)  
on  (1e)  
on  (1f)  
on  (1g)  
and initial conditions  
in  (1h)  
in  (1i)  
subject to the constraints over each and for every time in  
(1j)  
(1k)  
In the problem statement, is the Cauchy stress tensor, with C the fourthorder elasticity tensor; is the fluid volumetric flux in the fracture domain according to Darcy’s law witherspoon1980validity —assuming laminar flow—with the fluid pressure gradient, the fluid viscosity (constant), and the isotropic fracture hydraulic conductivity modeled as in GarKarTch16 :
(2) 
with the conductivity related to two irregular surfaces that are in contact kamenov2013laboratory ; denotes the relative displacement across , where and are the normal and tangential components, respectively, and and are the restrictions of on and ; is the limit value provided by the static Coulomb criterion, with and the cohesion and friction angle, respectively. Since we employ a static Coulomb criterion, the tangential velocity in (1k) is replaced with the tangential displacement increment wohlmuth2011variationally with respect to the previously converged timestep.
In our framework, we assume to be fixed with no propagation. The domain is partitioned into three portions, where the following contact conditions occur:

stick on : the fracture is closed () and the traction vector is unknown;

slip on : the fracture is closed in the normal direction ( and is unknown), but a slip displacement between and is allowed for, with ;

open on : the fracture is fully open and a free relative displacement is allowed for, with .
For additional details regarding the governing formulation, we refer the reader to KikOde88 ; Lau03 ; Wri06 ; fr2020alg .
2.1 Discrete weak form
In the solution to the model problem (1), the traction used as a primary variable plays the role of Lagrange multipliers. Denoting with the appropriate inner product of scalar, vector or tensor functions in the spatial domain , we introduce the finitedimensional subspaces , and :
(3a)  
(3b)  
(3c) 
and the discrete approximations of :
(4) 
where, as before, the pedices and denote the components of a vector function along the normal and tangential direction with respect to a fracturelocal reference frame on every . In (4), , , and denote the number of discrete displacement, traction and pressure unknowns. The weak form of (1) reads fr2020alg : find such that
(5a)  
(5b)  
(5c)  
where is with homogeneous conditions along , is the time step size, and is a weighted inner product representing the classical twopoint flux approximation (TPFA) scheme. This is introduced to allow a unified presentation of the coupled finite element/finite volume model EymGalHer00 ; EymGalHer07 ; Age_etal10 . In particular, we have:
(6) 
where and represent the set of edges included in and , respectively; and are the two adjacent cells and ; and is the harmonic average of onesided transmissibility and associated to and . Finally, and collect the boundary conditions to be prescribed on and , respectively. For further details, we refer the reader to fr2020alg .
To solve the problem (5), we transform the variational inequality (5b) into a variational equality. For this purpose, we apply an activeset algorithm, as described in nocedal2006numerical ; antil2018frontiers ; fr2020alg , which allows to identify the subdivision into stick/slip/open regions for every . At a given step of the activeset algorithm, the stick/slip/open regions of each fracture are fixed and the inequality (5b) becomes:
(7) 
with a coefficient needed to ensure the dimensional consistency of the equation. Introducing in (5a), (7) and (5c) the finitedimensional bases of , and yields the following system of nonlinear discrete residual equations:
(8) 
which is solved by a NewtonKrylov method. In (8), the algebraic vectors , and collect the coefficients , and of the discrete displacement, traction and pressure fields in (4) and is the activeset counter. After convergence of the NewtonKrylov method at the th step of the activeset 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 fr2020alg , i.e., firstorder continuous finite elements for displacements and facecentered piecewiseconstant elements for tractions and pressures, as schematically represented in Figure 1. To model the fractures, we use a DFM approach with a conforming mesh GarKarTch16 , hence any 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), with equal to three times the number of 3D finite element nodes, equal to the number of 2D faces discretizing the fracture network, and equal to . Displacements are represented in the global reference system and tractions are represented in a facebased local reference frame. This approach is intrinsically unstable, as it does not fulfill the infsup condition wohlmuth2011variationally . In this work, we use the global algebraic stabilization proposed in fr2020alg , 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:
(9a)  
(9b)  
where and are the stabilizing bilinear forms for the traction and pressure field, respectively. In particular, we have
(10) 
with denoting the jump of a quantity across the generic internal edge and is a positive definite secondorder tensor providing the appropriate scaling. The discrete formulation of is fully provided in fr2020alg . The contribution is computed as the normal projection of with respect to the surface .
At a given activeset iteration , 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 , we have to:
(11)  
The submatrices in the block Jacobian read:
(12a)  
(12b)  
(12c)  
(12d)  
(12e)  
(12f)  
(12g)  
The partial derivatives appearing in (12) are reported in (fr2020alg, , 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 block matrix of equation (11). The global matrix is large, sparse, and nonsymmetric, with properties that change with the evolution of the stick/slip/open regions in the fracture network. A representative evolution of the nonzero pattern of during a full simulation is shown in Figure 2. The features that follow are worth summarizing.

The first block row of includes the contributions arising from the linear momentum balance of the 3D domain . 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, is the classical symmetric positive definite (SPD) elastic stiffness matrix, while and are tall rectangular blocks collecting a surface measure of the fracture elements and transferring tractions and pressures to the 3D body as applied forces.

In the second block row of , varies as the stick/slip/open fracture regions evolve through the activeset algorithm, in both the entry values and the nonzero pattern (Figure 2). If all the fracture elements are in stick mode, we have that , otherwise the frictional law derivatives appear and .

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

The third block row of includes the contributions arising from the fluid mass balance on the fracture network. The coupling between fluid flow and fracture mechanics is controlled by . In particular, when all fracture elements are in stick mode, and is reducible with a symmetric saddlepoint matrix as leading block. Otherwise, contributions from the flux derivative with respect to the displacements appear, i.e., entries depend on the current pressure solution (Figure 2). By distinction with and , there is no simple relationship between and , in both the entry values and the nonzero pattern. Denoting with the matrixtomatrix operator returning a zero row if the corresponding element index belongs to and the original row if the element index belongs to , can be written as:
(13) where collects the contributions from the flux derivatives with respect to the displacements.

is the sum of the standard transmissibility matrix arising from the TPFA discretization in the 2D domain and the stabilization contribution. As such, it is SPD with the 5point stencil of a 2D discrete Laplacian. Moreover, has a block diagonal structure for all nonintersecting fractures. Observe also that traction and pressure fields are always decoupled.
From the observations above, it appears that matrix changes nature with the evolution of the fracture conditions, moving from a reducible matrix with a symmetric saddlepoint leading block to a general nonsymmetric 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.
3 Preconditioning framework
A preconditioner of is a nonsingular operator whose application to a vector resembles as much as possible the action of . The exact application of to some vector provides the vector such that:
(14) 
with , , and natural subvectors of , respectively. The objective is to approximate the solution to the multiphysics system (14) by exploiting the physicsbased variable partitioning. The system is first reduced to a singlephysics problem, and then prolonged back to the full multiphysics space. According to the selected sequence of reductions, different algorithms may arise.
3.1 Method no. 1: tpu 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 and 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 is singular, so a regular surrogate is needed. A block diagonal approximation can be used instead, where each block is the local stabilization matrix computed for each fracture element. Denoting with such a blockdiagonal approximation, we have:
(15a)  
(15b)  
With (15a) and (15b), the first equation of (14) becomes:
(16) 
which is a singlephysics equilibrium equation on the 3D domain where the elimination of fracture tractions and pressures introduces fictitious stiffness contributions. The matrix at the lefthand side of (16) is the Schur complement :
(17) 
Solution to (16) provides , which, introduced into equations (15), yields the final vector . The multiphysics reduction order performed in this case is tractionpressuredisplacement (tpu) and is schematically summarized in Figure 3.
The computation and inversion of in (17) cannot be performed exactly. The Schur complement is explicitly approximated by :
(18) 
where is a diagonal surrogate for . The inverse 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 brandt2014bootstrap ; dambra2018bootcmatch ; dambra2019improving ; paludetto2019novel . In this work, we use an aggregationbased multigrid as the reference AMG operator. Specifically, the application of is approximated by GAMG may2016extreme , the stateoftheart aggregation based multigrid provided by the PETSc package petscuserref .
The construction and application of the resulting preconditioning operator with the tpu approach is summarized in Algorithms 1 and 2. The operator gives a matrix with the diagonal blocks of , while the operator applies the selected AMG preconditioner of to the vector . Since is generally much larger than and , the cost for applying the exact inverse of and is negligible with respect to the AMG algorithm for . Hence, the latter can be roughly assumed as the cost per iteration for the application.
From an algebraic viewpoint, the preconditioning operator arising from the tpu approach can be written as an inexact block LDU factorization of . Using the permutation matrix :
(19) 
where
is the identity matrix in
andthe zero matrix of proper size, the block LDU factorization reads:
(20) 
with:
and  (21) 
Hence, the final algebraic expression of is:
(22) 
Remark 3.1
The multiphysics reduction approach proposed herein can be equivalently recast in other ways as well. Since we use a twofold approximation for , i.e., exact in , and inexact in , can be regarded as a member of the mixed constraint preconditioner class Bergamaschi2008 ; ferjangam08 ; Ferronato2010 . Similarly, the upper and lower block triangular factors in (22) play the role of decoupling operators for the original multiphysics problem and are the outcome of the generalpurpose algebraic procedure defined in Ferronato2019 . Finally, can be also regarded as an example of application in a block nonsymmetric context of the multigrid reduction framework Bui2018 ; Bui2020 , where fracture and body variables play the role of fine and coarse nodes, respectively, and replaces in matrix .
Let us introduce the matrices:
(23a)  
(23b)  
which can be regarded as a matrix measure of the quality of the approximations and introduced in . The following result holds.
Proposition 3.2.
The eigenvalues
of the preconditioned matrix with the tpu approach are either 1, with multiplicity , or such that:(24) 
with , , and
(25) 
for any compatible matrix norm.
Proof.
Recalling equations (19) and (22) and introducing the error matrices (23), the preconditioned matrix with the tpu approach reads:
(26) 
which has unitary eigenvalues. The remaining eigenvalues are those of the matrix obtained by dropping the second block row and column from (26):
(27) 
with the identity matrix of order and :
(28) 
The eigenvalues of satisfy the bound (24), thus closing the proof. ∎
Remark 3.3
Proposition 3.2 shows that the distance of from the identity and the approximation quality of are key factors for the overall performance of . While is fixed, notice, however, that of equation (18) has algebraic properties that change with the fracture state and the evolution of the stick/slip/open regions throughout the activeset algorithm. In stick mode, the contribution is symmetric positive semidefinite and . Hence, is SPD. Also in slip mode the contribution is positive definite and , so remains positive definite, though slightly nonsymmetric. With open elements, however, depends on the current pressure solution and no theoretical considerations can be made in general. In these conditions, is an indefinite nonsymmetric matrix.
3.2 Method no. 2: tup approach
An alternative multiphysics 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:
(29) 
where the matrix:
(30) 
is the firstlevel Schur complement. From a physical viewpoint, is an elasticity matrix with fictitious stiffness contributions arising along the fractures from the traction elimination. Then, a second reduction is needed by computing from (29) and introducing it in the third equation of (14):
(31) 
The matrix at the lefthand side of equation (31) is the secondlevel Schur complement:
(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 multiphysics reduction order is tractiondisplacementpressure (tup).
The computation of can be performed exactly, but its inverse has to be approximated. Since its nature is the same as that of of equation (17), we can effectively use an AMG operator, such as GAMG. We denote with the operator that approximately applies . By distinction, cannot be computed exactly. Recalling the physical interpretation of and , we can approximate the contribution