A Scalable Multigrid Reduction Framework for Multiphase Poromechanics of Heterogeneous Media

TL;DR

This paper introduces a scalable algebraic multigrid reduction framework tailored for efficiently solving the large, coupled linear systems arising in multiphase poromechanics simulations involving heterogeneous media.

Contribution

It presents a novel algebraic multigrid reduction approach that effectively decouples complex PDE systems for improved solver performance in multiphase poromechanics.

Findings

Framework is flexible for various scenarios.

Achieves efficiency and scalability for large problems.

Outperforms traditional black-box solvers.

Abstract

Simulation of multiphase poromechanics involves solving a multi-physics problem in which multiphase flow and transport are tightly coupled with the porous medium deformation. To capture this dynamic interplay, fully implicit methods, also known as monolithic approaches, are usually preferred. The main bottleneck of a monolithic approach is that it requires solution of large linear systems that result from the discretization and linearization of the governing balance equations. Because such systems are non-symmetric, indefinite, and highly ill-conditioned, preconditioning is critical for fast convergence. Recently, most efforts in designing efficient preconditioners for multiphase poromechanics have been dominated by physics-based strategies. Current state-of-the-art "black-box" solvers such as algebraic multigrid (AMG) are ineffective because they cannot effectively capture the strong coupling between the mechanics and the flow sub-problems, as well as the coupling inherent in the multiphase flow and transport process. In this work, we develop an algebraic framework based on multigrid reduction (MGR) that is suited for tightly coupled systems of PDEs. Using this framework, the decoupling between the equations is done algebraically through defining appropriate interpolation and restriction operators. One can then employ existing solvers for each of the decoupled blocks or design a new solver based on knowledge of the physics. We demonstrate the applicability of our framework when used as a "black-box" solver for multiphase poromechanics. We show that the framework is flexible to accommodate a wide range of scenarios, as well as efficient and scalable for large problems.

Equations

$$-\nabla \cdot \mathbf{\sigma }=\rho \mathbf{g}$$

$${\dot{m}}_w+\nabla \cdot {\mathbf{w}}_w=q_w$$

$${\dot{m}}_{nw}+\nabla \cdot {\mathbf{w}}_{nw}=q_{nw}$$

$${\mathcal{F}}_u$$

$${\mathcal{F}}_s$$

$${\mathcal{F}}_p$$

$$F({\mathbf{x}}_n)$$

$$A({\mathbf{x}}_n^{(k)})\Delta \mathbf{x}=-F({\mathbf{x}}_n^{(k)}),$$

$${\mathbf{x}}_n^{(k+1)}={\mathbf{x}}_n^{(k)}+\lambda \Delta \mathbf{x},$$

$$A=\left(\begin{array}{ccc}{A}_{uu}& A_{us}& A_{up}\\ A_{su}& A_{ss}& A_{sp}\\ A_{pu}& A_{ps}& A_{pp}\end{array}\right).$$

$$A{M}^{-1}\Delta \mathbf{y}=-F({\mathbf{x}}_n^{(k)}),$$

$$A=\left(\begin{array}{cc}{A}_{FF}& A_{FC}\\ A_{CF}& A_{CC}\end{array}\right)=\left(\begin{array}{cc}{I}_{FF}& 0\\ A_{CF}{A}_{FF}^{-1}& I_{CC}\end{array}\right)\left(\begin{array}{cc}{A}_{FF}& 0\\ 0& S\end{array}\right)\left(\begin{array}{cc}{I}_{FF}& A_{FF}^{-1}{A}_{FC}\\ 0& I_{CC}\end{array}\right),$$

$$P=\left(\begin{array}{c}-A_{FF}^{-1}{A}_{FC}\\ I_{CC}\end{array}\right),R=\left(\begin{array}{cc}-A_{CF}{A}_{FF}^{-1}& I_{CC}\end{array}\right).$$

$$A^{-1}=P(RAP{)}^{-1}R+Q(Q^{T}AQ{)}^{-1}{Q}^T,$$

$$0=I-A^{-1}A$$

$$=(I-P(RAP{)}^{-1}RA)(I-Q(Q^{T}AQ{)}^{-1}{Q}^{T}A)$$

$$=(I-Q(Q^{T}AQ{)}^{-1}{Q}^{T}A)(I-P(RAP{)}^{-1}RA),$$

$$\tilde{P}=\left(\begin{array}{c}{W}_p\\ I_{CC},\end{array}\right),\tilde{R}=\left(\begin{array}{cc}{W}_r& I_{CC}\end{array}\right).$$

$$W_r\approx -A_{CF}{A}_{FF}^{-1},W_p\approx -A_{FF}^{-1}{A}_{FC}.$$

$$W_r=0,W_p=-D_{FF}^{-1}{A}_{FC},$$

$$I-M_{MGR}^{-1}A$$

$$S=A_{CC}-\mathbf{G}(A_{CF}{D}_{FF}^{-1}{A}_{FC}).$$

$$A$$

$$A_{CC}=\left(\begin{array}{cc}{A}_{ss}& A_{sp}\\ A_{ps}& A_{pp}\end{array}\right).$$

$$S_1$$

$$=\left(\begin{array}{cc}{\tilde{A}}_{ss}& {\tilde{A}}_{sp}\\ {\tilde{A}}_{ps}& {\tilde{A}}_{pp}\end{array}\right).$$

