# Domain decomposition preconditioners for efficient parallel simulations of single-phase flow in three-dimensional fractured porous media with a very large number of fractures

Pierre Jolivet, Michel Kern, Frédéric Nataf, Géraldine Pichot, Daniel Zegarra Vasquez

<figure><img src="image_1.png" /></figure>

To cite this version:

Pierre Jolivet, Michel Kern, Frédéric Nataf, Géraldine Pichot, Daniel Zegarra Vasquez. Domain decomposition preconditioners for efficient parallel simulations of single-phase flow in three-dimensional fractured porous media with a very large number of fractures. 2025. hal-05029676v2

HAL Id: hal-05029676

https://inria.hal.science/hal-05029676v2

Preprint submitted on 11 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Domain decomposition preconditioners for efficient parallel
simulations of single-phase flow in three-dimensional fractured
porous media with a very large number of fractures

Pierre Jolivet<sup>‡</sup>

Michel Kern*†

Frédéric Nataf§*

Géraldine Pichot*†

Daniel Zegarra Vasquez*†

**Abstract**

Three-dimensional fractured-porous media are classically described by the discrete fracture matrix (DFM) model, where the rock matrix remains three-dimensional, while the fracture network is of co-dimension one. This study deals with efficient solutions of the linear system arising in single-phase flow in large DFMs using a mixed-hybrid finite element method. It demonstrates the effectiveness of two-level domain decomposition methods based on the GenEO framework for preconditioning iterative solvers applied to such flow problems. Due to the mixed-dimensional nature of the geometry, a specific decomposition strategy is required. The coarse spaces for these preconditioners are built by locally solving spectral problems. This paper presents simulations with very large DFMs, with up to hundreds of thousands of fractures, achieving a low iteration count and short total computation time. Furthermore, this work confirms the scalability and robustness of these preconditioners.

**Keywords**

single-phase flow, heterogeneous fractured porous media, discrete fracture matrix model, conforming meshes, mixed-hybrid finite element method, domain decomposition methods.

**Introduction**

The efficient simulation of single-phase flow in large-scale fractured porous rocks is a challenging computational problem. This paper builds on previous work by some of the authors [55] and aims to enhance computational performance for such simulations.

In the subsurface, fractures are ubiquitous, appearing as discontinuities in the rock matrix medium in the form of narrow zones. There exist many relevant geometric representations of fractures. In this work, we choose an explicit representation called discrete fracture matrix (DFM) [18, 14, 22], where the fractures are modeled as lower dimensional objects [3, 60, 8]. The rock matrix is a three-dimensional domain and fractures, which, in this work are assumed conductive, are considered as two-dimensional flat surfaces inside a three-dimensional space. The distribution of the fractures follows probability laws derived from field observations [25, 17, 32, 31].

In a DFM representation, Darcy-like averaged equations govern flow in both the rock matrix and the fracture network, with a coupling term that describes exchanges between the two media [60, 63]. This results in mixed-dimensional partial differential equations (PDEs), coupling the two-dimensional fracture flow with three-dimensional surrounding flow in the rock matrix.

The mesh of the fractured porous medium is generated in two stages, we first build a surface simplicial mesh of the fracture network, then we generate the volume mesh constrained by the set

<footnote>‡LIP6, CNRS UMR 7606, Sorbonne Université, 4 place Jussieu, 75005 Paris, France.</footnote>

<footnote>*Inria, 48 rue Barrault, 75013 Paris, France.</footnote>

<footnote>†CERMICS, École Nationale des Ponts et Chaussées, 6 avenue Blaise Pascal, 77420 Champs-sur-Marne, France.</footnote>

<footnote>§LJLL, CNRS UMR 7598, Sorbonne Université, 4 place Jussieu, 75005 Paris, France.</footnote>

<footnote>E-mails: pierre@joliv.et, michel.kern@inria.fr, frederic.nataf@sorbonne-universite.fr, geraldine.pichot@inria.fr, daniel.zegarra-vasquez@inria.fr.</footnote>
of fracture faces (conforming meshes). Because the fractures often intersect at very small angles, the surface and the volume mesh will usually contain low-quality elements [55].

In this work, the set of PDEs is discretized with the mixed-hybrid finite element method (MHFEM) [24, 39, 16, 47] of the lowest order, resulting in a large, sparse, symmetric, and positive definite (SPD) linear system. The matrix of this linear system turns out to be highly ill-conditioned. Two main causes of this ill-conditioning have been identified in [55], namely the presence of low-quality elements [29] and strong hydraulic conductivity contrasts. In fact, condition numbers as high as $2 \cdot 10^{16}$ have been observed [55].

The use of black-box solvers has been investigated in [55]. There, it was shown that direct solvers would not be viable for the configurations investigated in this paper, as their memory requirements would be prohibitive. Iterative methods such as the conjugate gradient (CG) [45], preconditioned with purely algebraic methods, e.g., algebraic multigrid (AMG), were found to display very slow convergence. While there exist high quality AMG implementations, such as BoomerAMG [48], GAMG [38] or MueLu [61], it has been observed that multigrid methods have difficulties with hybridized discretization [23, 42, 72] and also with meshes containing low-quality elements [29]. There exist also ad-hoc methods for the specific case of flow in fractured porous media, such as geometric multigrid (for saddle point problems, which is not the case in this work) [10] or approximate block factorization [9].

In this paper, we focus on another family of preconditioners based on domain decomposition (DD) methods. In the family of DD preconditioners, there exist both one-level preconditioners and two-level (and potentially multi-level) methods [70, 34, 66]. Among one-level preconditioners, with or without overlap, one finds additive Schwarz (ASM), restricted ASM (RAS) and Neumann-Neumann methods. Two-level methods are needed to obtain scalability when the number of subdomains is increased. Without being exhaustive, we can list Composyx [27, 41, 1], BDDC [33, 59], spectral coarse spaces as in [62, 35] as well as HPDDM based on a splitting method for the coarse space [2], here called block_splitting, and HPDDM based on a GenEO coarse space [49, 68, 51, 50], here called GenEO.

Even if DD methods have been studied for single-phase flow in three-dimensional fractured porous media with a low number of fractures [4, 58, 7], to the best of our knowledge, no study has been carried out with large number of fractures (tens to hundreds of thousands).

In this work, we focus on the additive Schwarz method (ASM) and the restricted additive Schwarz (RAS) as one-level methods and HPDDM (both splitting and GenEO) as two-level methods. The main goal of this paper is to assess the performance of these preconditioners through massively parallel simulations, by analyzing the number of iterations, the clock time, scalability with respect to the number of MPI processes (which is taken equal to the number of subdomains), and robustness with respect to the geometry and the heterogeneous properties of the DFM. DD preconditioners require additional information beyond the linear system itself: a decomposition, describing how the unknowns are distributed among the MPI processes. The DD methods used in this work are of the overlapping type. The decomposition starts with a partitioning that can be initiated either from the matrix of the linear system or, more naturally, directly from the geometrical domain mesh. The most widely used partitioning software are METIS [52, 53] and SCOTCH [30, 65, 64]. In this paper, we utilize METIS. From this partitioning, the overlap is built.

On the one hand, the ASM and RAS methods, as well as block_splitting, do not require any geometrical information as the decomposition can be done from the matrix alone. On the other hand, GenEO requires additional information, in the form of local subdomain Neumann matrices and a global/local mapping. Providing these data makes it possible to obtain a theoretical guarantee on the condition number. Each Neumann matrix must be obtained from a local assembly over its subdomain mesh, with one layer of overlap. In the specific case of DFMs, care must be taken to build the overlap both from the porous and the fracture network meshes, see § 3.2.1.

The main contributions of the paper are as follows:
* the first large-scale simulations of fractured porous media, handling up to 697k fractures and 243 million degrees of freedom on 6,825 cores in less than 4 minutes, even in the presence of strong permeability and transmissivity contrasts;

* a procedure for constructing an overlapping decomposition tailored to fractured porous media meshes;

* a comprehensive comparison of the GenEO domain decomposition method with state-of-the-art algebraic preconditioners.

This paper is organized as follows. Section 1 recalls the single-phase flow problem and the associated linear system arising from the mixed-hybrid finite element discretization. Section 2 reviews the theory of domain decomposition preconditioners, with a particular focus on the two-level GenEO coarse space. Section 3 demonstrates that the MHFEM discretization of the single-phase flow problem in fractured porous media is compatible with the theoretical framework of GenEO. We also propose an algorithm to decompose the DFM with overlap in order to construct the Neumann matrices and mappings required by GenEO. Section 4 presents the mesh and linear system characteristics of the large-scale DFMs considered in this study. Section 5 reports the numerical results obtained using algebraic multigrid (AMG) and DD preconditioners on the test cases introduced in Section 4. A performance comparison on a medium-sized test case shows that the most efficient preconditioner is GenEO. This preconditioner is then further evaluated in terms of weak and strong scalability on the same medium-sized test case. Finally, the robustness of GenEO with respect to both the geometry of the fracture network and to permeability/transmissivity contrasts between the rock matrix and the fractures is analyzed across a complete set of DFMs of increasing complexity.

# 1 Single-phase flow problem mathematical setting

## 1.1 Discrete fracture matrix notation

We denote by $\Omega$ the domain representing a permeable fractured porous rock. The fracture network $\Omega^f$ embedded in $\bar{\Omega}$ is modeled as a collection of fractures represented as flat polygons: $\Omega^f = \cup_{\ell \in \mathcal{I}} \gamma_\ell$. The cardinal of $\mathcal{I}$ is the number of fractures and is denoted by $N^f$. The rock matrix domain is defined as $\Omega^m = \Omega \setminus \bar{\Omega}^f$.

The fractures $\gamma_\ell$, $\ell \in \mathcal{I}$, partition the rock matrix domain $\Omega^m$ into connected components $\omega_\alpha$, such that $\Omega^m = \cup_{\alpha \in \mathcal{C}} \omega_\alpha$, where $\mathcal{C}$ is the index set of connected components.

Figure 1 illustrates the notation used here in the case of a DFM with 3 fractures and two connected components.

<figure><img src="image_3.png" /><figcaption>Figure 1: Illustration of a DFM $\Omega$ with three fractures, $\Omega^f = \gamma_1 \cup \gamma_2 \cup \gamma_3$ and two connected components, $\Omega^m = \omega_1 \cup \omega_2$</figcaption></figure>
The fractures may intersect along line segments, we denote by $\Sigma$ the set of such fracture intersections. For each $\sigma \in \Sigma$, let $\mathcal{I}_\sigma$ denote the index set of fractures that contain $\sigma$ as an intersection.

The presentation given here is sufficient for the purpose of this paper, which is mainly numerical. We note that a more precise definition of the fractures and the subdomains is needed for the mathematical analysis, we refer the reader to [54] for more details.

## 1.2 Physical governing equations

We consider a single-phase, incompressible, and time-independent flow problem, as detailed for instance in [60, 20, 21, 13, 14, 7, 15]. The unknowns are:

* the fluid flow velocity in the rock matrix domain $u^m = (u_\alpha^m)_{\alpha \in C}$, obtained from the fluid velocity $u_\alpha^m$ computed in each subdomain $\omega_\alpha$, in m/s;

* the fluid flow velocity in the fracture network $u^f = (u_\ell^f)_{\ell \in \mathcal{I}}$, obtained from the fluid velocity $u_\ell^f$ computed in each fracture $\gamma_\ell$, in m²/s;

* the hydraulic head in the rock matrix $p^m = (p_\alpha^m)_{\alpha \in C}$ obtained from the hydraulic head $p_\alpha^m$ computed in each subdomain $\omega_\alpha$, in m;

* the hydraulic head in the fracture network $p^f = (p_\ell^f)_{\ell \in \mathcal{I}}$ obtained from the hydraulic head $p_\ell^f$ computed in each subdomain $\gamma_\ell$, in m.

In the rock matrix, the hydraulic conductivity is a tensor denoted by $\mathcal{K}^m = (\mathcal{K}_\alpha^m)_{\alpha \in C}$ (units m/s). In the fracture network, the transmissivity is also a tensor denoted by $\mathcal{K}^f = (\mathcal{K}_\ell^f)_{\ell \in \mathcal{I}}$ (units m²/s). We assume that both tensors are uniformly symmetric, and positive definite.

The source term in the rock matrix (resp. in the fracture network) is denoted by $f^m = (f_\alpha^m)_{\alpha \in C}$ (resp. $f^f = (f_\ell^f)_{\ell \in \mathcal{I}}$), whose unit is s⁻¹ (resp. m/s).

We need to look at traces of functions in $\Omega^m$ along a fracture $\gamma_\ell$. The functions can be either continuous or discontinuous across $\gamma_\ell$. Therefore, we have to distinguish between the two sides of the fracture $\gamma_\ell$, which we denote informally by $\pm$.

For a smooth function $p^m = (p_\alpha^m)_{\alpha \in C}$ defined in $\Omega^m$, we denote its traces on the two sides $\pm$ of a fracture $\gamma_\ell$ by $\gamma_\ell^{m, \pm} = p^m|_{\gamma_\ell}$.

For a smooth vector field $u^m = (u_\alpha^m)_{\alpha \in C}$ defined in $\Omega^m$, we denote its normal trace on the side $\pm$ of a fracture $\gamma_\ell$ by $\gamma_\ell^{m, \pm} = u_\alpha^m|_{\gamma_\ell} = u_\ell^m \cdot n^\pm$, where $n^\pm$ (resp. $n^-$) is the unit outward normal to $\gamma_\ell$ from the + side (resp. from the - side). We then define the jump of $u^m$ across the fracture $\gamma_\ell$ as

$$ [[u^m]]|_{\gamma_\ell} = \gamma_\ell^{m,+} u_\alpha^m|_{\gamma_\ell} + \gamma_\ell^{m,-} u_\ell^m|_{\gamma_\ell}. $$

In each porous connected component $\omega_\alpha$, $\alpha \in C$, and in each fracture $\gamma_\ell$, $\ell \in \mathcal{I}$, the flow problem is written as follows:

$$ \nabla \cdot u_\alpha^m = f_\alpha^m \quad \text{in } \omega_\alpha, \qquad \nabla_\tau \cdot u_\ell^f = [[u^m]]|_{\gamma_\ell} + f_\ell^f \quad \text{in } \gamma_\ell, \quad (1.1a) $$

$$ u_\alpha^m = -\mathcal{K}_\alpha^m \nabla p_\alpha^m \quad \text{in } \omega_\alpha, \qquad u_\ell^f = -\mathcal{K}_\ell^f \nabla_\tau p_\ell^f \quad \text{in } \gamma_\ell. \quad (1.1b) $$

Equations (1.1a) express the conservation of mass in the rock matrix and in the fractures. The additional term given by the jump $[[u^m]]|_{\gamma_\ell}$ in the fracture models the exchange between the rock matrix and fracture.

Equations (1.1b) express Darcy's law, both in the porous domain in the fractures, linking the velocity to the gradient of the hydraulic head, through the permeability (in the porous matrix) or the transmissivity (in the fractures).

We next state the transmission conditions that have to be satisfied across the fractures, and at fracture-fracture intersections.