$$S_1$$

$$S_2={\tilde{A}}_{pp}-{\tilde{A}}_{ps}{D}_{ss}^{-1}{\tilde{A}}_{sp},$$

Full text

1 Introduction

Modeling subsurface systems requires an understanding of many different physical processes, including multiphase fluid flow and transport, and geomechanical deformations. These processes are often tightly coupled in a “two-way” fashion: for example, the flow process can have a large influence on the mechanical process and in turn can be affected by the feedback of the produced mechanical response. To simulate these processes, one needs to solve a set of coupled, nonlinear, time-dependent partial differential equations (PDEs) that govern the conservation of mass of the fluid phases and linear momentum of the solid-fluid mixture. For this system, fully-implicit time discretization is the widely preferred approach because it is unconditionally stable and allows for large time steps. However, using an implicit approach, one must solve a large, sparse, and ill-conditioned linear system at each nonlinear iteration. Robust and scalable solvers are therefore needed for large scale simulations on high performance computing platforms. This paper presents our efforts to design an efficient preconditioning strategy based on an algebraic framework that is flexible and capable of addressing the inherent ill-conditioning of a complicated multi-physics system.

In recent years, much of the work in developing preconditioning strategies for coupled problems has focused on so-called physics-based strategies. The key is to use knowledge of the specific physical processes involved to break the tightly coupled systems into smaller sub-problems whose properties are well-studied. For example, these sub-problems could take the form of an elliptic, hyperbolic, or parabolic PDE, to which appropriate techniques can be applied. For fully implicit simulation of complex multiphase flow and transport without mechanics, one of the most popular methods is the Constrained Residual Pressure (CPR) multistage preconditioning technique [42, 43]. For single-phase flow poromechanics, many block preconditioners have been developed [1, 5, 6, 23, 28, 45, 47]. In the context of multiphase poromechanics, one recent strategy [46] uses the fixed-stress partitioning [25, 31, 37, 47] of the mechanics and the flow parts combined with a CPR approach [42] for the multiphase flow system. In general, these physics-based preconditioners are among the most effective techniques available. However, designing a good strategy is both time-consuming and challenging as it requires extensive knowledge of the particular continuous model of interest. One alternative is to use a “black-box” approach, such as algebraic multigrid (AMG) [39, 41]. AMG techniques are among the most efficient and scalable methods for solving sparse linear systems. Unlike geometric multigrid, these methods do not need an explicit hierarchy of computational grids. However, they are originally designed for scalar elliptic PDEs, and their applicability is much more limited for strongly coupled systems of PDEs.

Recently, multigrid reduction (MGR) [33, 34], a variant of AMG, has been applied successfully to coupled systems of multiphase flow and transport with phase transitions [8, 44]. Drawing on this success, in this work we further develop MGR into a general multi-level framework for solving discrete systems coming from discretization of tightly coupled PDEs. We also introduce a new dropping strategy for computing the reduction onto the coarse-grid within the MGR V-cycle that effectively captures the coupling between mechanics and flow. The goal of this strategy is two-pronged: (1) to keep the coarse grid sparse during multi-level reduction, and (2) to make the coarse grid amenable to classical AMG. We show that with this new feature, MGR is effective as a general-purpose algebraic solver for multiphase poromechanics, and it also scales well with problem size. The rest of the paper is organized as follows. Section 2 and 3 introduce the governing equations of multiphase poromechanics and the discretization scheme. Section 4 describes the nonlinear solution algorithm. In section 5, we describe the MGR framework and how it is applied to solve the linear systems coming from the linearization. We show numerical results in section 6 to demonstrate the robustness and scability of the proposed preconditioner. We then end with some concluding remarks and directions for future work.

2 Problem Statement

We focus on a displacement-saturation-pressure formulation of immiscible two-phase flow through a deforming poroelastic medium [13]. We limit the discussion to quasi-static small-strain kinematics. Let the subscript $w$ and $nw$ denote the wetting and non-wetting fluid phase, respectively. Since the medium’s pore space is always fluid-filled, the fluid phase saturations must always sum to unity, i.e. $(s_{w}+s_{nw})=1$. Here, the wetting fluid phase saturation, denoted from now on by lower case $s$ without subscript, is used as a primary unknown. Capillary pressure, which is the difference between the phase pressure of the non-wetting phase and the wetting phase, is not considered—a frequent assumption in many practical engineering applications. Hence, we have $p_{w}=p_{nw}=p$.

For a given closed domain $\overline{\Omega}=\Omega\cup\Gamma\in\mathbb{R}^{3}$, with $\Omega$ an open set and $\Gamma$ its boundary, and time interval $\mathbb{I}=(0,t_{\max}]$, the strong form of the multiphase poromechanical initial/boundary value problem (IBVP) consists of finding the displacement vector field $\bm{u}:\overline{\Omega}\times\mathbb{I}\rightarrow\mathbb{R}^{3}$, the wetting fluid phase saturation $s:\overline{\Omega}\times\mathbb{I}\rightarrow\mathbb{R}^{3}$, and the pore pressure $p:\overline{\Omega}\times\mathbb{I}\rightarrow\mathbb{R}$ such that [13]:

[TABLE]

where

• •

$\bm{\sigma}=\left(\textbf{C}:\nabla^{s}\bm{u}-bp\bm{1}\right)$ is the total Cauchy stress tensor, with C the rank-4 elasticity tensor, $b$ the Biot coefficient, and $\bm{1}$ the rank-2 identity tensor;

• •

$\rho\bm{g}$ is a body force due to the self-weight of the multiphase mixture, with $\rho=((1-\phi)\rho_{s}+\phi\rho_{w}s+\phi\rho_{nw}(1-s))$ the density of the mixture, $\phi$ the porosity, $\rho_{s}$, $\rho_{w}$, $\rho_{nw}$ the density of the solid, the wetting, and the non-wetting fluid phase, respectively, and $\bm{g}$ the gravity vector;

• •

$m_{w}=(\phi\rho_{w}s)$ and $m_{nw}=(\phi\rho_{nw}(1-s))$ denote wetting and non-wetting fluid phase mass per unit volume;

• •

$\bm{w}_{w}=-(\rho_{w}\lambda_{w}\bm{\kappa}\cdot\nabla\Phi_{w})$ and $\bm{w}_{nw}=-(\rho_{nw}\lambda_{nw}\bm{\kappa}\cdot\nabla\Phi_{nw})$ are wetting and non-wetting fluid phase mass fluxes [2], with $\lambda_{\ell}=k_{r{\ell}}/\mu_{\ell}$ the phase mobility, $\mu_{\ell}$ the phase viscosity, $k_{r\ell}$ the phase relative permeability factor, $\bm{\kappa}$ the absolute permeability tensor, $\Phi_{\ell}=(p-\rho_{\ell}\bm{g}\cdot\bm{x})$ the phase potential, $\bm{x}$ the position vector in $\mathbb{R}^{3}$, $\ell=\{w,nw\}$;

• •

$q_{w}$ and $q_{nw}$ are mass source/sink per unit volume terms for the wetting and the non-wetting fluid phase, respectively;

• •

$\nabla$, $\nabla^{s}$, and $\nabla\cdot$ are the gradient, symmetric gradient, and divergence operator, respectively;

• •

the superposed dot, $\dot{(\bullet)}$, indicates the derivative of quantity $(\bullet)$ with respect to time .

For the application of boundary conditions, let us introduce two disjoint partitions of the domain boundary such that $\Gamma=\overline{\Gamma_{u}^{D}\cup\Gamma_{u}^{N}}=\overline{\Gamma_{f}^{D}\cup\Gamma_{f}^{N}}$. Without loss of generality, consider homogeneous displacement boundary conditions $\bm{u}=\bm{0}$ on $\Gamma_{u}^{D}\times\mathbb{I}$ and homogeneous flux conditions $\bm{w}_{w}\cdot\bm{n}=\bm{w}_{nw}\cdot\bm{n}=0$ on $\Gamma_{f}^{N}\times\mathbb{I}$, along with a prescribed total traction conditions $\bm{\sigma}\cdot\bm{n}=\bm{t}^{N}$ on $\Gamma_{u}^{N}\times\mathbb{I}$ and pressure/saturation conditions $p=p^{D}$ and $s=s^{D}$ on $\Gamma_{f}^{D}\times\mathbb{I}$, where $\bm{n}$ denotes the outer normal vector for $\Gamma$. More complicated boundary conditions may be introduced as needed with modest changes to the discretization below. The formulation is completed by appropriate: (i) initial conditions for $\bm{u}$, $s$, and $p$, and (ii) equations of state and constitutive equations to specify the following dependencies: $\phi=\phi(\bm{u},p)$, $\rho_{\ell}=\rho_{\ell}(p)$, $\mu_{\ell}=\mu_{\ell}(p)$, and $k_{r\ell}=k_{r\ell}(s)$, with $\ell=\{w,nw\}$. For additional details on the adopted poromechanical model we refer the reader to [46].

3 Discretization

Several space discretization methods for the multiphase poromechanical IBVP have been proposed in the literature—see, e.g., [29, 22, 46, and references therein]. In this work, the discrete form of (1) is obtained by combining a finite element (FE) method for the mechanical subproblem with a finite volume (FV) approach for the multiphase flow and transport subproblem. This choice is quite common when modeling nonlinear hydromechanical processes in subsurface formations characterized by highly heterogeneous hydrogeological properties, e.g. high-contrast permeability fields typically encountered in practical reservoir simulation [37, 26, 32, 19, 38].

Let $\bm{H}_{0}^{1}(\Omega)$ denote the Sobolev space of vector functions satisfying displacement homogeneous Dirichlet conditions over $\Gamma_{u}^{D}$ and whose first derivatives belong to $L^{2}(\Omega)$, with $L^{2}(\Omega)$ the space of square integrable functions in $\Omega$; let $\bm{\mathcal{U}}^{h}\subset\mathbf{H}_{0}^{1}(\Omega)$, $\mathcal{S}^{h}\subset L^{2}(\Omega)$, $\mathcal{P}^{h}\subset L^{2}(\Omega)$ denote finite-dimensional subspaces associated with a conforming triangulation $\mathcal{T}^{h}$ of the domain into nonoverlapping hexahedral cells; and let $\widehat{w}_{\ell}^{\varepsilon}$ denote a conservative numerical flux approximating the $\ell$ fluid phase mass flux across face $\varepsilon$ in $\mathcal{E}^{h}$, namely the set of faces in $\mathcal{T}^{h}$, such that $\widehat{w}_{\ell}^{\varepsilon}\approx-\int_{\varepsilon}\bm{w}_{\ell}\cdot\bm{n}_{\varepsilon}\;\mathrm{d}A$, with $\bm{n}_{\varepsilon}$ a unit normal vector defining the global face orientation.