For each intersection $\sigma \in \Sigma$, let us consider the associated intersecting fractures $\gamma_\ell$, $\ell \in \mathcal{I}_\sigma$. For each $\ell \in \mathcal{I}_\sigma$, the fracture $\gamma_\ell$ can have either one side or two sides at $\sigma$. These sides are enumerated using a generic index $j \in \mathcal{I}_{\sigma, \ell}$, where $\mathcal{I}_{\sigma, \ell}$ consists either of one or two elements.
For a fracture $\gamma_\ell$ and a side $j, j \in \mathcal{I}_{\sigma,\ell}$, we denote the trace of the hydraulic head of $p_\ell^\mathcal{f}$ by $\gamma_j^{\mathcal{f},\sigma} p_\ell^\mathcal{f}$ and the normal trace of the vector field $u_\ell^\mathcal{f}$ by $\gamma_{n,j}^{\mathcal{f},\sigma} u_\ell^\mathcal{f}$.

We assume the continuity of the hydraulic heads between the rock matrix and the fractures [8]:

$$ \gamma^{\mathcal{m},+} p_{|\gamma_\ell}^{\mathcal{m}} = \gamma^{\mathcal{m},-} p_{|\gamma_\ell}^{\mathcal{m}} = p_\ell^\mathcal{f}, \forall \ell \in \mathcal{I}. \quad (1.2) $$

Similarly, we assume that, for all $\sigma \subset \Sigma$,