Precisely, our space discretization employs: (i) continuous piecewise trilinear finite elements for $\bm{\mathcal{U}}^{h}$, (ii) piecewise constant functions for $\mathcal{S}^{h}$ and $\mathcal{P}^{h}$, and (iii) a linear two-point flux approximation (TPFA) combined with a first-order upwinding strategy for the numerical flux $\widehat{w}_{\ell}^{\varepsilon}$ [2]. Using a fully-implicit time marching scheme, with the subscript $n$ indicating the discrete time level, the fully discrete mixed FE/FV variational statement of (1) is: find $\{\bm{u}^{h}_{n},s^{h}_{n},p^{h}_{n}\}\in\boldsymbol{\mathcal{U}}^{h}\times\mathcal{S}^{h}\times\mathcal{P}^{h}$ such that for all $\{\bm{\eta}^{h},\psi^{h},\chi^{h}\}\in\boldsymbol{\mathcal{U}}^{h}\times\mathcal{S}^{h}\times\mathcal{P}^{h}$

[TABLE]

where $n\in\{1,2,\ldots\}$. The compact notation $(\bullet,\bullet)$ denotes the $L^{2}$-inner product of scalar, vector, or rank-2 tensor functions in $L^{2}(\Omega)$, $[L^{2}(\Omega)]^{3}$, or $[L^{2}(\Omega)]^{3\times 3}$, as appropriate. In (2b)-(2c), $\Delta t_{n}=(t_{n}-t_{n-1})$ is the timestep size; $\mathcal{E}_{f}^{h,N}$ is the set of faces belonging to the boundary $\Gamma_{f}^{N}$; and $\llbracket\bullet\rrbracket_{\varepsilon}$ indicates the jump of a quantity $(\bullet)$ across $\varepsilon$. For an internal face $\varepsilon$ shared by cells $K$ and $L$, $\bm{n}_{\varepsilon}$ pointing from $K$ to $L$, $\llbracket\bullet\rrbracket_{\varepsilon}=((\bullet)_{|L}-(\bullet)_{|K})$, with $(\bullet)_{|K}$ and $(\bullet)_{|L}$ the restriction of $(\bullet)$ on $K$ and $L$, respectively. For a boundary lying on $\Gamma_{f}^{D}$, $\bm{n}_{\varepsilon}$ coincides with the outer normal to the domain boundary and the jump expression simply reads $\llbracket\bullet\rrbracket_{\varepsilon}=-(\bullet)_{|K}$.

Finally, introducing in (2) the expressions $\bm{u}_{n}^{h}=\sum_{i}u_{i,n}\bm{\eta}_{i}^{h}$, $s_{n}^{h}=\sum_{j}s_{j,n}\psi_{j}^{h}$, and $p_{n}^{h}=\sum_{k}p_{k,n}\chi_{k}^{h}$, with $\{\bm{\eta}_{i}^{h}\}$, $\{\psi_{j}^{h}\}$, and $\{\chi_{k}^{h}\}$ bases for $\bm{\mathcal{U}}^{h}$, $\mathcal{S}^{h}$, and $\mathcal{P}^{h}$, respectively, a standard Galerkin approach yields a system of nonlinear discrete equations

[TABLE]

Here, vector $\mathbf{x}_{n}$ contains the nodal displacement ($u_{i,n}$), cell-centered saturation $s_{i,n}$ and cell-centered pressure $p_{i,n}$ coefficients that are used to expand $\bm{u}^{h}_{n}$, $s^{h}_{n}$, and $p^{h}_{n}$ in terms of the respective basis functions at time level $n$.

4 Newton-Krylov Solver

The nonlinear system (6) is solved by means of Newton’s method, with a backtracking strategy added for robustness. The solution at time $t_{n}$ is computed as follows. Given an initial guess $\mathbf{x}_{n}^{0}$, for $k=0,1,\ldots$, until convergence

[TABLE]

where $A(\mathbf{x}_{n}^{(k)})=(\partial F/\partial\mathbf{x}_{n})(\mathbf{x}_{n}^{(k)})$ is the Jacobian matrix associated with the nonlinear residual function $F$, and $\lambda\in(0,1]$ is an appropriately chosen line-search parameter. For convenience of notation, we omit from now on to specify that $A$ is evaluated at $\mathbf{x}_{n}^{(k)}$. Clearly, at each nonlinear iteration $k$, the solution of a linear system with $A$ is required.

The linearization of Eq. 6 produces a Jacobian system with an inherent $3\times 3$ block structure

[TABLE]

This system has size proportional to the number of vertices (three displacement degrees of freedom per vertex) and cells (one saturation and one pressure degree of freedom per cell) in the computational mesh. For detailed expressions of the sub-blocks in $A$, we refer the reader to [46]. Briefly, we emphasize the properties of the three diagonal blocks that motivate choices in designing the preconditioning operator described in Section 5. Specifically:

• •

$A_{uu}$ is the elasticity block and has the structure of a discrete elliptic operator;

• •

$A_{ss}$ is the saturation block that, in the abscence of capillarity effects, has the structure of a discrete time-dependent hyperbolic problem;

• •

$A_{pp}$ is the pressure block that, similar to the elastic block, has the structure of a discrete elliptic operator.

In this work, the linear system with matrix $A$ is solved iteratively with generalized minimal residual (GMRES) [35], a Krylov subspace method designed for nonsymmetric systems. Since Krylov methods’ practical convergence depends on the availability of an effective preconditioner, we introduce the preconditioning operator $M$ and replace the linear system in Eq. 7 with the right preconditioned system,

[TABLE]

where $\Delta\mathbf{x}=M^{-1}\Delta\mathbf{y}$. In the following section, we describe an algebraic method to construct $M$ given a matrix $A$ with the structure specified in Eq. 8.

5 Multigrid Reduction

The idea of MGR has been around for a long time, tracing back to the work of Ries and Trottenberg [33, 34]. Recently, it has gained more attention through the work on multigrid reduction in time by Falgout et al. [15, 16]. MGR has also been applied successfully for problems in reservoir simulation and multiphase flow in porous media with phase transitions [8, 44]. A major advantage of the MGR approach is that it is an algebraic method and unlike geometric multigrid, it can be applied to general geometries and grid types. In this section, we first summarize the approach for the case of two-level reduction and then present the general multi-level reduction algorithm.

5.1 Two-grid Reduction Scheme

For a matrix $A$ of size $N\times N$, we define a partition of the row indices of the matrix into C-points and F-points. The C-points play a role analogous to the points on a coarse grid, and the F-points belong to the set that is the complement of the C-points. It is important to note that this partitioning is different from the one normally used in standard multigrid methods, in which the F-points correspond to all points on the fine grid, i.e. the set of F-points contains the set of C-points. In multigrid reduction, the C-points and F-points belong to non-overlapping sets. Following [15], using such CF-splitting we have

[TABLE]

where $I_{CC}$ and $I_{FF}$ are identity matrices and $S=A_{CC}-A_{CF}A_{FF}^{-1}A_{FC}$ is the Schur complement. We can define the ideal interpolation and restriction operators by

[TABLE]

Additionally, define the injection operator as $Q=\begin{pmatrix}I_{FF}\\ 0\end{pmatrix}$. Then since $A_{FF}=Q^{T}AQ$ and $S=RAP$, it is simple to derive that

[TABLE]

and

[TABLE]

where the equivalence occurs since $RAQ=Q^{T}AP=0$. This identity defines the two-level multigrid method with the ideal Petrov-Galerkin coarse-grid correction $(RAP)^{-1}$ and the F-relaxation $Q(Q^{T}AQ)^{-1}Q^{T}$: (i) Equation 13 is the additive MGR identity and (ii) Eqs. 14 and 15 are multiplicative identities with pre-smoothing and post-smoothing F-relaxation, respectively. However, constructing ideal interpolation and restriction operators is impractical. Similarly, computing the coarse-grid correction exactly is expensive, so we need to approximate these operators. In practice, MGR methods use a scalable solver such as AMG for the coarse-grid solve, and replace the ideal restriction and prolongation $R$ and $P$ with

[TABLE]

where

[TABLE]

There are many ways to construct these approximations. One simple choice is to use an injection operator for restriction and a Jacobi approach for interpolation

[TABLE]

where $D_{FF}=\text{diag}(A_{FF})$. Then the coarse grid operator $A_{h}=\tilde{R}A\tilde{P}$ can also be considered as an approximation to the Schur complement $S$. Besides the choices in Eq. 18, one can also choose to use Jacobi approach for restriction, that is $W_{r}=-A_{CF}D_{FF}^{-1}$. Another option is to construct $A_{FF}^{-1}$ using incomplete factorizations (ILU) or sparse approximate inverse techniques, such as sparse approximate inverse (SPAI) [21], factored sparse approximate inverse (FSAI) [18], or minimal residual (MR) [11]. Although these methods could provide a better approximation to $A_{FF}^{-1}$, and therefore better approximations for the restriction and interpolation operators, they tend to make these operators dense. The resulting coarse grid also becomes dense and unamenable to AMG. One can certainly apply a dropping strategy to keep such $\tilde{P}$ and $\tilde{R}$ sparse, but in practice, the potential improvement in performance using approximate inverse methods is usually offset by the cost to construct the approximation, which makes simple methods like Jacobi more appealing.

In general, we define the MGR operator with either pre-smoothing or post-smoothing F-relaxation by

[TABLE]

where $M_{CC}=(\tilde{R}A\tilde{P})$ is the coarse-grid correction and $M_{FF}^{-1}$ is the F-relaxation smoother. Additionally, similar to AMG methods, one can also apply a global smoothing step that extends to all the unknowns, not just the F-points. For the global smoother $M_{glo}^{-1}$, various methods including (block) Jacobi, (block) Gauss-Seidel, or ILU, can be used. The inclusion of this step can help eliminate error modes that both the F-relaxation and coarse-grid correction may have missed. The application of the two-grid MGR scheme consisting of a global smoother and an F-relaxation followed by a coarse-grid correction can be summarized as shown in Algorithm 1.