$$ \gamma_j^{\mathcal{f},\sigma} p_\ell^\mathcal{f} = \gamma_{j',\sigma}^{\mathcal{f},\sigma} p_{\ell'}^\mathcal{f}, \text{ for } (\ell, \ell') \in \mathcal{I}_\sigma^2, \text{ and } j \in \mathcal{I}_{\sigma,\ell}, j' \in \mathcal{I}_{\sigma,\ell'}, \quad (1.3a) $$

$$ \sum_{\ell \in \mathcal{I}_{\sigma}} \sum_{j \in \mathcal{I}_{\sigma, \ell}} \gamma_{n,j}^{\mathcal{f},\sigma} u_{\ell}^{\mathcal{f}} = 0. \quad (1.3b) $$

Equation (1.3a) and equation (1.3b) apply at each fracture-fracture intersection $\sigma \subset \Sigma$, the former expressing hydraulic head continuity, and the latter enforcing mass conservation [71]. Here we assume that no flow takes place along the intersections (this is a modeling choice, in contrast to the choice made in [18]).

Finally, boundary conditions must be added to the system (1.1). These conditions, imposed on the boundary of the domain, are prescribed for both the porous and the fracture flow problems, and may be either Dirichlet (imposed hydraulic head) or Neumann (imposed flux). We note that for embedded fractures (with part of their boundary not touching the boundary of the porous domain) we impose a zero-flux condition at the embedded fracture endings (“tips”).

## 1.3 Mixed-hybrid finite-element method (MHFEM) discretization

We discretize the problem (1.1)–(1.3b) with the MHFEM of the lowest order [24]. We recall that MHFEM is a variant of the mixed finite element where the flux continuity between adjacent elements is enforced through Lagrange multipliers. The application of MHFEM to flow in fractured porous media is detailed in [47, 54]. Here, we only give some notations and the resulting linear system.

We first describe the fracture mesh. Let $\mathcal{K}_h$ be a mesh of the fracture network, composed of $\mathcal{N}(\mathcal{K}_h)$ triangles, and let $\mathcal{E}_h$ be the set of edges of $\mathcal{K}_h$, with the number of edges denoted by $\mathcal{N}(\mathcal{E}_h)$. We define $\mathcal{E}_h^{\text{in}}$ to be the set of inner edges, (not intersection and not on the boundary), $\mathcal{E}_h^{\text{inter}}$ to be the set of intersection edges, $\mathcal{E}_h^D$ the set of Dirichlet edges and $\mathcal{E}_h^N$ the set of Neumann edges. We have $\mathcal{E}_h = \mathcal{E}_h^{\text{in}} \cup \mathcal{E}_h^{\text{inter}} \cup \mathcal{E}_h^D \cup \mathcal{E}_h^N$. The number of elements in each of these sets is denoted respectively by $\mathcal{N}(\mathcal{E}_h^{\text{in}})$, $\mathcal{N}(\mathcal{E}_h^{\text{inter}})$, $\mathcal{N}(\mathcal{E}_h^D)$ and $\mathcal{N}(\mathcal{E}_h^N)$. Let $\mathcal{E}_h \setminus \mathcal{E}_h^D$ the set of non-Dirichlet edges of cardinal $\mathcal{N}(\mathcal{E}_h \setminus \mathcal{E}_h^D)$. In the same way, let $\mathcal{T}_h$ be a mesh of the porous domain composed of $\mathcal{N}(\mathcal{T}_h)$ tetrahedral elements. Let $\mathcal{F}_h$ be the set of faces of $\mathcal{T}_h$, whose number is denoted by $\mathcal{N}(\mathcal{F}_h)$. We define $\mathcal{F}_h^{\text{in}}$ the set of inner faces, not on the boundary and not fractures (rock faces), $\mathcal{F}_h^{\text{frac}}$ the set of tetrahedral faces that are fracture faces (triangles, see below), $\mathcal{F}_h^D$ the set of tetrahedral faces that are Dirichlet faces and $\mathcal{F}_h^N$ the set of Neumann faces. We have $\mathcal{F}_h = \mathcal{F}_h^{\text{in}} \cup \mathcal{F}_h^{\text{frac}} \cup \mathcal{F}_h^D \cup \mathcal{F}_h^N$. The number of inner faces is denoted by $\mathcal{N}(\mathcal{F}_h^{\text{in}})$, the number of fracture faces is denoted by $\mathcal{N}(\mathcal{F}_h^{\text{frac}})$ and the number of Dirichlet and Neumann faces by $\mathcal{N}(\mathcal{F}_h^D)$ and by $\mathcal{N}(\mathcal{F}_h^N)$ respectively. The set $\mathcal{F}_h \setminus \mathcal{F}_h^D$ of non-Dirichlet faces is of cardinal $\mathcal{N}(\mathcal{F}_h \setminus \mathcal{F}_h^D)$. The set $\mathcal{F}_h^{\text{in},N}$ of non-fracture and non-Dirichlet faces is of cardinal $\mathcal{N}(\mathcal{F}_h^{\text{in},N})$.

The mixed-hybrid formulation is obtained from the mixed finite element formulation by relaxing the flux-continuity requirement across the faces, or the zero-sum flux at fracture intersection and by enforcing these constraints weakly by means of a Lagrange multiplier on the faces or the tetrahedra and the edges of the triangles. We will use the lowest order Raviart–Thomas (–Nédélec) mixed finite element spaces to approximate the velocity on each element. Their precise definition can be found in most finite element texts, see for example [16, 39]) and it is also recalled in [54]. They will be denoted by $RT_0(K)$ in 2D and by $RTN_0(T)$ in 3D. We also denote by $P_0(X)$ the space of piecewise constant functions on mesh elements $X$, where $X$ stands for tetrahedra, triangles or edges.
We now define the following discontinuous spaces:

$$V_h^{\text{hyb,m}} = \{v_h^m \in L^2(\Omega^m), v_h^m|_T \in \text{RTN}^0(T), \forall T \in \mathcal{T}_h\},$$

$$V_h^{\text{hyb,f}} = \{v_h^f \in L^2(\Omega^f), v_h^f|_K \in \text{RT}^0(K), \forall K \in \mathcal{K}_h\},$$

$$S_h^m = \{s_h^m \in L^2(\Omega^m), s_h^m|_T \in P_0(T), \forall T \in \mathcal{T}_h\},$$

$$S_h^f = \{s_h^f \in L^2(\Omega^f), s_h^f|_K \in P_0(K), \forall K \in \mathcal{K}_h\},$$

$$M_h^m = \{\mu_h^m \mid \mu_h^m|_F \in P_0(F), \forall F \in \mathcal{F}_h \setminus \mathcal{F}_h^{\text{frac}}\},$$

$$M_h^f = \{\mu_h^f \mid \mu_h^f|_E \in P_0(E), \forall E \in \mathcal{E}_h\}.$$

and their product versions:

$$V_h^{\text{hyb}} = V_h^{\text{hyb,m}} \times V_h^{\text{hyb,f}}, \quad S_h = S_h^m \times S_h^f, \quad M_h = M_h^m \times M_h^f, \qquad (1.4)$$

We also need subspaces that include Dirichlet boundary conditions:

$$M_{h,D}^{\mathfrak{m}} = \{\mu_h^{\mathfrak{m}} \in M_h^{\mathfrak{m}}, \mu_h^{\mathfrak{m}}|_F = 0, \forall F \in \mathcal{F}_h^{\mathfrak{D}}\}, \quad M_{h,D}^{\mathfrak{f}} = \{\mu_h^{\mathfrak{f}} \in M_h^{\mathfrak{f}}, \mu_h^{\mathfrak{f}}|_E = 0, \forall E \in \mathcal{E}_h^{\mathfrak{D}}\}.$$

and we let $M_{h,D} = M_{h,D}^{\mathfrak{m}} \times M_{h,D}^{\mathfrak{f}}$.

We may now state the discrete mixed-hybrid approximation of the flow problem (1.1)–(1.2)–(1.3). We look for a triple $(u_h, p_h, \lambda_h) \in V_h^{\text{hyb}} \times S_h \times M_{h,D}$ such that:

$$a(u_h, v_h) - b(v_h, p_h) + c(\lambda_h, v_h) = L_a(v_h), \quad \forall v_h \in V_h^{\text{hyb}}, \qquad (1.5a)$$

$$b(u_h, s_h) = L_b(s_h), \quad \forall s_h \in S_h, \qquad (1.5b)$$

$$c(u_h, \mu_h) = 0, \quad \forall \mu_h \in M_{h,D}, \qquad (1.5c)$$

where the bilinear forms $a$, $b$, and $c$ are defined respectively by

$$a(u_h, v_h) = \sum_{T \in \mathcal{T}_h} \int_T (\mathcal{K}^{\mathfrak{m}})^{-1} u_h^{\mathfrak{m}} \cdot v_h^{\mathfrak{m}} + \sum_{K \in \mathcal{K}_h} \int_K (\mathcal{K}^{\mathfrak{f}})^{-1} u_h^{\mathfrak{f}} \cdot v_h^{\mathfrak{f}}, \quad \forall (u_h, v_h) \in V_h^{\text{hyb}}^2,$$

$$b(v_h, s_h) = \sum_{T \in \mathcal{T}_h} \int_T \operatorname{div} v_h^{\mathfrak{m}}|_T s_h^{\mathfrak{m}}|_T + \sum_{K \in \mathcal{K}_h} \int_K \left( \operatorname{div} v_h^{\mathfrak{f}}|_K - [[v_h^{\mathfrak{m}}]]|_K \right) s_h^{\mathfrak{f}}|_K, \quad \forall (v_h, s_h) \in V_h^{\text{hyb}} \times S_h,$$

$$c(v_h, \mu_h) = \sum_{F \in \mathcal{F}_h^{\text{in},N}} \int_F [[v_h^{\mathfrak{m}}]]|_F \mu_h^{\mathfrak{m}}|_F + \sum_{E \in \mathcal{E}_h} \int_E \sum_{K \in \mathcal{K}_h, E \subset \partial K} (\gamma^f v_h^f)|_E \mu_h^f|_E = 0, \quad \forall (v_h, \mu_h) \in V_h^{\text{hyb}} \times M_{h,D},$$

and $L_a$ and $L_b$ are the linear forms associated with the Dirichlet boundary conditions and the source terms.

It is shown in [54] that the discrete problem in (1.5) has a unique solution, and that this solution is the same as that of a related discrete mixed problem (for the sake of space, we have omitted the mixed formulation and its discretization). Actually, this is a well-known consequence of the fact that the condition $c(u_h, \mu_h) = 0, \forall \mu_h \in M_{h,D}$ expresses the contiguity of the fluxes, either for the rock matrix or for the fracture variables (again, see [54] for details).

We now describe the linear system that is obtained after bases and degrees of freedom are chosen for the spaces defined in (1.4).

* The spaces for the fluxes $V_h^{\text{hyb,m}}$ and $V_h^{\text{hyb,f}}$, and thus $V_h^{\text{hyb}}$, will be provided with the natural bases obtained from the product of the local Raviart-Thomas(-Nédélec) spaces. Their
degrees of freedom are the fluxes of each face of each tetrahedron (resp. each edge of each triangle) for the matrix and the fractures respectively:

$$Q^{\mathfrak{m}} = (q_{T,f}^{\mathfrak{m}})_{(T \in \mathcal{T}_h, f \in \{1,2,3,4\})} \in \mathbb{R}^{4N(\mathcal{T}_h)}, \quad Q^{\mathfrak{f}} = (q_{K,e}^{\mathfrak{f}})_{(K \in \mathcal{K}_h, e \in \{1,2,3\})} \in \mathbb{R}^{3N(\mathcal{K}_h)}.$$

For $T \in \mathcal{T}_h, f \in \{1,2,3,4\}$ (resp. $K \in \mathcal{K}_h, e \in \{1,2,3\}$), we let $(w_{T,f}^{\mathfrak{m}})$ (resp. $w_{K,e}^{\mathfrak{f}}$) be the basis functions for $V_h^{\text{hyb},\mathfrak{m}}$ (resp. $V_h^{\text{hyb},\mathfrak{f}}$). Then

$$u_h^{\mathfrak{m}} = \sum_{\substack{T \in \mathcal{T}_h \\ 1 \le f \le 4}} q_{T,f}^{\mathfrak{m}} w_{T,f}^{\mathfrak{m}}, \quad u_h^{\mathfrak{f}} = \sum_{\substack{K \in \mathcal{K}_h \\ 1 \le e \le 3}} q_{K,e}^{\mathfrak{f}} w_{K,e}^{\mathfrak{f}}.$$

* The spaces for the pressure are piecewise constant functions over each tetrahedron or triangles respectively:

$$\mathbf{P}^{\mathfrak{m}} = (p_{T}^{\mathfrak{m}})_{T \in \mathcal{T}_h} \in \mathbb{R}^{N(\mathcal{T}_h)}, \quad \mathbf{P}^{\mathfrak{f}} = (p_{K}^{\mathfrak{f}})_{K \in \mathcal{K}_h} \in \mathbb{R}^{N(\mathcal{K}_h)}.$$

The discrete functions for the pressure are constant on each tetrahedral or triangular element, so that, for $p_h = (p_h^{\mathfrak{m}}, p_h^{\mathfrak{f}}) \in S_h$:

$$p_h^{\mathfrak{m}} = \sum_{T \in \mathcal{T}_h} p_T^{\mathfrak{m}} \mathbb{1}_T, \quad p_h^{\mathfrak{f}} = \sum_{K \in \mathcal{K}_h} p_K^{\mathfrak{f}} \mathbb{1}_K.$$

* The spaces for the multipliers are piecewise constant on the inner faces of the tetrahedral mesh, and on the edges of the triangular mesh respectively:

$$\mathbf{\Lambda}^{\mathfrak{m}} = (\lambda_F^{\mathfrak{m}})_{F \in \mathcal{F}_h^{\text{in},N}} \in \mathbb{R}^{N(\mathcal{F}_h^{\text{in},N})}, \quad \mathbf{\Lambda}^{\mathfrak{f}} = (\lambda_E^{\mathfrak{f}})_{E \in \mathcal{E}_h \setminus \mathcal{E}_h^D} \in \mathbb{R}^{N(\mathcal{E}_h \setminus \mathcal{E}_h^D)}.$$

The multipliers are constant over the faces or the edges of the mesh, so that

$$\lambda_h^{\mathfrak{m}} = \sum_{F \in \mathcal{F}_h^{\text{in},N}} \lambda_F^{\mathfrak{m}} \mathbb{1}_F, \quad \lambda_h^{\mathfrak{f}} = \sum_{E \in \mathcal{E}_h} \lambda_E^{\mathfrak{f}} \mathbb{1}_E.$$

An illustration of the discrete unknowns with two tetrahedral elements $T_1$ and $T_2$ (rock) and one triangle $K$ (fracture) is displayed on Figure 2.

<figure><img src="image_4.png" /><figcaption>Figure 2: Mixed-hybrid finite element unknowns (Raviart–Thomas of the lowest order) for flow in fractured-porous media.</figcaption></figure>

As usual for hybridized formulations, it is possible to eliminate the flux unknowns $Q^{\mathfrak{m}}$ and $Q^{\mathfrak{f}}$ at the element level, as well as the element pressure $\mathbf{P}^{\mathfrak{m}}$ to obtain a reduced system. Note, however, that the face pressure $\mathbf{P}^{\mathfrak{f}}$ cannot be eliminated, as they are involved in the coupling terms. After reordering the unknowns to group matrix and fracture unknowns together, the reduced linear system for the pressure and multiplier unknowns is as follows:

$$\begin{pmatrix} A^{\mathfrak{m},\mathfrak{m}} & A^{\mathfrak{m},\mathfrak{f}} & 0 \\ (A^{\mathfrak{m},\mathfrak{f}})^T & A^{\mathfrak{f},\mathfrak{f}} + D^{\mathfrak{f}} & -R^{\mathfrak{f}} \\ 0 & -(R^{\mathfrak{f}})^T & M^{\mathfrak{f}} \end{pmatrix} \begin{pmatrix} \mathbf{\Lambda}^{\mathfrak{m}} \\ \mathbf{P}^{\mathfrak{f}} \\ \mathbf{\Lambda}^{\mathfrak{f}} \end{pmatrix} = \begin{pmatrix} \tilde{\mathcal{F}}^{\mathfrak{m}} + \mathcal{Y}^{\mathfrak{m}} \\ \tilde{\mathcal{F}}^{\mathfrak{m},\mathfrak{f}} + V^{\mathfrak{m},\mathfrak{f}} + \mathcal{F}^{\mathfrak{f}} \\ \mathcal{Y}^{\mathfrak{f}} \end{pmatrix}, \quad (1.6)$$
where the matrix can be shown to be SPD. The structure shows the coupling between the upper-left 2 × 2 block that represent the rock matrix unknowns and the lower-right 2 × 2 block that represents the fracture unknowns.

The resulting system 1.6 will be denoted by

$$A\Lambda = F. \qquad (1.7)$$

After solving the global linear system (1.7), the hydraulic heads $P^m$ and the fluxes $Q$ can be computed in a post-processing step, element by element. The objective is to solve the global linear system (1.7) efficiently.

# 2 Domain decomposition preconditioners

As seen in [55], when solving the linear system (1.7) with a large number of fractures $N^f$, iterative linear solvers such as CG preconditioned with algebraic methods tend to struggle in terms of clock time and number of iterations. We therefore investigate the applicability and efficiency of one-level and two-level domain decomposition (DD) preconditioners to solve (1.7).

## 2.1 Domain decomposition framework

We provide here a general overview, highlighting the key components of DD methods. A more comprehensive description is available in [34]. We assume that we are solving a problem coming from a discretized PDE posed on a mesh $T_h$ of a domain $\Omega$, leading to a set of discrete unknowns, called the degrees of freedom (dofs) and resulting in a linear system $Ax = b$. We assume that the matrix $A$ is sparse and SPD. The vector of unknowns $x$ is in one-to-one correspondence with the dofs of the problem, which belong to the set $\mathcal{I}^{\text{dofs}}$.

For the application of DD methods, we partition the domain $\Omega$ in $N$ non-overlapping subdomains $\Omega_i$, $1 \le i \le N$, and we let $\mathcal{T}_{h,i}$ be the corresponding submesh of $\Omega_i$. In practice, such subdomains can be created thanks to specialized partitioning software. From this geometrical partition, we create a partition of the set of dofs $\mathcal{I}^{\text{dofs}} = \cup_{i=1}^N \mathcal{I}_i^{\text{dofs}}$ with a specific procedure detailed in § 3.2.1.

In a second step the subdomains $\Omega_i$, $1 \le i \le N$, are enlarged so as to obtain overlapping subdomains $\Omega_i^\delta$, $1 \le i \le N$. We work here with an overlap obtained by adding the neighboring elements to those in $\Omega_i$, $1 \le i \le N$ (referenced here as $\delta = 1$). For each enlarged $\Omega_i^\delta$, we take the associated dofs, denoted by $\mathcal{I}_i^{\text{dofs},\delta}$. The restriction of a vector $u \in \mathbb{R}^{\mathcal{N}(\mathcal{I}^{\text{dofs}})}$ to a subdomain $\Omega_i^\delta$, $1 \le i \le N$, can be expressed as $R_i u$, where $R_i$ is a rectangular $\mathcal{N}(\mathcal{I}_i^{\text{dofs},\delta}) \times \mathcal{N}(\mathcal{I}^{\text{dofs}})$ Boolean matrix; the extension operator is its transpose $R_i^T$.

Now, we define the matrices $\{\tilde{R}_i\}_{i=1}^N$ that are the same operators as $\{R_i\}_{i=1}^N$ except that coefficients corresponding to the dofs in $\{\mathcal{I}_i^{\text{dofs},\delta} \setminus \mathcal{I}_i^{\text{dofs}}\}_{i=1}^N$ are set to 0. These matrices define a Boolean partition of unity [26, 50]:

$$\sum_{i=1}^{N} \tilde{R}_{i}^{T} R_{i} = I_{\mathcal{N}(\mathcal{I}^{\text{dofs}})},$$

where $I_{\mathcal{N}(\mathcal{I}^{\text{dofs}})}$ is the identity matrix of size $\mathcal{N}(\mathcal{I}^{\text{dofs}})$.

## 2.2 One-level domain decomposition methods

We start with one-level DD methods [34]. In the additive Schwarz method (ASM) each subdomain has its own local problem. These local problems are solved independently, and the solutions are combined to form the preconditioner $M_{\text{ASM}}^{-1}$:

$$M_{\text{ASM}}^{-1} = \sum_{i=1}^{N} R_{i}^{T} A_{i}^{-1} R_{i}, \qquad (2.1)$$
where $A_i$ is defined by $A_i = R_iAR_i^T$.

The restricted additive Schwarz (RAS) preconditioner [26] is a non-symmetric version of the ASM preconditioner where local subdomain solutions are restricted to the non-overlapping parts of the domain when forming the global solution:

$$M_{\text{RAS}}^{-1} = \sum_{i=1}^{N} \tilde{R}_{i}^{T} A_{i}^{-1} R_{i}. \qquad (2.2)$$

It has been observed that RAS improves performance with respect to the number of iterations if we compare it to ASM [36], but at the cost of losing symmetry. Thus, such a preconditioner cannot be applied to CG, but can be used with the generalized minimal residual method (GMRES) [67].

## 2.3 Two-level domain decomposition methods

The efficiency of one-level DD preconditioners is known to deteriorate as the number of subdomains increases [34]. In order to make the DD preconditioner more scalable, a key approach is to help one-level methods by adding a second level containing a coarse space correction. In a two-level method, the solution process is split into two components:

* a fine level, where local solutions are computed on overlapping subdomains thanks to ASM or RAS (see Equations (2.1) and (2.2)), leading to a preconditioner component denoted by $M_{\text{fine}}^{-1}$;

* a coarse space correction, where a global correction ensures faster convergence, especially for low-frequency error components that cannot be easily resolved by local subdomain solvers, leading to a preconditioner component denoted by $M_{\text{coarse}}^{-1}$.

The final two-level DD preconditioner is obtained by combining the fine operator (e.g., ASM (2.1) or RAS (2.2)) and the coarse operator $M_{\text{coarse}}^{-1}$. We consider here two possibilities, deflated or balanced:

$$M_{\text{deflated}}^{-1} = M_{\text{coarse}}^{-1} + M_{\text{fine}}^{-1}(I - A M_{\text{coarse}}^{-1}), \qquad (2.3a)$$

$$M_{\text{balanced}}^{-1} = M_{\text{coarse}}^{-1} + (I - M_{\text{coarse}}^{-1}A)M_{\text{fine}}^{-1}(I - A M_{\text{coarse}}^{-1}). \qquad (2.3b)$$

Other possibilities exist, see [69] for a comprehensive overview. Note that in order to use CG, one needs to use $M_{\text{balanced}}^{-1}$ with ASM (so that the preconditioned system remains symmetric). When the fine level method is RAS, we use $M_{\text{deflated}}^{-1}$ together with GMRES.

We now present two methods for constructing the coarse operator $M_{\text{coarse}}^{-1}$: one referred to hereafter as GenEO, and the other as block_splitting.

### 2.3.1 GenEO

One of the most efficient ways to build the coarse operator $M_{\text{coarse}}^{-1}$ is by using spectral information, as demonstrated in the GenEO framework [49, 68, 51, 50], the Generalized Eigenvalue problems (GEVP) on the Overlap. The GEVP solved by GenEO requires going back to the geometric domain $\Omega$, by partitioning the latter in subdomains $\Omega_i$, $1 \le i \le N$, such that $\Omega = \cup_{i=1}^{N} \Omega_i$, and such that we build an extension, denoted by $\Omega_i^\delta$, where one layer of overlap is added to $\Omega_i$. These subdomains $\Omega_i^\delta$, $1 \le i \le N$, are used to build local Neumann matrices $A_i^\delta$, by assembling the flow problem on $\Omega_i^\delta$, which corresponds to a problem with Neumann boundary conditions. The local GEVP on the $i$th subdomain

$$A_i^\delta y_{i,j} = \rho_{i,j} \tilde{R}_i R_i^T A_i R_i \tilde{R}_i^T y_{i,j} \qquad (2.4)$$

computes the $\mu_i \in \mathbb{N}^*$ smallest eigenpairs $\{y_{i,j}, \rho_{i,j}\}_{j=1}^{\mu_i}$, up to a desired threshold $\nu$ (i.e., $\rho_{i,j} \le \nu$). In PETSc/HPPDM implementation, SLEPc [44] is used to solve the GEVP (2.4) on each subdomain $\Omega_i^\delta$ concurrently. Then the eigenvectors are assembled into a local deflation dense matrix $W_i$ on the $i$th subdomain:

$$W_i = [\tilde{R}_i(R_i)^T y_{i,1} \quad \dots \quad \tilde{R}_i(R_i)^T y_{i,\mu_i}]. \qquad (2.5)$$
The global deflation matrix P is defined as

$$P = [(R_1)^T W_1 \quad \dots \quad (R_N)^T W_N]. \qquad (2.6)$$

Finally, the GenEO coarse operator $M_{\text{coarse}}^{-1}$ is defined thanks to a Galerkin product between A and P:

$$M_{\text{coarse}}^{-1} = P(P^T A P)^{-1} P^T. \qquad (2.7)$$

It is worth noting that the need for Neumann matrices in GenEO prevents the implementation of this method using only the entries of the global matrix A.

When choosing $M_{\text{balanced}}^{-1}$ (see (2.3b)) with ASM as fine level, we have the following result on the condition number of the preconditioned version of the linear system $Ax = b$:

**Theorem 2.1** ([34, Theorem 7.23]). The GenEO coarse space guarantees that the condition number $\kappa$ of the two-level preconditioned operator $M_{\text{balanced}}^{-1}A$, with $M_{\text{fine}}^{-1} = M_{\text{ASM}}^{-1}$, is bounded from above:

$$\kappa(M_{\text{balanced}}^{-1}A) \le k_0(1 + k_1\nu^{-1}), \qquad (2.8)$$

where $k_0 \in \mathbb{N}^*$ is the maximum number of neighbors of a subdomain including itself, $k_1 \in \mathbb{N}^*$ is the maximum multiplicity of a degree of freedom and $\nu \in ]0, 1[$ is the user-chosen parameter for the GEVP (2.4).

**Remark 2.2.** The parameter $\nu$ lets the user control the trade-off between the cost and quality of the preconditioner. The inequality 2.8 shows that a small value of $\nu$ means that few eigenvalues will be computed and the condition number may increase, while a larger value of $\nu$ means that more eigenvalues will be computed, likely increasing the quality of the preconditioner. In the experiments reported in Section 5, the value of $\nu$ is set to 0.1.

**Remark 2.3.** The proof of Theorem 2.1 is fully algebraic except for [34, Lemma 7.13], which depends on the existence of a variational formulation of the discrete problem with a bilinear form that is additive and local. This is the only condition that needs to be checked before applying GenEO in our specific case. See Section 3.2.3 for the application to flow in fractured-porous media.

**Remark 2.4.** In practice, one may prefer RAS at the fine level over ASM, as it tends to result in smaller iteration counts, albeit without theoretical control.

### 2.3.2 block_splitting

Another two-level DD preconditioner, called block_splitting, exists which, unlike GenEO, does not require the Neumann matrices $A_i^\delta$. block_splitting is a local block splitting method using lumping in the overlap, as described in detail in [2]. The advantage of this method is that it creates a coarse space, thus helping the fine-level part of the preconditioner, as GenEO would. The main difference is that the left-hand side operator replacing $A_i^\delta$ in equation (2.4) does not necessarily require any geometrical information to construct it. Thus, an algebraic partitioning of A is enough to obtain the two-level preconditioner block_splitting, following similar steps to those in GenEO. However, the construction of such an algebraic operator does not guarantee a sharp bound on the condition number of the preconditioned linear system, unlike GenEO with Neumann matrices, see Theorem 2.1, unless A is diagonally dominant [2].

## 3 Application to single-phase flow in fractured porous media

In this section, we go back to the linear system (1.7) arising in our application. The matrix A is very large, sparse, SPD and, as shown in [55], very ill-conditioned. We aim at solving the system with a Krylov method preconditioned thanks to DD methods.
## 3.1 ASM, RAS, and block_splitting preconditioners

To build the one-level ASM and RAS preconditioners, as well as the two-level block_splitting, only a decomposition of the matrix A is required. To generate a partition of A, one may use a classical partitioning software, like METIS [52, 53] or SCOTCH [65, 64]. Therefore, the application of ASM, RAS, and block_splitting to our application is rather straightforward.

## 3.2 GenEO preconditioner

The use of GenEO in the context of fractured porous flow is not straightforward and there are two specific steps to be performed. First, GenEO needs Neumann matrices for the GEVP (2.4). The main difficulty is not the assembly phase itself, since it is very close to what is done for the global problem. The key difficulty is to create a coherent decomposition of the DFM with overlap. Indeed, in our case, $\Omega$ is made of two different components: $\Omega^f$, meshed conformingly with triangles $K_h$, and $\Omega^m$, meshed conformingly with tetrahedra $T_h$ taking $K_h$ as constraints. One of the main contributions of this work is to give an algorithm in Section 3.2.1 that describes subdomain generation with a one-layer overlap for three-dimensional DFM. We show that, even if the triangles in $K_h$ are essentially “special” faces from the tetrahedra in $T_h$, there are particular configurations we have to consider in order to successfully obtain subdomains with one layer of overlap (see Figure 4). Second, it is required to check that the condition explained in Remark 2.3 is fulfilled. This is done in Section 3.2.3.

### 3.2.1 DFM decomposition

The general idea for the decomposition of a 2D-3D DFM is to first build a partition of the domain into $N$ subdomains, $\Omega = \cup_{i=1}^N \Omega_i$. Each subdomain is defined as $\Omega_i = \Omega_i^m \cup \Omega_i^f$, $1 \le i \le N$, with $\Omega_i^m$ (resp. $\Omega_i^f$) the subdomain $i$ corresponding to the rock matrix medium (resp. fractured medium). Then, we build an overlap for each subdomain $\Omega_i$, denoted by $\Omega_i^\delta = \Omega_i^{m,\delta} \cup \Omega_i^{f,\delta}$, $1 \le i \le N$, with $\Omega_i^{m,\delta}$ (resp. $\Omega_i^{f,\delta}$) the subdomain $\Omega_i^m$ (resp. $\Omega_i^f$) with overlap. Here, we propose to build the subdomains and their overlap from the DFM mesh.

We denote by $\mathcal{T}_{h,i}$ the tetrahedral mesh composing the rock matrix subdomain $\Omega_i^m = \cup_{T \in \mathcal{T}_{h,i}} T$. We denote by $\mathcal{T}_{h,i}^\delta$ the tetrahedral mesh composing the overlap of the subdomain $\Omega_i^m$. Then, each subdomain with overlap $\Omega_i^{m,\delta}$, $1 \le i \le N$, writes $\Omega_i^{m,\delta} = \cup_{T \in (\mathcal{T}_{h,i} \cup \mathcal{T}_{h,i}^\delta)} T$.

Similarly, we denote by $K_{h,i}$ triangular mesh composing the fractures subdomain $\Omega_i^f = \cup_{K \in \mathcal{K}_{h,i}} K$. We denote by $\mathcal{K}_{h,i}^\delta$ the triangular mesh composing the overlap of the subdomain $\Omega_i^f$. Then, each subdomain with overlap $\Omega_i^{f,\delta}$ writes $\Omega_i^{f,\delta} = \cup_{K \in (\mathcal{K}_{h,i} \cup \mathcal{K}_{h,i}^\delta)} K$.

Algorithm 1 describes a workflow that decomposes the DFM domain into overlapping subdomains. Figure 3 displays an example of a fractured porous medium decomposed in 4 subdomains.

<caption>**Algorithm 1** Decomposition with overlap of the DFM domain based on the mesh.</caption>

1: [Find the elements in $\mathcal{T}_{h,i}$] Partition the volume mesh $\mathcal{T}_h$ with METIS to build $\mathcal{T}_{h,i}$, $1 \le i \le N$.

2: **for** $i \leftarrow 1$ to $N$ **do**

3: [Find the elements in $\mathcal{T}_{h,i}^\delta$] Add in $\mathcal{T}_{h,i}^\delta$ the tetrahedral elements not in $\mathcal{T}_{h,i}$ but sharing a face with elements in $\mathcal{T}_{h,i}$.

4: [Find the elements in $K_{h,i}$] Add in $K_{h,i}$ all the fracture faces (triangles) in $\mathcal{T}_{h,i}$.

5: [Find the elements in $K_{h,i}^\delta$] Add in the overlap $K_{h,i}^\delta$ all the triangles that are fracture faces in $\mathcal{T}_{h,i}^\delta$. Also add in $K_{h,i}^\delta$ all the triangles not already in $K_{h,i}^\delta$ but sharing an edge with triangles in $K_{h,i}$.

6: **end for**

Let us explain in more details the different stages of Algorithm 1.

Stage 1 is standard, we use METIS to build the partition $\mathcal{T}_{h,i}$, $1 \le i \le N$.
<figure><img src="image_5.png" /><figcaption>Figure 3: Domain decomposition in 4 subdomains of a fractured porous medium containing 55 fractures. The tetrahedral elements in the overlap $\mathcal{T}_{h,i}^{\delta}$ of the porous mesh are colored in red, the triangles in the overlap $\mathcal{K}_{h,i}^{\delta}$ of the fractures mesh in pink. We associated each subdomain $i$, $i = 1, \dots, 4$, with a given color for clarity.</figcaption></figure>

Then, we perform a loop over the different subdomains $i$, $1 \le i \le N$, to identify the elements in $\mathcal{T}_{h,i}^{\delta}$, $\mathcal{K}_{h,i}$, and $\mathcal{K}_{h,i}^{\delta}$. Within the loop, Stage 3 and Stage 4 are straightforward. Stage 5 finds the triangles in $\mathcal{K}_{h,i}^{\delta}$ in two steps. These two steps are illustrated on Figure 4 with the notation $\mathcal{F}(\mathcal{T}_h) = \{F \in \mathcal{F}_h, F \text{ is a face of } T, T \in \mathcal{T}_h\}$. First, it collects the fractures faces in the overlap $\mathcal{T}_{h,i}^{\delta}$ (Figure 4a). Notice there might be triangles belonging to fractures that end in the overlap (Figure 4b). Second, it collects the triangles that cannot be reached in $\mathcal{T}_{h,i}^{\delta}$ but are reachable through the boundary edges of $\mathcal{K}_{h,i}$ (Figure 4c).

<figure><img src="image_6.png" /><figcaption>(a) $K_2 \in \mathcal{F}(\mathcal{T}_{h,1}^{\delta})$ with a neighbor $K_1 \in \mathcal{K}_{h,1}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) $K_1 \in \mathcal{F}(\mathcal{T}_{h,1}^{\delta})$ with no neighbor in $\mathcal{K}_{h,1}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) $K_2 \notin \mathcal{F}(\mathcal{T}_{h,1}^{\delta})$.</figcaption></figure>