Balancing the quality of the approximation to the Schur-complement and the convergence of the coarse-grid solve is key to the success of MGR. One extreme is to design a coarse grid that is perfectly suitable for AMG. Assuming, for example, that the block $A_{CC}$ comes from a scalar elliptic PDE and $A_{CC}$ is SPD, then one can choose $W_{p}=W_{r}=0$ and the coarse grid becomes $RAP=A_{CC}$. In this case, the convergence of the coarse grid solve is optimal, but the approximation of the Schur-complement far from ideal, since the coarse grid neither takes into account any information from the F-points nor the coupling between the C and F points. At the other extreme, one can use the exact Schur-complement as the coarse grid by choosing $W_{r}=-A_{CF}A_{FF}^{-1}$ and $W_{p}=0$. However, because of the exact inversion of $A_{FF}^{-1}$, the coarse grid is dense. Furthermore, since the F-points and C-points actually represent equations obtained from the discretization of different continuous physical models, capturing the coupling between them on the coarse grid can lead to loss of ellipticity, which can make the coarse-grid solve with AMG ineffective. Thus, finding a good approximation of the Schur-complement that is still amenable to AMG methods is essential.

Remark 5.1**.**

*The appeal of the MGR approach is that it provides a general framework for choosing the coarse/fine grids, the interpolation and restriction operators, and the solvers for the F-relaxation and coarse-grid correction. As an example, it was shown in [8, 44] that one can recast any CPR-AMG strategy [9, 20, 27, 30, 36, 40, 48] or block preconditioner [7] used in reservoir simulation as a particular variant of the two-grid MGR reduction scheme by appropriately defining the different components of the algorithm, namely prolongation, restriction and smoothing operators. *

5.2 A general multi-level MGR algorithm

One can replace the coarse grid solve in Algorithm 1 with a two-level MGR scheme and apply the method recursively to obtain a multi-level MGR algorithm. The general application of the MGR V-cycle with global smoothing is summarized in Algorithm 2, where the hierarchy of coarse grid operators, i.e. $A_{l+1}=\tilde{R}_{l}A_{l}\tilde{P}_{l}$, is assumed to be computed for each level $l$.

Based on Algorithm 2, W- and F-cycle versions of the MGR algorithm can also be defined [34]. Note that the Schur-complement $S$ is approximated by the triple product $RAP=A_{CC}-A_{CF}D_{FF}^{-1}A_{FC}$ in the classical two-grid reduction scheme in Algorithm 1. Even though we have introduced a sparse approximation to $S$ by replacing $A_{FF}$ with its diagonal $D_{FF}$, i.e. $A_{FF}^{-1}\approx D_{FF}^{-1}$, in a multi-level reduction scheme, the coarse grid can still become dense or unsuitable for standard AMG because the correction term $A_{cor}=A_{CF}D_{FF}^{-1}A_{FC}$ involves a matrix-matrix product.

In this work, we develop a dropping strategy for $A_{cor}$ to keep the coarse grid sparse as well as suitable for AMG. One approach is to drop all entries of $A_{cor}$ that are smaller than a prescribed tolerance. Here, we use a different strategy based on a maximum number of non-zero values per row. Specifically, we choose to keep only $N_{\text{max}}$ entries with largest absolute values on each row. To preserve at least some information of the first level of reduction, however, we always keep the diagonals of the sub-blocks in $A_{cor}$. For instance, in a three-level reduction scheme, in the first-level reduction, $A_{cor}$ has $2\times 2$ block structure, and applying maximum dropping (i.e. using an extremely large tolerance or $N_{\text{max}}=0$), $A_{cor}$ is still a $2\times 2$ block matrix, whose sub-blocks are diagonal matrices. Using this dropping strategy results in a non-Galerkin coarse grid

[TABLE]

where $\mathbf{G}$ is a sparsifying operator that performs one of the aforementioned dropping strategies. So far we only assume that a CF-splitting of the rows is given. How to choose such a splitting is dependent on the problem and it is up to the user to make the decision. However, as a general principle, it is usually a good idea to choose a CF-splitting so that the final coarse grid corresponds to the variable associated with an elliptic equation, e.g. pressure, since we want to solve the coarse grid using an efficient method such as standard AMG. In the next section, we show how to choose an appropriate CF-splitting at each level of reduction for our multiphase poromechanical problem.

5.3 MGR for Multiphase Poromechanics

We propose a three-level MGR reduction scheme to precondition the Jacobian matrix (8). For the first level of reduction, we aim at decoupling the mechanics sub-problem from the flow. Therefore, we assign all the displacement unknowns as F-points while both saturation and pressure unknowns are labeled as C-points. This leads to the following partitioning

[TABLE]

Then $A_{FF}\equiv A_{uu}$ and the coarse grid $A_{CC}$, which corresponds to the flow sub-problem, has the $2\times 2$ block structure

[TABLE]