<caption>Figure 4: Sketch of a mesh to display the different configurations of triangles in the overlap $\mathcal{K}_{h,1}^{\delta}$.</caption>

Overall, our decomposition of the fractures is guided by that of the porous matrix. Alternative decomposition strategies could be explored within the HPDDM framework to achieve a better balance between subdomains. This could be done by constraining the graph partitioner to produce subdomains that are as homogeneous as possible while balancing the number of unknowns. Such an approach was applied to porous media in [43], using a cost function inspired by work on algebraic multigrid methods. Extending this strategy to DFM models would require accounting for the fractures and is beyond the scope of this article.

### 3.2.2 Generation of Neumann matrices and mappings

In addition to the global matrix $A$, GenEO requires the following structures: the local Neumann matrices $A_i^{\delta}$, $1 \le i \le N$, the matrices $R_i$ and $\tilde{R}_i$, $1 \le i \le N$.

The Neumann matrix $A_i^{\delta}$, $1 \le i \le N$, is obtained by assembling the problem on the local
subdomain $\Omega_i^\delta$ ($\Omega_i$ with one layer of overlap). As mentioned in Section 2, $A_i^\delta$ corresponds to a discretization of the problem endowed with Neumann boundary conditions.

From the decomposition presented in Section 3.2.1, one can build maps from global unknowns to unknowns local to each subdomain, which are essentially data structures encoding $R_i$ and $\tilde{R}_i$.

### 3.2.3 Compatibility of the discrete bilinear form

As noted in Remark 2.3, the only condition that needs to be checked in order to apply Theorem 2.1 to flow in fractured-porous media is the existence of a local and additive bilinear form. In the MHFE context of this work, it has been shown in [54] that the matrix in the linear system (1.7) does come from a mesh-dependent bilinear form. We recall its expression:

$$a_{\Omega}(\Lambda, \hat{\Lambda}) = \int_{\Omega^{\mathfrak{m}}} (\mathcal{K}^{\mathfrak{m}})^{-1} u_h^{\mathfrak{m}}(\Lambda) \cdot u_h^{\mathfrak{m}}(\hat{\Lambda}) + \int_{\Omega^{\mathfrak{f}}} (\mathcal{K}^{\mathfrak{f}})^{-1} u_h^{\mathfrak{f}}(\Lambda) \cdot u_h^{\mathfrak{f}}(\hat{\Lambda}), \quad (3.1)$$

where the mapping $\Lambda \to u_h(\Lambda)$ is defined through the solution of local version of the MHFE problem. We note that, on each element, $u_h(\Lambda)$ only depends on the values of $\Lambda$ on the same element, see [54] for details.

Since all mappings involved are local, we may split the discrete bilinear form (in its integral definition (3.1)) as a sum over elements, with the contributions to each term coming only from data that is local to the element. We define the subdomain version of the bilinear form as:

$$a_{\Omega_i^\delta}(\Lambda, \hat{\Lambda}) = \int_{\Omega_i^{\mathfrak{m},\delta}} (\mathcal{K}^{\mathfrak{m}})^{-1} u_h^{\mathfrak{m}}(\Lambda) \cdot u_h^{\mathfrak{m}}(\hat{\Lambda}) + \int_{\Omega_i^{\mathfrak{f},\delta}} (\mathcal{K}^{\mathfrak{f}})^{-1} u_h^{\mathfrak{f}}(\Lambda) \cdot u_h^{\mathfrak{f}}(\hat{\Lambda}).$$

It is now sufficient to follow the same steps as in the proof of [34, Lemma 7.13] to obtain the following inequality:

$$\sum_{i=1}^{N} a_{\Omega_i^\delta}(\Lambda, \Lambda) \leq k_1 a_{\Omega}(\Lambda, \Lambda),$$

where $k_1 \in \mathbb{N}^*$ is the maximum multiplicity of the subdomain intersection. With this inequality, we have verified that the discretization with the MHFEM of the single-phase flow problem in fractured-porous media is compliant with the theoretical framework of GenEO.

## 3.3 Software implementation

Flow simulations presented in this work are performed with a combination of different software. Work is in progress to merge all three following steps.

**First**, the assembly phase, and the possible construction of Neumann matrices, are performed with the MATLAB code nef-flow-fpm that implements the mixed-hybrid finite element method, as described in [54]. In particular, what is new with respect to [55] is that nef-flow-fpm is able to decompose DFMs and generate Neumann matrices following the methods described in Section 3.2.1. We recall that nef-flow-fpm is able to store the linear system in PETSc [12, 11] binary format.

**Second**, the linear system (1.7) is solved with PETSc built-in solvers. The PETSc implementation of the DD methods studied in this work has been thoroughly explained in [50].

**Third**, the solution obtained with PETSc is post-processed with nef-flow-fpm, to compute the velocity fields and perform various checks such as local and global mass conservation.
# 4 Discrete fracture matrix test cases

In this work, we follow the construction of test cases introduced in [55], where the rock matrix is modeled as a cube and sub-cubes of varying sizes are extracted to create a set of test cases with increasing complexity. Hereafter, the geometries are labeled according to the cube's side length, which takes the values 60, 80, 100, 130, and 160. The test case with the cube of side length 60 is the smallest considered in this work, containing approximately 39k fractures, while the largest one, with a size length of 160, includes up to 697k fractures. Compared to the naming convention in [55], we omit the "FPM" (fractured porous media) label, since all simulations presented here are exclusively within fractured porous media. In this section, we describe the mesh properties and the characteristics of the linear systems (1.7) that arise from the MHFEM discretization.

## 4.1 Mesh properties

The mesh of the fractured porous medium, modeled as a DFM, is generated in two consecutive stages described in [55].

1. The surface mesh for the DFN is generated with MODFRAC [19, 56]. As in [54], we generate two types of meshes according to the mesh size: "geo" means that the mesh step is automatically set by MODFRAC according to the geometry and "fine" means that we manually chose a finer mesh step.

2. The volume mesh is generated with GHS3D [40], using the DFN mesh as input so that the surface and volume meshes are conforming.

We use the same mesh quality criteria $\delta(K)$ for triangles and $\delta(T)$ for tetrahedra as in [55]. The quality threshold for MODFRAC is chosen equal to $10^{-4}$, and we do not impose a quality threshold for GHS3D.

Table 1 gives the geometry and mesh properties presented geometries in this work with MODFRAC and GHS3D. Table 2 gives the mesh quality for the above examples with MODFRAC and GHS3D.

<table><caption>Table 1: Geometry and mesh properties of the DFMs: cube size, number of fractures, number of fracture-fracture intersections, mesh name, number of tetrahedra, number of faces (total and fracture), number of triangles (fracture) and number of edges.</caption><thead><tr><th>Geometry</th><th>N<sup>g</sup></th><th># fracture-fracture intersections</th><th>Mesh name</th><th>N(T<sub>h</sub>)</th><th>N(F<sub>h</sub>)</th><th>N(K<sub>h</sub>)</th><th>N(E<sub>h</sub>)</th></tr></thead><tbody><tr><td>[L60]</td><td>39,456</td><td>66,181</td><td>[L60geo]</td><td>7,876,897</td><td>15,825,076</td><td>1,790,168</td><td>2,751,907</td></tr><tr><td>[L60]</td><td>39,456</td><td>66,181</td><td>[L60fineA]</td><td>16,197,529</td><td>32,501,063</td><td>2,737,569</td><td>4,236,628</td></tr><tr><td>[L60]</td><td>39,456</td><td>66,181</td><td>[L60fineB]</td><td>29,827,028</td><td>59,807,738</td><td>4,024,218</td><td>6,214,235</td></tr><tr><td>[L60]</td><td>39,456</td><td>66,181</td><td>[L60fineC]</td><td>59,549,925</td><td>119,337,725</td><td>6,321,429</td><td>9,719,152</td></tr><tr><td>[L80]</td><td>90,447</td><td>155,362</td><td>[L80geo]</td><td>16,544,945</td><td>33,206,750</td><td>3,825,145</td><td>5,860,352</td></tr><tr><td>[L100]</td><td>174,315</td><td>305,455</td><td>[L100geo]</td><td>30,454,172</td><td>61,082,397</td><td>7,058,863</td><td>10,801,781</td></tr><tr><td>[L130]</td><td>377,269</td><td>674,705</td><td>[L130geo]</td><td>59,869,820</td><td>120,004,957</td><td>14,122,726</td><td>21,603,994</td></tr><tr><td>[L160]</td><td>696,839</td><td>1,260,231</td><td>[L160geo]</td><td>102,593,939</td><td>205,561,258</td><td>24,385,046</td><td>37,327,575</td></tr></tbody></table>

Table 2 shows that the surface mesher MODFRAC does not generate any elements with quality below the user-specified minimum of $10^{-4}$. As the number of fractures increases, the mean quality values $\delta(K)$ and $\delta(T)$ of the "geo" meshes exhibit a slight decrease and the number of tetrahedral elements with a quality $\delta(T)$ below 0.1 slightly increases. The minimum 3D element quality can be as low as $2.56 \cdot 10^{-9}$, as observed for [L160geo]. Nonetheless, low-quality elements constitute only a small fraction, less than 2%, of the total number of tetrahedral elements. Despite these low-quality elements, we demonstrate in this work that such meshes remain suitable for performing flow simulations.
<table><caption>Table 2: Quality of the 2D cells $K \in K_h$ and 3D cells $T \in T_h$ of the DFM meshes.</caption><thead><tr><th>Mesh name</th><th>mean δ(K)</th><th>min<sub>K</sub> δ(K)</th><th>min<sub>K</sub> |K|</th><th>mean δ(T)</th><th>min<sub>T</sub> δ(T)</th><th>min<sub>T</sub> |T|</th><th>median δ(T)</th><th>N(T<sub>h</sub>) s.t. δ(T) &lt; 0.1</th><th>% N(T<sub>h</sub>) s.t. δ(T) &lt; 0.1</th></tr></thead><tbody><tr><td>[L60geo]</td><td>0.85</td><td>1.05e-04</td><td>1.39e-09</td><td>0.62</td><td>1.69e-06</td><td>1.53e-12</td><td>0.65</td><td>88.9k</td><td>1.13%</td></tr><tr><td>[L60fineA]</td><td>0.90</td><td>1.15e-04</td><td>2.01e-08</td><td>0.71</td><td>3.54e-06</td><td>5.33e-11</td><td>0.74</td><td>83.7k</td><td>0.52%</td></tr><tr><td>[L60fineB]</td><td>0.93</td><td>2.62e-04</td><td>2.01e-08</td><td>0.75</td><td>1.07e-06</td><td>5.08e-11</td><td>0.77</td><td>85.3k</td><td>0.29%</td></tr><tr><td>[L60fineC]</td><td>0.95</td><td>1.60e-04</td><td>2.01e-08</td><td>0.76</td><td>8.88e-07</td><td>5.33e-11</td><td>0.79</td><td>99.2k</td><td>0.16%</td></tr><tr><td>[L80geo]</td><td>0.84</td><td>1.11e-04</td><td>3.58e-09</td><td>0.61</td><td>6.90e-08</td><td>1.36e-13</td><td>0.63</td><td>216k</td><td>1.30%</td></tr><tr><td>[L100geo]</td><td>0.84</td><td>9.01e-04</td><td>5.59e-09</td><td>0.61</td><td>6.41e-07</td><td>4.05e-13</td><td>0.63</td><td>448k</td><td>1.47%</td></tr><tr><td>[L130geo]</td><td>0.83</td><td>1.04e-04</td><td>9.45e-09</td><td>0.60</td><td>1.90e-08</td><td>5.51e-13</td><td>0.63</td><td>965k</td><td>1.61%</td></tr><tr><td>[L160geo]</td><td>0.82</td><td>1.00e-04</td><td>1.43e-08</td><td>0.59</td><td>2.56e-09</td><td>1.01e-12</td><td>0.62</td><td>1.89M</td><td>1.84%</td></tr></tbody></table>

## 4.2 Global linear systems with permeameter boundary conditions

For this work, we consider the same hydrogeological settings as in [55]. We recall that we use the following permeameter boundary conditions (BCs): a pressure difference is imposed between two opposite faces of the cube (the same for the rock matrix and the fracture network) and zero-flux on the other sides. Furthermore, zero source terms, both in the rock matrix and the fracture network are considered. We also recall that we have two settings for the hydraulic conductivity in the rock matrix and the transmissivity in the fractures:

* “1v100”, where the hydraulic conductivity in the rock matrix is chosen equal to $K_\alpha^m = I_3 \text{ m/s}$, for all connected components $\omega_\alpha$, and where the fracture tangential transmissivity is chosen equal to $K_\ell^f = 100 I_2 \text{ m}^2/\text{s}$, for each fracture $\gamma_\ell$ (the same for all fractures);

* “heter”, where the hydraulic conductivity in the rock matrix is chosen equal to $K_\alpha^m = 10^{-8} I_3 \text{ m/s}$, for all connected components $\omega_\alpha$, and where the fracture tangential transmissivity $K_\ell^f$ takes one value per fracture $\gamma_\ell$ in the range $[10^{-6}, 20] \text{ m}^2/\text{s}$,

where $I_d$ is the identity matrix of size $d$.

Each test case is labeled by combining the name of a mesh from Table 1 in Section 4.1 and either “1v100” or “heter” depending on the permeability/transmissivity setting.

<table><caption>Table 3: Properties of the linear systems: number of dofs and number of non-zero elements.</caption><thead><tr><th>Test cases</th><th># dofs</th><th>nnz</th></tr></thead><tbody><tr><td>[L60geo-1v100] and [L60geo-heter]</td><td>18,525,914</td><td>134,224,135</td></tr><tr><td>[L60fineA-1v100]</td><td>36,658,783</td><td>263,406,425</td></tr><tr><td>[L60fineB-1v100]</td><td>65,907,706</td><td>471,436,936</td></tr><tr><td>[L60fineC-1v100]</td><td>128,881,985</td><td>918,288,659</td></tr><tr><td>[L80geo-1v100] and [L80geo-heter]</td><td>38,986,205</td><td>282,941,795</td></tr><tr><td>[L100geo-1v100] and [L100geo-heter]</td><td>71,758,072</td><td>521,157,724</td></tr><tr><td>[L130geo-1v100] and [L130geo-heter]</td><td>141,421,137</td><td>1,028,204,671</td></tr><tr><td>[L160geo-1v100] and [L160geo-heter]</td><td>242,616,910</td><td>1,764,733,010</td></tr></tbody></table>

The properties of the linear systems built with nef-flow-fpm can be found in Table 3. As in [55], the maximum ratio of heterogeneity is also denoted here by $K_{\text{ratio}}^{m,f}$. For the “1v100” test cases, it is equal to 100. For the “heter” test cases, even if the exact value depends on the chosen geometry, it is always close to $10^7$. We emphasize that the linear systems presented in Table 3 are highly ill-conditioned. Although computing the (estimated) condition number is not carried out in this work, it has been shown in [55] that fractured-porous media single-phase flow linear systems tend to be very ill-conditioned, even with a small number of fractures. In order to give an idea, the
condition number of [L60geo-1v100] was estimated in [55] to be 5.4 · 10¹². For a smaller “heter”
test case with a cube of 20, not used here, the condition number was observed to be 6.9 · 10¹⁶. We
expect the condition number of the test cases used in this study to be much worse, especially for
the largest [L160] test case with 697k fractures with 243M dofs.

5 Numerical results

In this section we present the results of a set of numerical experiments for flow in fractured
porous media. The simulations presented hereafter are run on the supercomputer Joliot-Curie
located at CEA's Very Large Computing Center (TGCC), using the Rome partition. That partition
consists of 2292 dual-processor AMD EPYC Rome computer nodes at 2.6 GHz with 64 cores per
processor, and 256 GB DDR4 memory/node.

Due to the large number of simulations that were run, the parser launcher pruners [57] is used.
pruners generates the input parameter files, launches the simulations and organizes the results in
a convenient folder hierarchy. The parameters used for testing the preconditioners with PETSc
are specified in A. We report wall-clock times for the solver phase, and we report separately the
times for the preconditioner setup and for the Krylov method iterations. The times are obtained
from PETSc log files. They are respectively the “PCSetUp”, and “KSPSolve” log event in PETSc
lingo. For the remainder of this work, we indicate these two times for every simulation. The total
solver time is essentially the sum of these times.

We note that the relative tolerance for all the Krylov methods presented hereafter is set to
10⁻¹⁴. The main reason for such a tight tolerance is to make sure that computed solutions respect
the flux continuity condition between two elements as much as possible, which is, as explained
in [55], a recurrent local mass conservation issue for single-phase flow simulations in fractured
porous media with a large number of fractures. The maximum number of iterations for iterative
methods is set to 10,000.

**5.1 Methodology**

We carry out a sequence of four numerical experiments. The first compares the performance
of several preconditioners available in PETSc, demonstrating that GenEO-RAS outperforms the
other preconditioners tested. The remaining three investigate the behavior of GenEO-RAS under
various configurations.

**COMPARISON** Run different parallel preconditioners, compare their performance, in terms and number of iterations and wall-clock time, and pick the best one. In this study, the tested preconditioners are:

* Algebraic Multi-Grid (AMG) methods: GAMG and BoomerAMG. A brief description of both methods, as well as a list of the parameters used, can be found in A.2.

* One-level DD: ASM and RAS. A list of the parameters passed to this preconditioner can be found in A.3;

* Two-level GenEO-ASM and GenEO-RAS. A list of the parameters passed to this pre-
conditioner can be found in A.4.

* Two-level block_splitting-ASM and block_splitting-RAS. A list of the parameters passed to this preconditioner can be found in A.5.

This first test will show that the combination of GMRES with GenEO-RAS as a preconditioner
gives the best performance. We next examine the performance of this method in more detail.

The following two tests are designed to study the scalability of GenEO-RAS.
**STRONG SCALING** We fix a test case (see Section 5.3), so that the size of the linear system is kept constant, and increase the number of MPI processes. Because Neumann matrices are needed for GenEO-RAS, the number of subdomains created by METIS in nef-flow-fpm must match the number of MPI processes.

**WEAK SCALING** Here we increase the number of MPI processes, while keeping the number of dofs per subdomain approximately constant. In practice, the size of the linear system is increased by refining the mesh multiple times and assembling the corresponding system with nef-flow-fpm. For each refined mesh, the number of MPI processes is chosen such that the number of degrees of freedom per subdomain remains nearly uniform across all simulations.

The fourth test aims at assessing the robustness of GenEO-RAS with respect to increasing the size of the DFM, and also the heterogeneity ratio $K_{ratio}^{\text{gen,lf}}$ in a way consistent with hydrogeological modeling.

**ROBUSTNESS WITH RESPECT TO DFM COMPLEXITY** We increase the size of the sub-cubes so as to get DFM of increasing complexity and we consider, for each DFM, both “1v100” and “heter” properties. For a fair comparison, we select the number of MPI processes, so that the number of degrees of freedom per subdomain remains nearly uniform across all simulations.

## 5.2 COMPARISON test

For the COMPARISON test, we choose the medium-sized DFM cases [L60geo-1v100] and [L60geo-heter], featuring 39k fractures and 18.5M degrees of freedom (dofs). A detailed description of the mesh is provided in Table 1, and the properties of the associated linear systems are given in Table 3 of Section 4.

Since the linear system is square, sparse, and SPD, the first solver tested is naturally CG. However, whenever RAS is used, either as a one-level preconditioner or as the fine-level component of HPDDM, the preconditioned system is no longer SPD, so that GMRES is employed, as it is suitable for both symmetric, and non-symmetric matrices. Note that we also tested GMRES in the symmetric case, possibly expecting better wall-clock times and/or number of iterations. The restart size for GMRES is set to 90 for the remainder of this work. We recall that GMRES relative residuals are decreasing with respect to the number of iterations, making the last iteration the one with the smallest relative residual value. This is not necessarily true for CG.

We also recall that the number of MPI processes is set equal to the number of subdomains. For this COMPARISON test, this number is set to 522.

We remind the reader that both HPDDM methods need to solve an auxiliary linear system for the coarse space correction (see (2.7)). The associated linear system is much smaller than the original system, but as we will see can still become fairly large. This system is solved with a direct solver. In HPDDM, the parallel direct solver MUMPS [5] is used. As sparse direct solvers do not scale as well as iterative solvers with a fined-grain workload, the coarse linear system is solved on a smaller number of processes, and HPDDM requires that number as input. In this section, we fixed that number to 4 independently of the number of subdomains. Furthermore, while the calculation of the number of deflation vectors is automatic in GenEO, block_splitting requires that number as an input. We chose that number equal to the maximum number of deflation vectors computed by GenEO. In practice, for [L60geo-1v100], GenEO computed 107 deflation vectors at most for a given subdomain, and for [L60geo-heter], GenEO computed 111 deflation vectors at most for a given subdomain.
### 5.2.1 Results for the test case [L60-1v100]

Table 4 presents, for the test case [L60-1v100], the number of iterations and the wall-clock time for CG and GMRES preconditioned with GAMG, BoomerAMG, ASM, RAS, block_splitting. Of course, RAS is omitted when the Krylov space method is CG.

<table><caption>Table 4: COMPARISON test - [L60geo-1v100] - time and iterations for CG and GMRES with various preconditioners.</caption><thead><tr><th rowspan="2">Preconditioner</th><th colspan="3">CG</th><th colspan="3">GMRES</th></tr><tr><th># it.</th><th>PCSetUp</th><th>KSPSolve</th><th># it.</th><th>PCSetUp</th><th>KSPSolve</th></tr></thead><tbody><tr><td>GAMG</td><td>2,303</td><td>00:02</td><td>02:27</td><td>&lt;!&gt;</td><td>&lt;!&gt; 00:02</td><td>&lt;!&gt; 12:23</td></tr><tr><td>BoomerAMG</td><td>4,621</td><td>00:01</td><td>01:50</td><td>&lt;!&gt;</td><td>&lt;!&gt; 00:01</td><td>&lt;!&gt; 05:41</td></tr><tr><td>ASM</td><td>635</td><td>00:02</td><td>00:23</td><td>1,396</td><td>00:02</td><td>01:00</td></tr><tr><td>RAS</td><td>—</td><td>—</td><td>—</td><td>924</td><td>00:02</td><td>00:41</td></tr><tr><td>block_splitting-ASM</td><td>67</td><td>00:45</td><td>00:25</td><td>65</td><td>00:46</td><td>00:19</td></tr><tr><td>block_splitting-RAS</td><td>—</td><td>—</td><td>—</td><td>72</td><td>00:47</td><td>00:20</td></tr><tr><td>GenEO-ASM</td><td>59</td><td>00:32</td><td>00:15</td><td>58</td><td>00:33</td><td>00:09</td></tr><tr><td>GenEO-RAS</td><td>—</td><td>—</td><td>—</td><td>38</td><td>00:33</td><td>00:06</td></tr></tbody></table>

The symbol "<!>" in Table 4 means that the simulation stopped at the maximum number of iterations, equal to 10k, not converging down to the requested tolerance value of $10^{-14}$. This only happened with GMRES preconditioned with GAMG and BoomerAMG. For GMRES-GAMG, the minimum residual reached is $5.5 \cdot 10^{-12}$, and for GMRES-BoomerAMG, the minimum residual reached is $3.3 \cdot 10^{-11}$.

Figure 5 shows for the test case [L60geo-1v100] the relative residual versus the number of iterations (truncated at 650 iterations) of CG and GMRES for these various preconditioners.

<figure><img src="image_9.png" /><figcaption>Figure 5: COMPARISON test - [L60geo-1v100] - (left) CG convergence and (right) GMRES convergence with various preconditioners - relative residual (log-scale) versus the number of iterations (linear scale, truncated at 650 iterations).</figcaption></figure>

In this [L60geo-1v100] case, we observe that:

* the AMG preconditioners converge with CG, however they both took more than 2,000 iterations, and they do not lead GMRES to converge;

* the best one-level preconditioner was ASM used with CG;

* the two-level DD preconditioners, used either with CG or GMRES, performed well, as they always converged in less than 72 iterations.
For comparison, solving [L60geo-1v100] with a direct solver, specifically Cholesky factorization, required 107 GiB of memory and took 6min26s [55]. The computation was performed using MUMPS [5, 6] with 48 MPI processes.

It is clear that, on this first test case, the best overall combination was GenEO-RAS used with GMRES, both in terms iteration counts (38) and total wall-clock time (39s).

This comparison demonstrates that, for low permeability/tranmissivity constraints, two-level DD preconditioners should be the preferred method in order to help the linear solvers to converge.

### 5.2.2 Results for the test case [L60-heter]

Table 5 presents, for the test case [L60-heter], the number of iterations and the wall-clock time for CG and GMRES preconditioned with GAMG, BoomerAMG, ASM, RAS, block_splitting, and GenEO. We again omit RAS when CG is used.