For the F-relaxation step, we need to solve the elasticity problem involving the elliptic operator $A_{uu}$. Here, we use one AMG V-cycle. Because of the vectorial nature of the elasticity operator, this is the most expensive part of the setup phase. Also, given the global system size, we ignore the first-level global relaxation step. Using the interpolation and restriction operators specified in Eq. 18 combined with the dropping strategy defined in Equation 21 yields the following first level coarse grid

[TABLE]

For our multiphase poromechanics problem, we use $N_{\text{max}}=4$ for $\mathbf{G}$. Again, we emphasize the flexibility of our framework as it allows for experimenting with different choices of $\mathbf{G}$. For example, choosing an appropriate $\mathbf{G}$, we can mimic the fixed-stress preconditioner developed in [46].

The second reduction step is essentially a CPR approach that is embedded within a multigrid reduction framework. Hence, we label saturation unknowns as F-points and pressure unknowns as C-points in $S_{1}$:

[TABLE]

Again, using interpolation and restriction operators in Eq. 18, we obtain the second-level Schur-complement

[TABLE]

where $D_{ss}=\text{diag}(\tilde{A}_{ss})$. In other physics-based approaches commonly used in reservoir simulation, one can seek to further sparsify the Schur-complement. For example, in a Quasi-IMPES reduction scheme, the block $\tilde{A}_{ps}$ is also replaced by its diagonal $D_{ps}=\text{diag}(\tilde{A}_{ps})$. This approximation ensures that the matrix sparsity pattern of $A_{pp}$ coming from the original finite volume stencil is preserved and the resulting Schur-complement is near elliptic. In the MGR approach, however, no further sparse approximation is needed for this level since the flow part is relatively small compared to the elasticity block and the coarse grid generated in Eq. 36 is still well-suited for AMG. At the second level, the F-relaxation involving the $\tilde{A}_{ss}$ block is done using a simple Jacobi relaxation. However, the second-level global smoothing step is required to reduce the error associated with the hyperbolic component of the flow subproblem. Indeed, the global smoothing plays a key role particularly in the later stage of the simulation when the pressure field approaches steady-state conditions and the multiphase flow and transport process transitions to an advection dominated regime. The need for a robust global smoother will become clear through numerical results presented in the next section.

Remark 5.2**.**

*Even though we formally present the multigrid reduction framework for the $3\times 3$ system in field-ordered form, in our implementation, the input matrix has interleaved ordering for saturation and pressure. This choice produces a sparsity pattern in which dense $2\times 2$ blocks appear for the first-level coarse grid in Eq. 30. Block versions of relaxation or incomplete factorization preconditioners are therefore appealing, as dense multiplication and inversion operations can be applied to the small blocks. *

Remark 5.3**.**

*Common strategies for the second-level global smoothing step include block relaxation methods (e.g. Jacobi, Gauss-Seidel) or incomplete factorizations (e.g. ILU(k), ILUT). In this work, we explore two options. The first option uses several sweeps of hybrid block Gauss-Seidel (HBGS). The second option uses one sweep of processor-local, pointwise ILU(k) [10]. *

6 Numerical Results

We perform numerical experiments to test the performance of the MGR preconditioner on two problems: (1) a weak scaling study for a simple synthetic configuration; and (2) a strong scaling study using a realistic, highly heterogeneous reservoir based on the SPE10 [12] example. Both examples have been designed as community benchmark problems and exhibit tight-coupling between displacement, pressure, and saturation fields. Problem specifications are described in detail in [46] and so are only briefly reported below.

In this study, the simulator is provided by Geocentric, which utilizes the deal.ii Finite Element Library [4] for discretization functionality. It also provides a direct interface with MGR, which is implemented as a separate solver in hypre [17]. All the numerical experiments were run on Quartz, a cluster at the Lawrence Livermore Computing Center with 1344 nodes containing two Intel Xeon E5-2695 18-core processors sharing 128 GiB of memory on each node with Intel Omni-Path interconnects between nodes. We use pure MPI-based parallelism.

For the elastic block $A_{uu}$, we use one V-cycle of BoomerAMG [24], with an unknown approach for a system of three PDEs, with one level of aggressive coarsening, one sweep of hybrid forward $l_{1}$-Gauss-Seidel [3] for the down cycle and one sweep of hybrid backward $l_{1}$-Gauss-Seidel for the up cycle. The coarsest grid is solved directly with Gaussian elimination. The MGR coarse-solve in Eq. 36 also uses BoomerAMG with the same smoother configuration, but for a scalar problem and a Hybrid Modified Independent Set (HMIS) coarsening strategy [14]. For the global smoother, we use one step of processor-local, pointwise ILU(1).

6.1 Staircase Benchmark

The configuration of the first test problem is illustrated in Figure 1. A highly-permeable channel winds its way in a “staircase” fashion through a lower-permeability host rock. A denser, wetting phase is injected through a well at the top corner, leading to a saturation plume driven by gravity and pressure that migrates towards a production well in the lower corner. The whole system is deformable and exhibits significant poromechanical coupling. Visualizations of the resulting pressure and deformation fields have been omitted for brevity. A detailed specification of mesh geometry, material properties, and boundary conditions can be found in [46].