<table><caption>Table 5: COMPARISON test - [L60geo-heter] - time and iterations for CG and GMRES with various preconditioners. Symbol `<!>` indicates lack of convergence and symbol `@@` means that GMRES suffered from a breakdown (the residual increased after a restart).</caption><thead><tr><th rowspan="2">Preconditioner</th><th colspan="3">CG</th><th colspan="3">GMRES</th></tr><tr><th># it.</th><th>PCSetUp</th><th>KSPSolve</th><th># it.</th><th>PCSetUp</th><th>KSPSolve</th></tr></thead><tbody><tr><td>GAMG</td><td>&lt;!&gt;</td><td>&lt;!&gt; 00:03</td><td>&lt;!&gt; 10:39</td><td>&lt;!&gt;</td><td>&lt;!&gt; 00:02</td><td>&lt;!&gt; 12:30</td></tr><tr><td>BoomerAMG</td><td>9,227</td><td>00:01</td><td>03:27</td><td>@@ 2,430</td><td>@@ 00:01</td><td>@@ 03:32</td></tr><tr><td>ASM</td><td>&lt;!&gt;</td><td>&lt;!&gt; 00:02</td><td>&lt;!&gt; 05:28</td><td>&lt;!&gt;</td><td>&lt;!&gt; 00:02</td><td>&lt;!&gt; 06:47</td></tr><tr><td>RAS</td><td>—</td><td>—</td><td>—</td><td>@@ 1,080</td><td>@@ 00:02</td><td>@@ 00:47</td></tr><tr><td>block_splitting-ASM</td><td>74</td><td>00:49</td><td>00:28</td><td>69</td><td>00:48</td><td>00:20</td></tr><tr><td>block_splitting-RAS</td><td>—</td><td>—</td><td>—</td><td>52</td><td>00:48</td><td>00:17</td></tr><tr><td>GenEO-ASM</td><td>65</td><td>00:35</td><td>00:16</td><td>60</td><td>00:35</td><td>00:10</td></tr><tr><td>GenEO-RAS</td><td>—</td><td>—</td><td>—</td><td>39</td><td>00:35</td><td>00:06</td></tr></tbody></table>

The symbol `<!>` in Table 5 means that the simulation stopped at the maximum number of iterations, equal to 10k, not converging down to the requested tolerance value of $10^{-14}$. This only happened with the preconditioners GAMG and ASM. For CG-GAMG, the minimum residual is $1.4 \cdot 10^{-9}$ obtained at iteration 9,919 and for GMRES-GAMG, the minimum residual is $4.4 \cdot 10^{-11}$. For CG-ASM, the minimum residual is $8.5 \cdot 10^{-8}$ obtained at iteration 9,891 and for GMRES-ASM the minimum residual is $9.9 \cdot 10^{-10}$.

The symbol `@@` in Table 5 means that GMRES suffered from a breakdown (the residual increase after a restart). This happened with RAS and BoomerAMG preconditioners. Even if over 1,000 iterations were performed, the iterative solver stopped at iteration 1,080, with relative residual equal to $5.0 \cdot 10^{-8}$. The same goes for BoomerAMG, where the iterative solver stopped at iteration 2,430 with relative residual equal to $5.9 \cdot 10^{-11}$. Thus, none of these two preconditioners made GMRES converge down to the requested tolerance value of $10^{-14}$.

Figure 6 shows the relative residual versus the number of iterations of CG and GMRES for these various preconditioners.

In the “heter” case, we observe that:

* GAMG preconditioner and the one-level DD ASM preconditioner stopped at the maximum number of iterations, equal to 10k, not converging down to the requested tolerance value of $10^{-14}$;

* BoomerAMG makes CG converge with 9,227 iterations;

* the one-level DD RAS and BoomerAMG make GMRES breakdown;

* block_splitting-ASM converges in less than 80 iterations, both with CG and GMRES;
<figure><img src="image_10.png" /><figcaption>Figure 6: COMPARISON test - [L60geo-heter] - (left) CG convergence and (right) GMRES convergence with various preconditioners - relative residual (log-scale) versus the number of iterations (linear scale, truncated at 650 iterations).</figcaption></figure>

* better performance are obtained with GMRES and block_splitting-RAS with only 52 iterations (25% less iterations than block_splitting-ASM) and in very reasonable total wall-clock time (around one minute);

* GenEO-RAS took only 39 iterations and a wall-clock time of 41s. GenEO-RAS saves 35% of iterations with respect to GenEO-ASM.

In summary, the best performance for the “heter” test case was obtained with the combination GMRES with GenEO-RAS. These results confirm that two-level DD methods are the preferred candidates in order to precondition such large and ill-conditioned linear systems.

### 5.2.3 Selection of the most efficient preconditioner

The inadequacy of AMG methods for such simulations is consistent with previous findings [55], and is most likely due to the discretization method, which is not well suited for AMG solvers without appropriate intergrid transfer operators [23, 42, 72]. Additionally, the presence of low-quality elements in the mesh, as discussed in [29], may further contribute to their poor performance.

The one-level DD methods ASM and RAS converge with around 1k iterations for the “1v100” case but fail to converge for strong $K_{ratio}^{m,f}$ ratio.

Despite the additional work required for the decomposition with overlap and the construction of Neumann matrices, the two-level DD preconditioner GenEO-RAS has shown the best performance, both in terms of iteration counts and total wall-clock time, for any $K_{ratio}^{m,f}$ configuration. Even though the theory described in Section 2 is only valid for ASM, in practice, we observe here that the number of iterations is lower with RAS than with ASM, as explained in [36].

From this COMPARISON test, we conclude that the most efficient preconditioner is GenEO-RAS, hereafter simply named as GenEO. It will be used with GMRES as GMRES was shown to perform better than CG, especially with GenEO. This is the combination that will be retained for the remaining tests.

## 5.3 STRONG SCALING test

For this STRONG SCALING test, we continue with the DFM [L60geo-1v100], featuring 39k fractures and 18.5M degrees of freedom (dofs) and presented in Tables 1 and 3. As announced at the end of the previous Section 5.1, we now only test GMRES preconditioned with GenEO (with RAS for the fine level). In this test, the linear system remains unchanged, while the number of MPI processes is increased. As a result, the number of degrees of freedom per subdomain is decreasing whereas the size of the global coarse space is increasing. We compare results obtained
using 64, 128, 256, 512, 1024, and 2048 MPI processes. The number of processes allocated to the direct solver for the coarse space is set to 2, 2, 3, 4, 6, and 7, respectively.

It is worth noting that, with GenEO, changing the number of MPI processes requires generating a new set of Neumann matrices for each configuration, even though the left- and right-hand sides of the linear system remain unchanged across simulations.

Table 6 presents the results obtained with GMRES preconditioned with GenEO. We report in particular, the number of iterations, the wall-clock time, as well as some statistics on the number of computed eigenvalues and finally the coarse-space size.

<table><caption>Table 6: STRONG SCALING test - [L60geo-1v100] - Time and iterations of GMRES preconditioned with GenEO, statistics on the number of computed eigenvalues and coarse-space size.</caption><thead><tr><th># MPI procs.</th><th># dofs / sbdm.</th><th># it.</th><th>PCSetUp mm:ss</th><th>KSPSolve mm:ss</th><th># ev / sbdm max</th><th>min</th><th>median</th><th>Coarse space size</th></tr></thead><tbody><tr><td>64</td><td>289,467</td><td>35</td><td>18:47</td><td>00:40</td><td>368</td><td>134</td><td>278</td><td>15.3k</td></tr><tr><td>128</td><td>144,734</td><td>36</td><td>04:37</td><td>00:17</td><td>249</td><td>64</td><td>155</td><td>20.6k</td></tr><tr><td>256</td><td>72,366</td><td>36</td><td>01:25</td><td>00:08</td><td>161</td><td>56</td><td>137</td><td>27.3k</td></tr><tr><td>512</td><td>36,183</td><td>38</td><td>00:32</td><td>00:05</td><td>115</td><td>33</td><td>54</td><td>35.6k</td></tr><tr><td>1024</td><td>18,092</td><td>40</td><td>00:16</td><td>00:05</td><td>79</td><td>15</td><td>49</td><td>46.9k</td></tr><tr><td>2048</td><td>9,046</td><td>39</td><td>00:13</td><td>00:05</td><td>58</td><td>8</td><td>26</td><td>60.8k</td></tr></tbody></table>

Figure 7 presents the relative residual versus the number of iterations of GMRES for this strong scalability study.

<figure><img src="image_11.png" /><figcaption>Figure 7: STRONG SCALING test - [L60geo-1v100] - GMRES preconditioned with GenEO-relative residual (log-scale) versus the number of iterations (linear scale).</figcaption></figure>

### 5.3.1 Strong scalability of GenEO

Table 6 shows that the number of iterations remains remarkably stable, between 35 and 40. This means that GenEO is numerically strongly scalable, as could be expected from Theorem 2.1, even though RAS is used instead of ASM. Note the number of iterations is below the restart parameter of GMRES.

### 5.3.2 Wall-clock time and choice of the number of dofs per subdomain

We now look at how the wall-clock time changes when the number of MPI processes is increased. The setup time is the dominant contribution. The emerging behavior is somewhat surprising: one observes that, as the number of processes is doubled, the ratios in wall-clock time are 4, 3, 2.3, 1.5, and 1.2, so that the speedups observed for the smaller processor counts (from 64 to 512
processes) are superlinear. This effect can likely be explained by the fact that the complexity of
the local factorization is also superlinear in the number of unknowns. When the local number of
dofs is divided by two, the runtime is divided by more than two (close to four, corresponding to
a quadratic complexity for the sparse factorization). On the other hand, the relative slowdown
for larger processes counts is expected for a fixed-size problem, as the number of communication
between processes grows with respect to the size of the local computations.

One can deduce a *sweet spot* beyond which increasing the number of subdomains, and thus the computational resources, is no longer beneficial. The configuration with 512 subdomains is close to the “optimal” choice, leading to a number of dofs per subdomain around 36k. In the following experiments, in Sections 5.4 and 5.5, this is the approximate number we use.

5.3.3 Wall-clock time and size of the coarse space

Table 6 also presents some statistics on the distribution of eigenvalues across the subdomains.

In all experiments presented in this paper we have favored a control on the condition number of the preconditioned system, which as shown above (and also in later results) leads to a stable number of iterations with respect to the discretization parameters and the number of subdomains. Since it may be difficult to estimate in advance how many eigenvalues should be retained, we have made use of a recent feature of GenEO that automatically estimates the number of eigenvalues less than the threshold $\nu$ introduced in Section 2.3.1 by computing the inertia of the matrix pencil, see [46, Thm. 4.5.8]. Consequently, the number of computed eigenvalues for each subdomain, and thus the size of the coarse space, cannot be determined in advance. We have not attempted to tune the value of $\nu$ for each test case, nor to investigate the trade-off that would be incurred by keeping a fixed number of eigenvalues, while accepting an increase in the number of iterations.

The size of the coarse space does increase with the number of MPI processes, increasing from approximately 15k to 60k, representing a growth factor of about 4, compared to a 32-fold increase in the number of processes. This is due to the fact that, as subdomains get smaller, they become less heterogeneous and contribute fewer eigenvectors. The influence of this variation on the wall-clock time comes from the fact the coarse space matrix is factorized using a direct solver, which is known not to scale ideally with problem size. However, the number of processes used for this factorization is kept small, and only increases from two to seven.

To conclude, GenEO demonstrates strong scalability with respect to the number of iterations. Furthermore, choosing around 36k dofs per subdomain ensures a *sweet spot* where the wall-clock time is well-balanced with the number of MPI processes, ensuring an efficient use of computational resources.

5.4 WEAK SCALING test

For this WEAK SCALING test (see Section 5.1), we investigate the performance of GenEO with the following test cases: [L60geo-1v100] with 18M dofs, [L60fineA-1v100] with 36M dofs, [L60fineB-1v100] with 65M dofs, and [L60fineC-1v100] with 128M dofs. Their main properties can be found in Tables 1 and 3. The number of MPI processes for each case is selected so that the number of degrees of freedom per subdomain remains approximately constant at 36k, as motivated by the discussion in Section 5.3.2. The test cases [L60geo-1v100], [L60fineA-1v100], [L60fineB-1v100], and [L60fineC-1v100] are called with the following values for the number of processes used for the factorization of the coarse matrix: 4, 8, 12, and 24 respectively (see A.4 for details).

Table 7 presents the results obtained with GenEO, in particular, the clock time and the number of iterations.

Figure 8 presents the relative residual versus the number of iterations of GMRES for this weak scalability study.

From Table 7 and Figure 8, we observe that GenEO is weakly scalable in terms of iterations. In fact, even if the problem size grows for the same geometry due to refinement, the number of
<table><caption>Table 7: WEAK SCALING test - time and iterations of GMRES preconditioned with GenEO, statistics on the number of computed eigenvalues and coarse-space size.</caption><thead><tr><th>Test case</th><th># MPI<br>procs.</th><th># it.</th><th>PCSetUp<br>mm:ss</th><th>KSPSolve<br>mm:ss</th><th># ev / sbdm<br>max</th><th>min</th><th>median</th><th>Coarse<br>space<br>size</th></tr></thead><tbody><tr><td>[L60geo-1v100]</td><td>522</td><td>38</td><td>00:33</td><td>00:06</td><td>107</td><td>25</td><td>90</td><td>36k</td></tr><tr><td>[L60fineA-1v100]</td><td>1,033</td><td>39</td><td>00:42</td><td>00:08</td><td>100</td><td>32</td><td>83</td><td>70k</td></tr><tr><td>[L60fineB-1v100]</td><td>1,857</td><td>39</td><td>00:50</td><td>00:10</td><td>102</td><td>27</td><td>76</td><td>116k</td></tr><tr><td>[L60fineC-1v100]</td><td>3,631</td><td>39</td><td>01:03</td><td>00:14</td><td>89</td><td>26</td><td>62</td><td>207k</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>Figure 8: Weak scaling - GMRES preconditioned “GenEO” - relative residual (log-scale) versus the number of iterations (linear scale) for the different test cases: [L60geo-1v100] with 18M dofs, [L60fineA-1v100] with 36M dofs, [L60fineB-1v100] with 65M dofs and [L60fineC-1v100] with 128M dofs.</figcaption></figure>

iterations remains stable, less than 39 iterations. There is no effect of the restart, set to 90, in these simulations.

Nonetheless, we note that the solver wall-clock time does not scale ideally. In this case, the increase in the setup times can be attributed to the increase in the size of the coarse space (from 36k to 207k). The number of processes used for the coarse factorization increased at about the same rate as the size of the coarse space (by a factor of approximately six), but the complexity of the factorization is superlinear (see also the discussion in Section 5.3.2), so that the overall time increases. The number of eigenvalues per subdomain remains more or less stable, so that the size of the coarse space is approximately proportional to the total number of subdomains.

## 5.5 ROBUSTNESS WITH RESPECT TO DFM COMPLEXITY test

We now assess the performance of GMRES preconditioned with GenEO on DFMs of increasing complexity. The considered DFMs go from [L60] with 40k fractures and 18M dofs to [L160] with 696k fractures and 242M dofs. We refer the reader to Tables 1 and 3 for detailed descriptions of all the test cases. Additionally, each test case is run both in the “1v100” and the “heter” configurations.

As motivated in Section 5.3.2, the number of MPI processes is chosen so that the number of dofs per subdomain is approximately 36k. Furthermore, the number of MPI processes used for the factorization of the coarse matrix is 4, 8, 16, 32, and 96 for [L60], [L80], [L100], [L130], and [L160] respectively.

Table 8 shows the result obtained with GenEO, in particular, the number of iterations and the wall-clock time.

Figure 9 presents the relative residual versus the number of iterations of GMRES for the study of the impact of the geometry on the solver convergence.
<table><caption>Table 8: ROBUSTNESS WITH RESPECT TO DFM COMPLEXITY test - Time and iterations of GMRES preconditioned with GenEO for various test cases, statistics on the number of computed eigenvalues and coarse-space size.</caption><thead><tr><th>Test case</th><th># MPI procs.</th><th># it.</th><th>PCSetUp mm:ss</th><th>KSPSolve mm:ss</th><th># ev / sbdm max</th><th>min</th><th>median</th><th>Coarse space size</th></tr></thead><tbody><tr><td>[L60geo-1v100]</td><td>522</td><td>38</td><td>00:34</td><td>00:05</td><td>107</td><td>25</td><td>90</td><td>36k</td></tr><tr><td>[L80geo-1v100]</td><td>1,098</td><td>38</td><td>00:46</td><td>00:08</td><td>123</td><td>32</td><td>72</td><td>79k</td></tr><tr><td>[L100geo-1v100]</td><td>2,022</td><td>40</td><td>01:10</td><td>00:13</td><td>122</td><td>27</td><td>73</td><td>153k</td></tr><tr><td>[L130geo-1v100]</td><td>3,985</td><td>41</td><td>01:46</td><td>00:21</td><td>128</td><td>26</td><td>89</td><td>319k</td></tr><tr><td>[L160geo-1v100]</td><td>6,825</td><td>41</td><td>02:06</td><td>00:44</td><td>134</td><td>35</td><td>111</td><td>577k</td></tr><tr><td>[L60geo-heter]</td><td>522</td><td>39</td><td>00:35</td><td>00:06</td><td>111</td><td>26</td><td>93</td><td>37k</td></tr><tr><td>[L80geo-heter]</td><td>1,098</td><td>45</td><td>00:45</td><td>00:10</td><td>125</td><td>30</td><td>70</td><td>81k</td></tr><tr><td>[L100geo-heter]</td><td>2,022</td><td>46</td><td>01:16</td><td>00:17</td><td>124</td><td>27</td><td>72</td><td>155k</td></tr><tr><td>[L130geo-heter]</td><td>3,985</td><td>50</td><td>01:58</td><td>00:30</td><td>133</td><td>29</td><td>88</td><td>325k</td></tr><tr><td>[L160geo-heter]</td><td>6,825</td><td>51</td><td>02:23</td><td>01:00</td><td>139</td><td>33</td><td>107</td><td>587k</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Figure 9: ROBUSTNESS WITH RESPECT TO DFM COMPLEXITY test - GMRES preconditioned with GenEO for two configurations, “1v100” (right) and “heter” (left) - relative residual (log-scale) versus the number of iterations (linear scale) for different test cases: [L60] with 40k fractures and 18M dofs, [L80] with 90k fractures and 38M dofs, [L100] with 174k fractures and 71M dofs, [L130] with 377k fractures and 141M dofs, and [L160] with 696k fractures and 242M dofs.</figcaption></figure>

### 5.5.1 ROBUSTNESS WITH RESPECT TO DFM COMPLEXITY - GenEO iteration counts

From Table 8, and also from Figure 9, we observe that, for all test cases presented here, GMRES preconditioned with GenEO converges with a relatively low number of iterations (less than 51, below the restart of GMRES) and within a reasonable wall-clock time (less than 4min). While the number of iterations tends to increase as the number of fractures grows, the iteration counts remain generally stable. The increase is somewhat more pronounced in the “heter” configuration, but going from 39 to 51 iterations over a 13-fold increase in the number of processes.

These results show that GMRES with GenEO is very robust with respect to the number of fractures, up to 697k fractures and 243M dofs, and also with respect to high contrast in permeability/transmissivity. We note that a slight increase in the number of iterations as the fracture network becomes more complex could be expected from the bound 2.8, as both parameters $k_0$ and $k_1$ will likely increase.
### 5.5.2 ROBUSTNESS WITH RESPECT TO DFM COMPLEXITY - GenEO wall-clock time

Table 8 shows that the wall-clock time for the setup phase increases as the DFM size and complexity grow. This trend is also observed, and actually the growth is more pronounced, for the solve phase. Nevertheless, the setup phase still accounts for most of the wall-clock time. We recall that the number of dofs increased by a factor of approximately 13 (and the number of MPI processes increased accordingly), while the time for the setup phase rose by a factor of about four in both the “1v100” and the “heter” cases. Once again, the cause of the increased wall-clock time is the construction and the factorization of the coarse matrix. On the one hand, the size of the coarse space increased by a factor of 16 (while the number of processes by a factor of 24) and the time required to factorize the corresponding matrix grows more than linearly with its size. On the other hand, the factorizations of the local Dirichlet and Neumann matrices are performed on matrices of (roughly) the same size in all cases, since we have kept the number of dofs per MPI process approximately constant.

It is remarkable that GMRES preconditioned with GenEO is capable of solving very large test cases, the most challenging being [L160geo-heter], with 697k fractures and 243M dofs, solved in less than four minutes with 6,825 MPI processes.

## Conclusion

In this work, our main objective was to efficiently simulate steady-state single-phase flow in very large DFM models. The critical stage is the solution of the highly ill-conditioned linear system, containing hundreds of millions of unknowns, that arises from the mixed-hybrid finite element (MHFEM) discretization of the flow model. We investigated Krylov subspace methods such as CG or GMRES. As preconditioners, we evaluated algebraic multigrid (AMG) methods, one-level domain decomposition (DD) methods such as additive Schwarz method (ASM) and restricted ASM (RAS), and two-level DD methods equipped with coarse space operators, one based on a splitting method, and the other based on the generalized eigenvalue problem on the overlap (GenEO). To enable the use of GenEO in the context of DFM flows, we proved the compatibility of the discrete bilinear form arising from the MHFEM with the GenEO theoretical framework. We also proposed a decomposition strategy that takes into account the mixed-dimensional geometry of DFM s and allows for the construction of the so-called Neumann matrices and mappings required by GenEO.

The numerical experiments demonstrated that the GenEO preconditioner outperforms all other tested preconditioners, both in terms of iteration count and wall-clock time. Moreover, the GMRES-RAS combination outperformed CG-ASM, making it the preferred choice. GMRES preconditioned with GenEO showed both weak and strong scalability. As the geometric complexity of the DFM s increases by adding more fractures, the GenEO preconditioner remains efficient, even under strong permeability/transmissivity contrasts. Thanks to GenEO, the solver time is now reasonable, even for very large-scale problems. In the largest test case, with 697k fractures, 243M degrees of freedom, and 6,825 MPI processes, GMRES preconditioned with GenEO maintained robust performance, completing the simulation in less than four minutes with 51 iterations. These simulations represent a milestone in the field of fractured porous media.

A promising direction for future work is to use such fine-scale simulations to derive reduced flow models and upscaled quantities, as in [28]. Another research avenue is to address more complex physical phenomena, including heat transfer, solute transport, multiphase flow, and mechanics.

## Acknowledgements

This work was granted access to the GENCI-sponsored HPC resources of TGCC@CEA under allocation AD010607519R2. The authors are grateful to Simon Legrand for his support with pruners, and to Raphaël Zanella for helping in the development of the different PETSc launchers.
The authors also thank Fractory, a joint laboratory between ITASCA Consultants SAS, CNRS, and the University of Rennes, for supplying the geometries with very large number of fractures and the heterogeneous transmissivity field.

# A Parameters for the different preconditioners

## A.1 Common parameters

The following parameters are used for all simulations:

- ksp_type: the iterative solver type, which can be either `cg` or `gmres`;

- ksp_max_it: the maximum number of iterations;

- ksp_gmres_restart: only used when the Krylov method is GMRES, here it is set to 90;

- ksp_pc_side: the preconditioning side, set to `left` for CG and `right` for GMRES;

- ksp_norm_type: for the convergence test, for CG we choose the natural norm type $\frac{\|b-Ax\|_{M-1}}{\|b\|_{M-1}}$, while for GMRES we choose the unpreconditioned one $\frac{\|b-Ax\|_{2}}{\|b\|_{2}}$;

## A.2 AMG

We briefly recall the main features of the two Algebraic Multigrid Methods that were tested.

GAMG is the native PETSc algebraic multigrid framework, and is used here in the default smoothed aggregation mode, with one level of aggressive coarsening. The smoother on each level is Chebyshev iterations with SOR preconditioning. The smoothed aggregation needs eigenvalue estimates, which are provided by a few iterations of CG. The coarse grid solver is block Jacobi iterations, with a direct solver in each subdomain.

Besides the parameters presented in A.1, the following parameters were used:

- pc_type: GAMG is called with the value `gamg`;

- mgLevels_pc_type: the chosen smoother to use in each level is `sor`;

- mgLevels_esteig_ksp_type: by default, GAMG calls Chebyshev smoothers; when we use GAMG as a preconditioner for a CG solver, passing `cg` to this option makes the eigenvalue estimation better, thus avoiding possible indefinite preconditioner messages; otherwise, passing `gmres` does not perform this extra step;

- pc_gamg_aggressive_coarsening: the number of aggressive coarsening (Modified Independent Set MIS-2) levels from finest is set to 1.

BoomerAMG is an interface from PETSc to the HYPRE algebraic multigrid solver [37]. It is a parallel implementation of the Ruge-Stüben algebraic multigrid. It is used here with Parallel Maximal Independent Sets (PMIS) coarsening and $l^1$-scaled Jacobi smoothers on each level. For BoomerAMG, from the HYPRE library, besides the parameters presented in A.1, the following parameters were used:

- pc_type: we have to call first `hypre`;

- pc_hypre_type: then we call `boomeramg` as the preconditioner;

- pc_hypre_boomeramg_relax_type_all: the relaxation method is set to `l1scaled-Jacobi`;

- pc_hypre_boomeramg_coarsen_type: the coarsening method is set to `Parallel MIS pmis`.

## A.3 ASM and RAS

Besides the parameters presented in A.1, the parameters used for the one-level domain decomposition method are as follows:

- pc_type: the one-level domain decomposition preconditioner is called with the value `asm`;

- pc_asm_type: in this work we have two possibilities for the one-level preconditioner, either `basic` for ASM or `restrict` for RAS;

- sub_pc_factor Mat_solver_type: local problems are factored with the direct solver library `mumps`;
-sub_pc_type: to be more specific, the factorization is done thanks to cholesky.

## A.4 GenEO

Besides the parameters presented in A.1, we also pass the following parameters:

-pc_type: the two-level domain decomposition preconditioner is called with the value hpddm;

-pcmani: set to 1 to inform the preconditioner that the decomposition into subdomains is supplied;

-pcmani: set to 1 to inform the preconditioner that the Neumann matrices are supplied;

-pcmani: similarly to A.3, we have two possibilities for the fine-level preconditioner, either basic for ASM or restrict for RAS;

-pcmani: similarly to A.3, the factorization in the fine level is done thanks to cholesky;

-pcmani: similarly to A.3, local problems are factorized with the direct solver library mumps;

-pcmani: this option makes use the same factorization method as in the fine level for the GEVP (2.4);

-pcmani: it is the $\nu$ threshold for the GEVP mentioned in theorem 2.1, set is equal to 0.1 in this work;

-pcmani: the maximum number of computed eigenpairs in the GEVP is automatically computed by GenEO for each subdomain $\Omega_i^\delta$;

-pcmani: the factorization of the coarse matrix done thanks to cholesky;

-pcmani: the factorization of the coarse matrix is performed with the direct solver library mumps;

-pcmani: this options filters out the values below a threshold in the coarse operator, which is set to 1e-14 in this work;

-pcmani: the number of MPI processes used for aggregating the coarse operator;

-pcmani: this is the choice we have to make in (2.3), between a deflated or balanced two-level preconditioner; in this work, for CG it is balanced, while for GMRES it is deflated.

## A.5 block_splitting

Besides the parameters presented in A.1, the parameters passed to this preconditioner are exactly the ones from A.4, except for the following:

-pc: the value is set to 0, so no Neumann matrix is used;

-pc: this option is not available for HPDDM-splitting, thus is not used; instead, -pcmani is the maximum number of computed eigenpairs in the GEVP for all subdomains $\Omega_i^\delta$ and has to be manually set by the user.

-pc: most importantly, we also add this option to activate the library related to the method presented in [2].

## References

[1] E. Agullo, L. Giraud, and L. Poirel. Robust preconditioners via generalized eigenproblems for hybrid sparse linear solvers. SIAM J. Matrix Anal. Appl., 40(2):417-439, 2019. doi:10.1137/17M1153765.

[2] H. Al Daas, P. Jolivet, and T. Rees. Efficient Algebraic Two-Level Schwarz Preconditioner for Sparse Matrices. SIAM J. Sci. Comput., 45(3):A1199-A1213, 2023. doi:10.1137/22M1469833.

[3] C. Albin. Deux outils mathématiques pour modéliser l'écoulement et le transport de polluants dans un milieu poreux fracturé. PhD thesis (in French), Paris IX Dauphine University (France), 2000. URL: https://bu.dauphine.psl.eu/fileviewer/view.php?doc=2000PA090024.
[4] C. Alboin, J. Jaffré, J. E. Roberts, X. Wang, and C. Serres. Domain Decomposition for Some Transmission Problems in Flow in Porous Media. In Z. Chen, R. E. Ewing, and Z.-C. Shi, editors, *Numerical Treatment of Multiphase Flows in Porous Media*, pages 22–34. Springer Berlin Heidelberg, 2000. doi:10.1007/3-540-45467-5_2.

[5] P. R. Amestoy, I. S. Duff, J. Koster, and J.-Y. L'Excellent. A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling. *SIAM J. Matrix Anal. Appl.*, 23(1):15–41, 2001. doi:10.1137/S0895479899358194.

[6] P. R. Amestoy, I. S. Duff, D. Ruiz, and B. Uçar. A parallel matrix scaling algorithm. In J. M. L. M. Palma, P. R. Amestoy, M. Daydé, M. Mattoso, and J. C. Lopes, editors, *High Performance Computing for Computational Science - VECPAR 2008*, pages 301–313, Berlin, Heidelberg, 2008. Springer Berlin Heidelberg. doi:10.1007/978-3-540-92859-1_27.

[7] L. Amir, M. Kern, Z. Mghazli, and J. E. Roberts. Intersecting fractures in porous media: mathematical and numerical analysis. *Applicable Analysis*, 102(5):1312–1334, 2023. doi:10.1080/00036811.2021.1981878.

[8] P. Angot, F. Boyer, and F. Hubert. Asymptotic and numerical modelling of flows in fractured porous media. *ESAIM: Mathematical Modelling and Numerical Analysis*, 43(2):239–275, 2009. doi:10.1051/m2an/2008052.

[9] P. F. Antonietti, J. De Ponti, L. Formaggia, and A. Scotti. Preconditioning Techniques for the Numerical Solution of Flow in Fractured Porous Media. *Journal of Scientific Computing*, 86(1):2, 2020. doi:10.1007/s10915-020-01372-0.

[10] A. Arrarás, F. J. Gaspar, L. Portero, and C. Rodrigo. Mixed-Dimensional Geometric Multigrid Methods for Single-Phase Flow in Fractured Porous Media. *SIAM J. Sci. Comput.*, 41(5):B1082–B1114, 2019. doi:10.1137/18M1224751.