Table 1 shows the results for a weak scaling study using the staircase example. We keep the number of degrees of freedom per core constant at 44k and increase the number of cores from 2 to 1024. The global problem size grows $8^{3}$ times from 88k to 42M. Due to the inherent nonlinearity, we observe an increase in the number of Newton iterations per time step as the mesh is refined. The average number of GMRES iterations per Newton step, however, only experiences a modest growth as desired. Even though we can use a more complex smoother in place of the hybrid $l_{1}$-Gauss-Seidel solves for the elasticity block and drive down the number of iterations, that will come at the expense of run-time performance. In general, we find that the $l_{1}$-Gauss-Seidel smoother strikes a good balance between iteration counts and run time. Similar to the number of iterations, the total run time, including both the setup and solve phases, also exhibits some growth, but again, the result is quite satisfactory even for large core counts. The increase in the run time can be attributed to communication costs in the MGR setup and solve phases, since the actual number of degrees of freedom per core is fairly small. Overall, however, the MGR framework provides a good platform for scalable performance.

6.2 SPE10-based Benchmark

We also perform a strong scaling study on a more realistic benchmark problem derived from the second model of the SPE10 Comparative Solution Project (Figure 2) [12]. The original SPE10 permeability and porosity fields are now treated as a poromechanical medium. These geostatistically generated fields exhibit both severe heterogeneity and anisotropy. In the current poromechanical benchmark, the reservoir itself is also embedded in a larger domain—with caprock and underburden—to provide more realistic boundary conditions. Water is injected through a central well, while fluids are produced from four wells at the corners of the domain. Mesh, material property, and boundary condition specifications are reported in [46]. Note that the well control conditions differ from the original SPE10 model to avoid well impacts on the linear solver. The treatment of well degrees-of-freedom within the linear solver is a critical issue, but is deliberately left out-of-scope for the current contribution. We remark, however, that the MGR approach provides a flexible framework to treat this additional complexity.

The resulting discrete problem has 16.7M degrees-of-freedom. We keep the problem size fixed and divide the work across an increasing number of compute cores. The results are shown in Table 2. Again, we observe only minor growth in the number of GMRES iterations with larger core counts. Similar to the weak scaling case, the reason for this growth is the use a hybrid $l_{1}$-Gauss-Seidel smoother in AMG solves for the elasticity block and the coarse grid. Good overall timing efficiency is also achieved up to 288 cores. For 576 cores, even though we still get good efficiency for the solve phase, there is a noticeable increase in the setup time because the problem size on each core becomes very small, i.e. about 17k total and less than 6k degrees of freedom for the elastic block and the coarse grid, respectively. Consequently, the majority of the time is spent in communication while not much computation is performed. However, the results still indicate that one can use the proposed framework with a large number of processors to efficiently reduce the long simulation time for challenging problems with highly heterogeneous media.

6.3 Effect of global smoother

As we have mentioned earlier, the performance of MGR is dependent on the effectiveness of the solvers for each component of the algorithm. In general, changing the configuration for one component, e.g. smoother choices for the F-relaxation or coarsening strategies for the coarse-grid AMG solve, would result in a different number of GMRES iterations. However, the effect could also be quite subtle and not manifest itself until the underlying property of the problem changes. For multiphase poromechanics simulations, early times are typically dominated by elliptic effects associated with the pressure and displacement fields, while at late times the hyperbolic effects associated with the saturation field become significant.

Here we explore the effectiveness of different global smoothers on the multiphase flow system as the simulation progresses. The first option uses three sweeps of HBGS, and the second option uses a single sweep of processor-local ILU(1). As one can see from Fig. 3(a), there is no apparent difference between the two smoothers until about 45 days of injection, when the number of iterations for the HBGS method increases sharply and continues to stay high. In contrast, ILU(1) is less sensitive. Even though the number of iterations also rises slightly around 85 days, it starts to decrease for the last period at the end of the simulation. We also plot the total time of HBGS(3) against ILU(1) in Fig. 3(b). It is clear that even though ILU(1) takes slightly more time in the beginning (mainly due to higher setup cost), the trade off is worthwhile thanks to its robustness, which leads to a modest reduction in total time for the whole simulation. This observation is confirmed by the widespread use of incomplete factorization smoothers in the reservoir simulation community.

7 Conclusion

In this work, we have presented an algebraic framework based on multigrid reduction for solving the linear system that comes from discretizing and linearizing the conservation equations governing multiphase flow coupled with poromechanics. This framework is flexible and allows us to construct different preconditioners based on different choices for CF-splitting strategies, interpolation and restriction operators, as well as solvers and smoothers. We have also developed a dropping strategy for computing the reduction onto the coarse-grid within the MGR V-cycle that captures the coupling between mechanics and flow and reduces the operator complexity at the same time. This results in an algebraic preconditioner that is robust and scalable for realistic and large-scale simulation of multiphase poromechanics.

Regarding future work, a number of improvements to the MGR framework could be explored. For example, constructing good approximations to the ideal interpolation and restriction operators that have low complexity remains a significant challenge. Also, it is unclear how one can choose an optimal coarse grid that is representative of the fine-grid system and at the same time still amenable to AMG in a multi-level reduction setting. Thus, better strategies for computing the non-Galerkin coarse grid introduced in this work are needed to improve robustness of the framework. Lastly, since MGR is designed to accommodate a wide range of coupled systems, we are looking into extending the approach to solve problems with non-isothermal flow and fractured media.