[11] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, H. Suh, S. Zampini, H. Zhang, H. Zhang, and J. Zhang. PETSc/TAO users manual. Technical Report ANL-21/39 – Revision 3.23, Argonne National Laboratory, 2025. doi:10.2172/2565610.

[12] S. Balay, W. D. Gropp, L. Curfman McInnes, and B. F. Smith. Efficient Management of Parallelism in Object Oriented Numerical Software Libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, *Modern Software Tools in Scientific Computing*, pages 163–202. Birkhäuser Press, 1997.

[13] I. Berre, W. M. Boon, B. Flemisch, A. Fumagalli, D. Gläser, E. Keilegavlen, A. Scotti, I. Stefansson, A. Tatomir, K. Brenner, S. Burbulla, P. Devloo, O. Duran, M. Favino, J. Hennicker, I-H. Lee, K. Lipnikov, R. Masson, K. Mosthaf, M. G. C. Nestola, C.-F. Ni, K. Nikitin, P. Schädle, D. Svyatskiy, R. Yanbarisov, and P. Zulian. Verification benchmarks for single-phase flow in three-dimensional fractured porous media. *Advances in Water Resources*, 147:103759, 2021. doi:10.1016/j.advwatres.2020.103759.

[14] I. Berre, F. Doster, and E. Keilegavlen. Flow in Fractured Porous Media: A Review of Conceptual Models and Discretization Approaches. *Transport in Porous Media*, 130(1):215–236, 2019. doi:10.1007/s11242-018-1171-6.

[15] S. Berrone, S. Pieraccini, and S. Scialò. Flow simulations in porous media with immersed intersecting fractures. *Journal of Computational Physics*, 345:768–791, 2017. doi:10.1016/j.jcp.2017.05.049.

[16] D. Boffi, F. Brezzi, and M. Fortin. *Mixed finite element methods and applications*, volume 44 of Springer Series in Computational Mathematics. Springer, Heidelberg, 2013. doi:10.1007/978-3-642-36519-5.

[17] E. Bonnet, O. Bour, N. E. Odling, P. Davy, I. Main, P. Cowie, and B. Berkowitz. Scaling of fracture systems in geological media. *Reviews of Geophysics*, 39(3):347–383, 2001. doi:10.1029/1999RG000074.

[18] W. M. Boon, J. M. Nordbotten, and I. Yotov. Robust discretization of flow in fractured porous media. *SIAM J. Numer. Anal.*, 56(4):2203–2233, 2018. doi:10.1137/17M1139102.

[19] H. Borouchaki, P. Laug, and P.-L. George. Parametric surface meshing using a combined advancing-front generalized delaunay approach. *International Journal for Numerical Methods in Engineering*, 49(1-2):233–259, 2000. doi:10.1002/1097-0207(20000910/20)49:1/2<233::AID-NME931>3.0.CO;2-G.

[20] K. Brenner, M. Groza, C. Guichard, G. Lebeau, and R. Masson. Gradient discretization of hybrid dimensional Darcy flows in fractured porous media. *Numerische Mathematik*, 134(3):569–609, 2016. doi:10.1007/s00211-015-0782-x.
[21] K. Brenner, J. Hennicker, R. Masson, and P. Samier. Gradient discretization of hybrid-dimensional Darcy flow in fractured porous media with discontinuous pressures at matrix-fracture interfaces. *IMA Journal of Numerical Analysis*, 37(3):1551–1585, 2017. doi:10.1093/imanum/drw044.

[22] K. Brenner, J. Hennicker, R. Masson, and P. Samier. Hybrid-dimensional modelling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions. *Journal of Computational Physics*, 357:100–124, 2018. doi:10.1016/j.jcp.2017.12.003.

[23] S. C. Brenner. A multigrid algorithm for the lowest-order Raviart–Thomas mixed triangular finite element method. *SIAM J. Numer. Anal.*, 29(3):647–678, 1992. doi:10.1137/0729042.

[24] F. Brezzi and M. Fortin. *Mixed and hybrid finite element methods*, volume 15 of Springer Series in Computational Mathematics. Springer-Verlag, New York, 1991. doi:10.1007/978-1-4612-3172-1.

[25] M. C. Cacas, E. Ledoux, G. de Marsily, B. Tillie, A. Barbreau, E. Durand, B. Feuga, and P. Peaudecerf. Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 1. the flow model. *Water Resources Research*, 26(3):479–489, 1990. doi:10.1029/WR026i003p00479.

[26] X.-C. Cai and M. Sarkis. A restricted additive Schwarz preconditioner for general sparse linear systems. *SIAM J. Sci. Comput.*, 21(2):792–797, 1999. doi:10.1137/S106482759732678X.

[27] L. M. Carvalho, L. Giraud, and G. Meurant. Local preconditioners for two-level non-overlapping domain decomposition methods. *Numerical Linear Algebra with Applications*, 8(4):207–227, 2001. doi:10.1002/nla.237.

[28] Eduardo S. Castro, Josué S. Barroso, Marcio A. Murad, and João N.C. Guerreiro. Reduced flow model and transmissibility upscaling in multi-layered faulted reservoirs. *Computers & Geosciences*, 180:105439, 2023. doi:10.1016/j.cageo.2023.105439.

[29] Y. Chen and G. N. Wells. Multigrid on unstructured meshes with regions of low quality cells, 2024. arXiv: 2402.12947, doi:10.48550/arXiv.2402.12947.

[30] C. Chevalier and F. Pellegrini. PT-SCOTCH: A tool for efficient parallel graph ordering. *Parallel Computing*, 34(6):318–331, 2008. Parallel Matrix Algorithms and Applications. doi:10.1016/j.parco.2007.12.001.

[31] P. Davy, R. Le Goc, and C. Darcel. A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling. *Journal of Geophysical Research: Solid Earth*, 118(4):1393–1407, 2013. doi: 10.1002/jgrb.50120.

[32] P. Davy, R. Le Goc, C. Darcel, O. Bour, J.-R. de Dreuzy, and R. Munier. A likely universal model of fracture scaling and its consequence for crustal hydromechanics. *Journal of Geophysical Research: Solid Earth*, 115(B10), 2010. doi:10.1029/2009JB007043.

[33] C. R. Dohrmann. An approximate BDDC preconditioner. *Numerical Linear Algebra with Applications*, 14(2):149–168, 2007. doi:10.1002/nla.514.

[34] V. Dolean, P. Jolivet, and F. Nataf. *An Introduction to Domain Decomposition Methods*. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015. doi:10.1137/1.9781611974065.

[35] Y. Efendiev, J. Galvis, R. Lazarov, and J. Willems. Robust domain decomposition preconditioners for abstract symmetric positive definite bilinear forms. *ESAIM Math. Model. Numer. Anal.*, 46(5):1175–1199, 2012. doi: 10.1051/m2an/2011073.

[36] E. Efstathiou and M. J. Gander. Why restricted additive Schwarz converges faster than additive Schwarz. *BIT*, 43:945–959, 2003. doi:10.1023/B:BITN.0000014563.33622.1d.

[37] R. D. Falgout, J. E. Jones, and U. M. Yang. The design and implementation of hypre, a library of parallel high performance preconditioners. In *Are Magnus Bruaset and Aslak Tveito, editors, Numerical Solution of Partial Differential Equations on Parallel Computers*, pages 267–294, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg. doi:10.1007/3-540-31619-1_8.

[38] GAMG: Geometric algebraic multigrid (AMG) preconditioner. https://petsc.org/release/manualpages/PC/GAMG/.

[39] G.N. Gatica. *A simple introduction to the mixed finite element method. Theory and applications*. SpringerBriefs in Mathematics. Springer, 2014. doi:10.1007/978-3-319-03695-3.
[40] P.-L. George, H. Borouchaki, F. Alauzet, A. Loseille, P. Laug, and L. Maréchal. *Maillage, modélisation géométrique et simulation numérique 2*. ISTE Group, 2018. URL: https://www.istegroup.com/fr/produit/maillage-modelisation-geometrique-et-simulation-numerique-2/.

[41] L. Giraud, A. Haidar, and L.T. Watson. Parallel scalability study of hybrid preconditioners in three dimensions. *Parallel Computing*, 34(6):363-379, 2008. doi:10.1016/j.parco.2008.01.006.

[42] J. Gopalakrishnan and S. Tan. A convergent multigrid cycle for the hybridized mixed method. *Numerical Linear Algebra with Applications*, 16(9):689-714, 2009. doi:10.1002/nla.636.

[43] Pascal Havé, Roland Masson, Frédéric Nataf, Mikolaj Szydlarski, Hua Xiang, and Tao Zhao. Algebraic domain decomposition methods for highly heterogeneous problems. *SIAM Journal on Scientific Computing*, 35(3):C284-C302, 2013. doi:10.1137/110842648.

[44] V. Hernandez, J. E. Roman, and V. Vidal. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems. *ACM Transactions on Mathematical Software*, 31(3):351-362, September 2005. doi: 10.1145/1089014.1089019.

[45] M.R. Hestenes and E. Stiefel. Methods of Conjugate Gradients for Solving Linear Systems. *Journal of Research of the National Bureau of Standards*, 49:409-435, 1952. doi:10.6028/jres.049.044.

[46] Roger A Horn and Charles R Johnson. *Matrix analysis*. Cambridge University Press, 2012. doi:10.1017/CBO9781139020411.

[47] H. Hoteit and A. Firoozabadi. An efficient numerical model for incompressible two-phase flow in fractured media. *Advances in Water Resources*, 31(6):891-905, 2008. doi:10.1016/j.advwatres.2008.02.004.

[48] hypre: High performance preconditioners. https://llnl.gov/casc/hypre, https://github.com/hypre-space/hypre.

[49] P. Jolivet, F. Hecht, F. Nataf, and C. Prud'homme. Scalable domain decomposition preconditioners for heterogeneous elliptic problems. In *Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis*, SC '13, New York, NY, USA, 2013. Association for Computing Machinery. doi:10.1145/2503210.2503212.

[50] P. Jolivet, J.E. Roman, and S. Zampini. KSPHPDDM and PCHPDDM: Extending PETSc with advanced Krylov methods and robust multilevel overlapping Schwarz preconditioners. *Computers and Mathematics with Applications*, 84:277-295, 2021. doi:10.1016/j.camwa.2021.01.003.

[51] P. Jolivet and P.-H. Tournier. Block iterative methods and recycling for improved scalability of linear solvers. In *Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis*, SC '16. IEEE Press, 2016.

[52] G. Karypis and V. Kumar. METIS: A software package for partitioning unstructured graphs, partitioning meshes, and computing fill-reducing orderings of sparse matrices. Technical report, Univ. of Minnesota, Dept. of Computer Sci. and Engr., 01 1997.

[53] G. Karypis, K. Schloegel, and V. Kumar. ParMETIS: Parallel graph partitioning and sparse matrix ordering library. Technical report, Univ. of Minnesota, Dept. of Computer Sci. and Engr., 01 1997.

[54] M. Kern, G. Pichot, and D. Zegarra Vasquez. Mathematical and numerical analysis of the mixed formulation of single phase flow in three-dimensional fractured porous media. Preprint, 2025. URL: https://inria.hal.science/hal-05029638.

[55] M. Kern, G. Pichot, and D. Zegarra Vasquez. Performance of algebraic preconditioners for large-scale simulations of single-phase flow in three-dimensional fractured porous media. Preprint, 2025. URL: https://inria.hal.science/hal-05029652.

[56] P. Laug and G. Pichot. Mesh Generation and Flow Simulation in Large Tridimensional Fracture Networks. In S. Carillo, C. Conti, D. Mansutti, F. Pitolli, and R. M. Spitaleri, editors, *IMACS Series in Computational and Applied Mathematics*, volume 22, pages 71-80, 2019. URL: https://hal.archives-ouvertes.fr/hal-02102811.

[57] S. Legrand, T. Martinez, and G. Pichot. pruners. *Git Repository*, 2022. URL: https://gitlab.inria.fr/pruners/pruners.

[58] S. B. Lunowa, I. S. Pop, and B. Koren. Linearized domain decomposition methods for two-phase porous media flow models involving dynamic capillarity and hysteresis. *Computer Methods in Applied Mechanics and Engineering*, 372:113364, 2020. doi:10.1016/j.cma.2020.113364.
[59] J. Mandel, B. Sousé, and C. R. Dohrmann. Multispace and multilevel BDDC. *Computing*, 83(2):55-85, 2008. doi:10.1007/s00607-008-0014-7.

[60] V. Martin, J. Jaffré, and J. E. Roberts. Modeling fractures and barriers as interfaces for flow in porous media. *SIAM J. Sci. Comput.*, 26(5):1667-1691, 2005. doi:10.1137/S1064827503429363.

[61] The MueLu Project Team. *The MueLu Project Website*, 2020 (accessed May 5, 2025). URL: https://trilinos.github.io/muelu.html.

[62] F. Nataf, H. Xiang, V. Dolean, and N. Spillane. A coarse space construction based on local Dirichlet to Neumann maps. *SIAM J. Sci. Comput.*, 33(4):1623-1642, 2011. doi:10.1137/100796376.

[63] B. Noetinger. A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks accounting for matrix to fracture flow. *Journal of Computational Physics*, 283:205-223, 2015. doi:10.1016/j.jcp.2014.11.038.

[64] F. Pellegrini. SCOTCH and PT-SCOTCH Graph Partitioning Software: An Overview. In U. Naumann and O. Schenk eds., editors, *Combinatorial Scientific Computing*, page 34. Chapman and Hall/CRC, 2012.

[65] F. Pellegrini and J. Roman. SCOTCH: A software package for static mapping by dual recursive bipartitioning of process and architecture graphs. In H. Liddell, A. Colbrook, B. Hertzberger, and P. Sloot (eds), editors, *High-Performance Computing and Networking. HPCN-Europe 1996. Lecture Notes in Computer Science.*, volume 1067, pages 493-498. Springer, Berlin, Heidelberg, 1996. doi:10.1007/3-540-61142-8_588.

[66] L. Poirel. *Algebraic Domain Decomposition Methods for Hybrid (direct/iterative) Solvers*. PhD thesis, University of Bordeaux (France), 2018. URL: https://theses.hal.science/tel-03555822v2.

[67] Y. Saad and M. H. Schultz. GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. *SIAM J. Sci. Comput.*, 7(3):856-869, 1986. doi:10.1137/0907058.

[68] N. Spillane, V. Dolean, P. Hauret, F. Nataf, C. Pechstein, and R. Scheichl. Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps. *Numerische Mathematik*, 126(4):741-770, 2014. doi:10.1007/s00211-013-0576-y.

[69] J.M. Tang, R. Nabben, C. Vuik, and Y.A. Erlangga. Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods. *J. Sci. Comput.*, 39(3):340-370, 2009. doi:10.1007/s10915-009-9272-6.

[70] A. Toselli and O. Widlund. *Domain decomposition methods—algorithms and theory*, volume 34 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 2005. doi:10.1007/b137868.

[71] M. Vohralík, J. Maryška, and O. Severýn. Mixed and nonconforming finite element methods on a system of polygons. *Applied Numerical Mathematics. An IMACS Journal*, 57(2):176-193, 2007. doi:10.1016/j.apnum.2006.02.005.

[72] T. Wildey, S. Muralikrishnan, and T. Bui-Thanh. Unified geometric multigrid algorithm for hybridized high-order finite element methods. *SIAM J. Sci. Comput.*, 41(5):S172-S195, 2019. doi:10.1137/18M1193505.