A Two-Stage Preconditioner for Multiphase Poromechanics in Reservoir Simulation

TL;DR

This paper introduces a two-stage Newton-Krylov solver with a novel block-preconditioning strategy that enhances scalability, robustness, and convergence for large-scale multiphase poromechanics simulations in reservoir modeling.

Contribution

The paper presents a new two-stage solution algorithm with a block-preconditioning approach that significantly improves the efficiency and scalability of multiphase poromechanics simulations.

Findings

Rapid convergence demonstrated on benchmark problems

Good parallel scalability achieved

Robust performance with heterogeneous materials

Abstract

Many applications involving porous media--notably reservoir engineering and geologic applications--involve tight coupling between multiphase fluid flow, transport, and poromechanical deformation. While numerical models for these processes have become commonplace in research and industry, the poor scalability of existing solution algorithms has limited the size and resolution of models that may be practically solved. In this work, we propose a two-stage Newton-Krylov solution algorithm to address this shortfall. The proposed solver exhibits rapid convergence, good parallel scalability, and is robust in the presence of highly heterogeneous material properties. The key to success of the solver is a block-preconditioning strategy that breaks the fully-coupled system of mass and momentum balance equations into simpler sub-problems that may be readily addressed using targeted algebraic methods. Numerical results are presented to illustrate the performance of the solver on challenging benchmark problems.

Figures

Tables

Table 1. Table 1 : Parameter values used for the numerical examples. Note that the units presented below are chosen for clarity but do not form a self-consistent unit system. They are converted to consistent values prior to use within the simulator.

Parameter	Units	Staircase Example	SPE10-based Example
Porosity:
Reservoir	–	0.20	SPE10 data ChrBlu01
Seal units	–	0.05	0.01
Permeability:
Reservoir	mD	1000	SPE10 data ChrBlu01
Seal units	mD	1	0.01
Relative perm:
Residual water sat.	–	0.2	0.2
Residual oil sat.	–	0.2	0.2
Water:
Reference density	kg/m3	1035	1035
Compressibility	MPa-1	$4.34\cdot {10}^{-4}$	$4.34\cdot {10}^{-4}$
Viscosity	cP	0.3	0.3
Oil:
Reference density	kg/m3	863	863
Compressibility	MPa-1	$1.98\cdot {10}^{-4}$	$1.98\cdot {10}^{-4}$
Viscosity	cP	3.0	3.0
Rock:
Young’s modulus	MPa	5000	5000
Biot coefficient	–	1	1
Grain density	kg/m3	2650	2650
Well control:
Injection $\Delta$BHP	MPa	$5$	$5$
Production $\Delta$BHP	MPa	$-5$	$-5$
Ramp time	day	1	1
Well radius	m	0.1524	0.1524
Skin factor	–	0	0
Time-stepping:
Initial $\Delta t$	day	0.1	0.1
Maximum $\Delta t$	day	1	5
End time	day	100	100
Solver tolerances:
Newton	–	10-5	10-5
Krylov	–	10-6	10-6

Table 2. Table 3 : Strong scaling performance for the SPE10-based example. The number of GMRES iterations is well-controlled and mostly insensitive to the parallel partitioning.

Cores	DoF / Core	Iterations			Setup Phase [s]		Solve Phase [s]	Overall [s]	Efficiency
		Newton per	GMRES per		Mechanics	Flow
		Timestep	Newton
		(avg.)	(avg.)		(once)	(avg.)	(avg.)
36	464,738	6.5	24.0		3.27	1.27	36.0	12,634	1.00
72	232,369	6.5	24.7		1.71	0.59	18.0	6,303	1.00
144	116,184	6.5	25.6		0.99	0.31	9.0	3,149	1.00
288	58,092	6.5	26.0		1.36	0.18	4.5	1,603	0.98
576	29,046	6.5	28.1		0.86	0.40	2.2	881	0.90

Equations

$$\mathbf{u}:\overline{\Omega }\times \mathbb{I}\to {\mathbb{R}}^3$$

$$\nabla \cdot \mathbf{\sigma }+\rho \mathbf{g}=0$$

$$\frac{\partial m_w}{\partial t}+\nabla \cdot {\mathbf{w}}_w=q_w^I-q_w^P$$

$$\frac{\partial m_{nw}}{\partial t}+\nabla \cdot {\mathbf{w}}_{nw}=q_{nw}^I-q_{nw}^P$$

$$\mathbf{u}=\bar{\mathbf{u}}$$

$$\mathbf{\sigma }\cdot {\mathbf{n}}_{\Gamma }=\bar{\mathbf{t}}$$

$$p=\bar{p}$$

$$s=\bar{s}$$

$${\mathbf{w}}_w\cdot {\mathbf{n}}_{\Gamma }={\bar{w}}_w$$

$${\mathbf{w}}_{nw}\cdot {\mathbf{n}}_{\Gamma }={\bar{w}}_{nw}$$

$$\mathbf{u}(\mathbf{x},0)={\mathbf{u}}_0(\mathbf{x})$$

$$s(\mathbf{x},0)=s_0(\mathbf{x})$$

$$p(\mathbf{x},0)=p_0(\mathbf{x})$$

$$\mathbf{\sigma }={\mathbb{C}}_{\text{dr}}:{\nabla }^s\mathbf{u}-bp1,$$

$${\mathbf{v}}_w=-{\lambda }_w\mathbf{\kappa }\cdot \nabla (p+{\rho }_{w}gz),$$

$$\Delta \varphi$$

$$\{{\mathbf{u}}_n^h,s_n^h,p_n^h\}\in {\mathcal{U}}^h\times {\mathcal{S}}^h\times {\mathcal{P}}^h$$

$${\mathcal{R}}_u=$$

$${\mathcal{R}}_s=$$

$${\mathcal{R}}_p=$$

$$\{\begin{array}{l}{\mathcal{U}}^h:=\left\{\mathbf{\eta }|\mathbf{\eta }\in [C^0(\overline{\Omega }){]}^3,\mathbf{\eta }=0\text{ on }{\Gamma }_u^{D,h},{\mathbf{\eta }}_{|K}\in [{\mathbb{Q}}_1(K){]}^3\forall K\in {\mathcal{T}}^h\right\},\\ {\mathcal{S}}^h:=\left\{\psi |\psi \in L^2(\Omega ),{\psi }_{|K}\in {\mathbb{P}}_0(K)\forall K\in {\mathcal{T}}^h\right\},\\ {\mathcal{P}}^h:=\left\{\chi |\chi \in L^2(\Omega ),{\chi }_{|K}\in {\mathbb{P}}_0(K)\forall K\in {\mathcal{T}}^h\right\},\end{array}$$

$${\mathbf{u}}_n^h(\mathbf{x})$$

$$F_{\ell ,n}^f$$

$${\varrho }_{\ell ,n}^f$$

$$q_{\ell ,n}^P$$

$$q_{\ell ,n}^I$$

$${\Phi }_{\ell ,n}^{K,W}$$

$$WI$$

$$\{\begin{array}{ll}{\mathsf{r}}_n^u({\mathsf{x}}_n)& =0,\\ {\mathsf{r}}_n^s({\mathsf{x}}_n,{\mathsf{x}}_{n-1})& =0,\\ {\mathsf{r}}_n^p({\mathsf{x}}_n,{\mathsf{x}}_{n-1})& =0,\end{array}$$

$${\mathsf{A}}^k\Delta \mathsf{x}=-{\mathsf{r}}^k,$$

$${\mathsf{x}}_n^k={\mathsf{x}}_n^{k-1}+\Delta \mathsf{x},$$

$${\left[\begin{array}{ccc}{\mathsf{A}}_{uu}& {\mathsf{A}}_{us}& {\mathsf{A}}_{up}\\ {\mathsf{A}}_{su}& {\mathsf{A}}_{ss}& {\mathsf{A}}_{sp}\\ {\mathsf{A}}_{pu}& {\mathsf{A}}_{ps}& {\mathsf{A}}_{pp}\end{array}\right]}^k\left[\begin{array}{c}\Delta \mathsf{u}\\ \Delta \mathsf{s}\\ \Delta \mathsf{p}\end{array}\right]$$

$${\left(\mathsf{A}{\mathsf{M}}^{-1}\right)}^k\Delta \mathsf{y}=-{\mathsf{r}}^k,$$

$$\Delta \mathsf{x}={\mathsf{M}}^{-1}\Delta \mathsf{y}.$$

$${\eta }_k=\gamma {\left(\frac{\Vert {\mathsf{r}}^k\Vert }{\Vert {\mathsf{r}}^{k-1}\Vert }\right)}^{\omega },$$

$${\eta }_k\leftarrow \max [{\eta }_k,\gamma {\eta }_{k-1}^{\omega }]\text{whenever}\gamma {\eta }_{k-1}^{\omega }>\mathrm{0.1.}$$

$$\displaystyle\boldsymbol{\mathsf{A}}=\left[\begin{array}[]{c:cc}\mathsf{A}_{\mathit{uu}}&\mathsf{A}_{\mathit{us}}&\mathsf{A}_{\mathit{up}}\\ \hdashline\mathsf{A}_{\mathit{su}}&\mathsf{A}_{\mathit{ss}}&\mathsf{A}_{\mathit{sp}}\\ \mathsf{A}_{\mathit{pu}}&\mathsf{A}_{\mathit{ps}}&\mathsf{A}_{\mathit{pp}}\end{array}\right]=\begin{bmatrix}\mathsf{A}_{\mathit{uu}}&\mathsf{A}_{\mathit{uf}}\\ \mathsf{A}_{\mathit{fu}}&\mathsf{A}_{\mathit{ff}}\end{bmatrix}.$$

$$\mathsf{A}$$

$${\mathsf{L}}^{-1}=\left[\begin{array}{cc}{\mathsf{A}}_{uu}^{-1}& \\ -{\mathsf{S}}_{ff}^{-1}{\mathsf{A}}_{fu}{\mathsf{A}}_{uu}^{-1}& {\mathsf{S}}_{ff}^{-1}\end{array}\right].$$

$${\mathsf{S}}_{ff}=\left[\begin{array}{cc}{\mathsf{A}}_{ss}& {\mathsf{A}}_{sp}\\ {\mathsf{A}}_{ps}& {\mathsf{A}}_{pp}\end{array}\right]-\left[\begin{array}{cc}{\mathsf{A}}_{su}{\mathsf{A}}_{uu}^{-1}{\mathsf{A}}_{us}& {\mathsf{A}}_{su}{\mathsf{A}}_{uu}^{-1}{\mathsf{A}}_{up}\\ {\mathsf{A}}_{pu}{\mathsf{A}}_{uu}^{-1}{\mathsf{A}}_{us}& {\mathsf{A}}_{pu}{\mathsf{A}}_{uu}^{-1}{\mathsf{A}}_{up}\end{array}\right].$$

Full text

††journal: arXiv preprint

\biboptions

sort&compress

\cortext

[cor1]Corresponding author

1 Introduction

Coupled simulations of multiphase fluid flow, transport, and geomechanics are useful for predicting performance of many subsurface systems. These simulations solve a set of nonlinear, time-dependent partial differential equations (PDEs) enforcing mass conservation of each fluid phase and linear momentum balance for the mixture. For this system, fully-implicit time discretization is widely preferred for its unconditional stability, allowing large timesteps. The major computational bottleneck, however, is that one must solve a sequence of large, ill-conditioned linear systems to advance the solution. Robust and scalable solvers are therefore needed to enable large-scale simulations on high performance computing platforms. In this paper, we propose an efficient preconditioning strategy for multiphase poromechanics that addresses the inherent ill-conditioning of the discrete system using a mixture of physics-based insight and flexible algebraic methods.

Over the last four decades, a rich literature focusing on the development of efficient iterative solvers for fully-implicit simulation of complex multiphase flow and transport—without mechanics—has developed within the reservoir simulation community. The most popular approach is based on Krylov subspace methods combined with a Constrained Pressure Residual (CPR) multistage preconditioning technique Wal83 ; WalKenLit85 . In essence, the CPR preconditioner uses separate stages to damp error modes associated with (i) the nearly-hyperbolic behavior of the saturation (or composition) variables and (ii) the nearly-elliptic character of the pressure variables. The standard reservoir simulator implementation LacVasWhe01 ; SchMasWen03 ; Cao_etal05 ; Stu_etal07 ; ZhoJiaThc13 ; Gri_etal14 ; LiuWanChe16 couples Incomplete LU (ILU) factorizations, effective for the hyperbolic part, with an algebraic multigrid (AMG) method, well-suited for the elliptic part. A variant of CPR preconditioning that replaces AMG by a multiscale finite volume (MSFV) or finite element (MSFE) solver was proposed in Cus_etal15 . Using a multigrid reduction (MGR) approach RieTroWin83 , the CPR method has been generalized in a broader multigrid framework and successfully applied to multiphase flow and transport with phase transitions in Wan_etal17 ; BuiWanOse18 . Recent solution alternatives to CPR can be found in Yan_etal18 ; SinPenWhe18 .

To date, the design of fully-implicit solvers for poromechanics has primarily targeted single-phase flow in fully or partially saturated media. Most efforts can be grouped in two broad categories: (i) preconditioned Krylov solvers BerFerGam07 ; WhiBor11 ; BerMar12 ; HagOsnLan12a ; WhiCasTch16 ; LeeMarWin17 ; Adl_etal18 , and (ii) multigrid solvers Gas_etal04 ; Luo_etal15 ; GasRod17 ; Luo_etal17 —with obvious overlap. In particular, the analysis of block preconditioners for the displacement-velocity-pressure formulation has been given increasing attention during recent years Lip02 ; Kuz_etal03 ; FerCasGam10 ; TurArb14 ; CasWhiFer16 . The extension to the case of multiphase flow and transport has typically been addressed by adopting sequential-implicit approaches SetMou98 ; ManLon02 ; Dai_etal02 ; Jea_etal_07 ; KimTchJua11a ; MikWhe13 ; KimTchJua13 ; DosNor15 ; CorMur18 where easier sub-problems are solved separately but iterated until coupled convergence. In this approach, a proper splitting strategy is necessary to guarantee unconditional convergence of the sequential iteration. The most established algorithms rely on the fixed-stress coupling scheme KimTchJua11a ; KimTchJua11b ; MikWhe13 .

In the present work, a preconditioning scheme for fully-implicit multiphase fluid flow, transport, and geomechanics is proposed. The algorithm relies on successive sparse Schur-complement approximations to reduce the coupled system into easier sub-problems for which algebraic methods can be applied. We build on fixed-stress partitioning SetMou98 ; KimTchJua11a ; MikWhe13 ; WhiCasTch16 and CPR (Wal83, ) concepts to construct a novel multi-stage preconditioning scheme well-suited to the problem at hand. The end result is a modular, easy-to-implement, and highly scalable solver.

The paper is organized as follows. Sections 2 and 3 introduce the mathematical model and discretization scheme. The nonlinear and linear solution algorithms are described in detail in Sections 4 and 5. Weak and strong scaling results are presented in Section 6 to demonstrate performance and robustness of the proposed preconditioner in a parallel setting. We then end with a few concluding remarks regarding future work.

2 Governing Equations

As a simple but representative model system, we consider a displacement-saturation-pressure formulation of multiphase poroelasticity (Cou04, ) assuming isothermal, quasi-static, two-phase flow of immiscible fluids. We limit the discussion to small-strain kinematics. Let $\overline{\Omega}=\Omega\cup\Gamma$ denote a closed domain in $\mathbb{R}^{3}$ occupied by a poroelastic medium, with $\Omega$ and $\Gamma$ an open set and its boundary, respectively. For the application of boundary conditions, let $\Gamma$ be partitioned as $\Gamma=\overline{\Gamma_{u}^{D}\cup\Gamma_{u}^{N}}=\overline{\Gamma_{f}^{D}\cup\Gamma_{f}^{N}}$, with $\Gamma_{u}^{D}\cap\Gamma_{u}^{N}=\Gamma_{f}^{D}\cap\Gamma_{f}^{N}=\emptyset$. Here, subscripts $u$ and $f$ denote solid displacement and fluid flow boundary conditions, respectively, while superscripts $D$ and $N$ indicate Dirichlet or Neumann conditions. Let $\mathbf{n}_{\Gamma}$ denote the outer normal vector. Wetting and non-wetting fluid phases are identified by subscript $w$ and $nw$, respectively. The requisite set of governing equations consists of one conservation law for linear momentum and two conservation laws for the mass balance of each fluid phase. Let $\mathbf{u}$, $s_{\ell}$, and $p_{\ell}$ denote the solid displacement vector, saturation, and pressure of fluid phase $\ell=\{w,nw\}$. Since the two fluids fill the voids, the saturations must satisfy the closure condition $\sum_{\ell=\{w,nw\}}s_{\ell}=1$. In the following, the wetting fluid phase saturation is selected as a primary unknown and will be denoted by lower case $s$ without subscript. Capillarity effects are not considered, leading to the simplification $p_{w}=p_{nw}=p$. This is a frequent assumption in many reservoir engineering applications, but we emphasize that it can have a significant impact on the nature of the governing equations and solver strategy.

Given a finite time interval $\mathbb{I}=(0,t_{\max}]$ the strong form of the initial/boundary value problem may be stated as follows:

[TABLE]

In Eq. (1b), the Cauchy total stress tensor is given by

[TABLE]

with $\mathbb{C}_{\text{dr}}$ the fourth-order drained elasticity tensor, $b$ the Biot coefficient, and $\boldsymbol{1}$ the second-order unit tensor; $\rho\mathbf{g}$ is a body force due to the self-weight of the multiphase mixture, with $\mathbf{g}$ the gravity vector; $\rho=[(1-\phi)\rho_{s}+\phi\rho_{w}s+\phi\rho_{nw}(1-s)]$ is the density of the mixture, with $\rho_{s}$, $\rho_{w}$, and $\rho_{nw}$ the density of the solid, the wetting, and the non-wetting fluid phase, respectively; $\phi$ is the porosity; and $\nabla^{s}$ and $\nabla\cdot$ are the symmetric gradient and divergence operator, respectively.

In Eqs. (1c), $m_{w}=(\phi\rho_{w}s)$, $\mathbf{w}_{w}=(\rho_{w}\mathbf{v}_{w})$, and $q^{\alpha}_{w}\geq 0$, $\alpha=\{I,P\}$, denote mass per unit volume, mass flux, and mass injection ($I$) and production ($P$) sources per unit volume for the wetting phase. Here, $\mathbf{v}_{w}$ is the volumetric flux. It follows the traditional multiphase flow extension of Darcy’s law MusMer36 ,

[TABLE]

with $\lambda_{w}=k_{rw}/\mu_{w}$ the phase mobility, $\mu_{w}$ the viscosity, $k_{rw}$ the relative permeability factor, $\boldsymbol{\kappa}$ the absolute permeability tensor, $g$ the gravitational acceleration, $z$ the elevation above a datum, and $\nabla$ the gradient operator. In Eqs. (1d), $m_{nw}=(\phi\rho_{nw}(1-s))$, $\mathbf{w}_{nw}=(\rho_{nw}\mathbf{v}_{nw})$, $q^{I}_{nw}\geq 0$, and $q^{P}_{nw}\geq 0$ are the counterpart of $m_{w}$, $\mathbf{w}_{w}$, $q^{I}_{w}$, and $q^{P}_{w}$ for the non-wetting fluid phase.

The formulation is completed with constitutive equations and equations of state to specify the following dependencies: $\phi=\phi(\mathbf{u},p)$, $\rho_{s}=\rho_{s}(p)$, $\rho_{\ell}=\rho_{\ell}(p)$, $\mu_{\ell}=\mu_{\ell}(p)$, $k_{r\ell}=k_{r\ell}(s)$, $\ell=\{w,nw\}$. The porosity change from a reference porosity $\phi_{0}$ is modeled in incremental form as Cou04

[TABLE]

where $\epsilon_{v}=\text{trace}(\nabla^{s}\mathbf{u})$ is the volumetric strain, and $K_{\mathit{dr}}$ is the drained (skeleton) bulk modulus. Specific forms for the remaining relative permeability and fluid compressibility relationships used in our simulations are detailed with the numerical examples in Section 6. The absolute permeability $\boldsymbol{\kappa}$ is here modeled as a constant field and exhibits no stress-dependence, though this additional form of solid-fluid coupling is of frequent interest in applications and often motivates the desire for a fully-coupled solver in the first place.

3 Discrete Formulation

A variety of choices are available for the spatial discretization of the PDEs governing multiphase poromechanics. Continuous Finite Element (FE) based formulations—the most popular technique employed in consolidation modeling, e.g. MurLou94 ; LewSch98 ; WhiBor08 ; Rod_etal16 —have been traditionally used for solving Biot’s equations. However, a discontinuous interpolation for the pressure field may be more appropriate in the presence of highly heterogeneous permeability fields. In addition, element-wise mass conservation can also be a fundamental requirement, in particular in geoscience applications. Therefore, formulations that combine FE approaches for mechanics with a Finite Volume (FV) method for fluid flow and transport are often adopted by modelers SetMou98 ; KimTchJua11a ; Pre14 ; GarKarTch16 . Some authors have recently started investigating face-centered and cell-centered finite-volume methods also for the mechanical subproblem DiP_etal11 ; Nor14 ; BofBotDiP16 ; Nor16 ; KeiNor17 to account for mechanical effects when flow processes are of predominant importance. Also, in the context of single-phase poromechanics, locally mass-conservative approaches have been proposed using mixed three-field (displacement-velocity-pressure) or four-field (stress tensor-displacement-velocity-pressure) formulations, respectively (e.g., Fri00, ; Lip02, ; PhiWhe07a, ; PhiWhe07b, ; JhaJua07, ; FerCasGam10, ; HagOsnLan12b, ; TurArb14, ; Yi14, ; Bot_etal17, ; BauRadKoc17, ; Hu_etal17, ; Rod_etal17, ; HonKrau18, ; Bot_etal18, ).

We partition the domain using a computational mesh $\mathcal{T}^{h}$. The mesh consists of disjoint, nonoverlapping cells such that $\Omega\approx\Omega^{h}=\bigcup_{i}K_{i}$ with $K_{i}\in\mathcal{T}^{h}$. A unique direction is associated with each mesh face $f$ through a unit normal vector $\mathbf{n}_{f}$. For an internal face we let $\mathbf{n}_{f}$ point in the direction of increasing global cell indices—i.e. for $f$ shared by cells $K$ and $L$, with $L$ greater than $K$, $\mathbf{n}_{f}=\mathbf{n}_{K}$ where $\mathbf{n}_{K}$ is the outward normal unit vector to cell $K$’s boundary. For a boundary interface we assume $\mathbf{n}_{f}=\mathbf{n}_{\Gamma}^{h}$. In the case of a non-planar face, a plane characterized by a mean normal $\mathbf{n}_{f}$ is associated to $f$.

The numerical solution of the system of governing equations (1) relies on a mixed FE-FV spatial discretization and a fully-implicit time-marching scheme (backward Euler) Hug00 ; EymGalHer00 ; Qua14 . Without loss of generality, we assume homogeneous prescribed displacement ($\bar{\mathbf{u}}=\mathbf{0}$) and wetting and non-wetting fluid phase mass fluxes ($\bar{w}_{w}=\bar{w}_{nw}=0$). Time integration is carried out by partitioning the time interval $\mathbb{I}$ into $n_{\Delta t}$ subintervals $\mathbb{I}_{n}=(t_{n-1},t_{n}]$, $n=1,\ldots,n_{\Delta t}$, with $\Delta t_{n}=(t_{n}-t_{n-1})$ the timestep size. The mass time derivatives in Eqs. (1c)-(1d) are discretized using discrete mass increments, while the remaining terms are approximated at the end-of-step time $t_{n}$. This leads to the following mesh-dependent, fully discrete, weak version of (1):

[TABLE]

using the discrete spaces,

[TABLE]

with $C^{0}(\overline{\Omega})$ and $L^{2}(\Omega)$ the space of continuous and square Lebesgue-integrable functions on $\overline{\Omega}$ and $\Omega$, respectively, $\mathbb{Q}_{1}(K)$ the space of trilinear polynomials in $K$, and $\mathbb{P}_{0}(K)$ the piecewise constant space in $K$. In (5c)-(5d), $F_{\ell,n}^{f}$ denotes a discrete approximation to the inter-cell mass flux for the $\ell$ fluid phase, such that $F_{\ell,n}^{f}\approx-\int_{f}\mathbf{w}_{\ell,n}\cdot\mathbf{n}_{f}\;\mbox{d}\Gamma$, with the summation taken over faces in $\mathcal{T}^{h}$ that do not belong to the (impervious) boundary $\Gamma_{f}^{N,h}$. The symbol $\llbracket\cdot\rrbracket_{f}$ denotes the jump of a quantity across a face $f$ in $\mathcal{T}^{h}$. For an internal face, $\llbracket\chi\rrbracket_{f}=({\chi}_{\left|L\right.}-{\chi}_{\left|K\right.})$, with ${\chi}_{\left|L\right.}$ and ${\chi}_{\left|K\right.}$ the restriction of $\chi$ on cells $K$ and $L$ sharing $f$, respectively. For a face belonging to the domain boundary, the jump expression reduces to $\llbracket\chi\rrbracket_{f}=-{\chi}_{\left|K\right.}$.

Discrete approximations for the displacement, pressure, and saturation field are given by

[TABLE]

with $\{\boldsymbol{\eta}_{i}\}$, $\{\psi_{j}\}$, and $\{\chi_{k}\}$ bases for $\boldsymbol{\mathcal{U}}^{h}$, $\mathcal{S}^{h}$, and $\mathcal{P}^{h}$, respectively. The nodal displacement components $\{u_{i,n}\}$, cell-centered wetting phase saturations $\{s_{j,n}\}$, and cell-centered pressures $\{p_{k,n}\}$ are collected in algebraic vectors $\mathsf{u}_{n}\in\mathbb{R}^{n_{u}}$, $\mathsf{s}_{n}\in\mathbb{R}^{n_{s}}$, and $\mathsf{p}_{n}\in\mathbb{R}^{n_{p}}$.

Remark 3.1*.*

It should be noted that the $\mathbb{Q}_{1}-\mathbb{P}_{0}-\mathbb{P}_{0}$ space adopted here is not automatically inf-sup stable. If the fluid and solid constituents are nearly incompressible, and if the permeability or timestep are very small, the resulting discrete system may allow for non-physical pressure modes, and pressure checkerboarding may be observed. To address this issue, we include additional pressure jump stabilization terms to the mass balance equations following the macroelement stabilization scheme proposed in KekSil92 . Techniques for stabilizing mixed discretizations of multiphase poromechanics are not well studied to date, however, and will be the subject of a separate contribution. Whenever sufficiently large timesteps or compressible constituents are used, spurious pressure modes are avoided, and such conditions are frequently satisfied in practice.

A linear two-point flux approximation (TPFA) scheme is adopted for the intercell numerical fluxes $F_{\ell,n}^{f}$. For an internal face $f$ shared by cells $K$ and $L$ the fluxes are computed as

[TABLE]

where the density $\rho^{\texttt{upw}}_{\ell,n}$ and mobility $\lambda^{\texttt{upw}}_{\ell,n}$ are evaluated using single-point upstream weighting (SPU) according to the sign of the potential $\Phi^{f}_{\ell,n}$, i.e. the so-called geometric part of the $\ell$-fluid phase mass flux (AziSet79, ). In (8), $\Upsilon^{f}$ is the transmissibility coefficient between cell $K$ and $L$, $z_{K}$ and $z_{L}$ are the elevation for cell $K$ and $L$ evaluated at the respective centroids, and $\varrho^{f}_{\ell}$ is the mass density averaged at face $f$ as ZhoTchMal11

[TABLE]

For a boundary face, $F_{\ell,n}^{f}$ is still computed using Eqs. (8)-(9), with the neighboring cell index $L$ being replaced by the face index. See A for details concerning the computation of the transmissibility coefficient $\Upsilon^{f}$.

Remark 3.2*.*

The linear TPFA scheme is the industry standard for reservoir simulation because of its robustness, ease of implementation, and monotonicity-preserving properties. However, the consistency of the two-point scheme is not guaranteed for arbitrary grid and permeability configurations. Therefore, the linear TPFA approach may lead to inaccurate results for highly distorted grids or full-permeability tensors. In such cases more sophisticated discretization methods like multipoint and/or nonlinear schemes are needed. For details on recent developments for heterogeneous anisotropic diffusion problems see Dro14 ; TerMalTch17 ; Sch_etal18 and references therein.

The source terms in Eqs. (5c)-(5d) are used to capture the effect of injection and production wells, modeled as line sources. These sources are governed by a suitable well model—typically derived by assuming radial flow in the immediate vicinity of the wellbore—relating well control parameters, such as bottomhole pressure, to mass flow rates through the wellbore (Pea78, ). In this work, we assume each well segment is vertical and connected to the centroid of a cell. In particular, for a cell $K$ connected to a well the source terms are computed as (CheHuaMa06, )

[TABLE]

with

[TABLE]

where $\bar{p}_{bh,n}$ is the bottom hole pressure, $\delta\left(\mathbf{x}-\mathbf{x}_{K}\right)$ is the Dirac function with $\mathbf{x}_{K}$ the cell centroid position vector, and $WI$ is the well index. In the present work all wells are assumed to operate under specified bottom-hole pressure conditions for simplicity. Clearly, the condition $\Phi^{K,W}_{\ell,n}<0$ must hold true to have fluid flow from the reservoir into a production well. Similarly, injection wells require $\Phi^{K,W}_{\ell,n}>0$. Note that injection and production wells differ in the treatment of the phase mobilities. In Eq. (10b), the sum of phase mobilities is used for injector wells so that saturation conditions are reflected in the wellbore. This is the standard approach adopted in reservoir simulators when no wellbore crossflow occurs Ecl13 . The well index depends on the well properties, local permeability, and cell dimensions. For a hexahedral $\boldsymbol{\kappa}$-orthogonal grid—i.e., when the principal axes of the permeability tensor are aligned with the grid—$WI$ is given as

[TABLE]

where $\kappa_{x}$, $\kappa_{y}$, and $\kappa_{z}$ are the diagonal (principal) components of $\boldsymbol{\kappa}$; $h_{x}$, $h_{y}$, and $h_{z}$ are the cell dimensions; $r_{w}$ is the well radius; and $s_{k}$ is the well skin factor.

4 Newton-Krylov Solver

The discrete form of the multiphase poromechanical problem (1) is obtained by introducing expressions (7), (8), and (10) into the weak form (5). This leads to a system of nonlinear equations at time $t_{n}$,

[TABLE]

for the latest solution vector $\boldsymbol{\mathsf{x}}_{n}=\{\mathsf{u}_{n},\mathsf{s}_{n},\mathsf{p}_{n}\}$. The solution $\boldsymbol{\mathsf{x}}_{n-1}$ is known from the previous timestep. Starting from an initial guess $\boldsymbol{\mathsf{x}}_{n}^{0}$, Newton’s method is used to drive the norm of the combined residual vector to below a specified relative tolerance, $\|\boldsymbol{\mathsf{r}}^{k}\|/\|\boldsymbol{\mathsf{r}}^{0}\|<\tau$ . At each iteration $k$, the Newton update is given by

[TABLE]

with Jacobian matrix $\boldsymbol{\mathsf{A}}^{k}=\partial\boldsymbol{\mathsf{r}}^{k}/\partial\boldsymbol{\mathsf{x}}_{n}$. In the current context, the linearization of three coupled governing equations leads to a Jacobian system characterized by an inherent $3\times 3$ block structure,

[TABLE]

Detailed expressions for the residual vectors and sub-matrices in $\boldsymbol{\mathsf{A}}$ are given in B.

The linear solution of the Jacobian system (15) at each Newton iteration is the most memory-demanding and time-consuming computational kernel in a multiphase poromechanics simulation. It is also the most difficult component to address in a scalable way on large parallel platforms. In this work, the system (15) is solved iteratively using a nonsymmetric Krylov subspace solver—specifically the generalized minimal residual (GMRES) method (SaaSch86, ). For large problems, practical convergence of a Krylov method is only possible provided that an effective preconditioner is available. Using right-preconditioning, we introduce a preconditioning operator $\boldsymbol{\mathsf{M}}^{-1}$ and work with the modified system

[TABLE]

In the next section, we describe an efficient and scalable two-stage preconditioner that can address both the ill-conditioning and coupled nature of the system at hand.

Remark 4.1*.*

For simplicity of presentation, we have restricted the discussion to a model system with two immiscible phases. Nevertheless, the approach devised here is naturally extendible to a larger number of phases and components. In this case, the saturation unknowns $\mathsf{s}_{n}$ should be understood to include all non-pressure, cell-based unknowns, e.g. phase compositions.

Remark 4.2*.*

Multiphase flow introduces strong nonlinearities, and a naïve implementation of Newton’s method may fail to converge when large timesteps are taken. More robust simulators include adaptive timestepping and/or globalization strategies to expand the neighborhood of convergence. Here, we apply a simple backtracking line search whenever a Newton step fails to satisfy the necessary descent conditions. We remark that smart strategies to address the intrinsic nonlinearities in these systems remain an active area of research.

Remark 4.3*.*

The system of residual equations (13) stem from the discretization of different governing equations that are not dimensionally consistent with one another. As a result, it is quite easy for the resulting discrete equations to be poorly scaled with respect to one another, causing a variety of problems at the linear and nonlinear solver level. A good choice of dimensionally consistent units can typically avoid many of these problems. As an extra safeguard, we also introduce a simple block-row scaling to ensure that the flow and mechanics equations have roughly equal magnitude before entering the solution routines.

Remark 4.4*.*

The Jacobian systems (15) is solved to a desired linear tolerance. Thus two tolerances must be selected: the Newton tolerance $\tau$ and the Krylov tolerance $\eta$. While a constant Krylov tolerance can be employed, this choice will often lead to oversolving when Newton iterates are far from the converged solution. That is, significant work will be invested in solving system (15) to a stringent tolerance, even though a much cruder solution is often acceptable for early Newton iterates. A more efficient strategy is to adapt the linear solver tolerance to reflect the observed nonlinear convergence profile—often referred to as an inexact Newton method. In this way, linear solver work is minimized in early Newton iterates. A stringent linear solver tolerance is then used in later iterates to maintain the rapid nonlinear convergence behavior of Newton’s method in a neighborhood of the solution. In particular, a good approach is the forcing strategy proposed in EisWal94 . Given parameters $\gamma\in(0,1]$, $\omega\in(1,2]$, and $\eta_{0}\in[0,1]$, set

[TABLE]

where $\eta^{k}$ is the relative Krylov tolerance adopted at Newton step $k$. In practice, an additional safeguard is also included to prevent the tolerance from becoming too small far away from the solution,

[TABLE]

The specific choices $\gamma=0.9$, $\omega=2$, and $\eta_{0}=10^{-1}$ work well in practice. As the primary focus of the current work is the effectiveness of the linear solver and proposed preconditioning strategy, however, in the numerical examples below we disable this adaptive tolerance and instead set a fixed linear solver tolerance of $\eta=10^{-6}$. This stringent tolerance leads to oversolving, but it makes the assessment of linear solver performance more straightforward.

5 Two-Stage Preconditioner

The preconditioner operator $\boldsymbol{\mathsf{M}}^{-1}$ for the Jacobian matrix $\boldsymbol{\mathsf{A}}$ is constructed using two ideas. The first component is a particular block factorization of the system matrix, which allows us to break the coupled multi-physics problem into simpler sub-problems. The second component is sparse Schur-complement approximation, allowing us to circumvent certain dense operators which would otherwise appear during the block factorization process.

5.1 Mechanics / Flow Partitioning

For clarity of presentation, it is helpful to first re-partition the Jacobian system into a $2\times 2$ form in which flow variables (saturations and pressures) are grouped and separated from the mechanics variables (displacements):

[TABLE]

A common approach for constructing preconditioning operators for multiphysics systems is to begin with a block factorization of the matrix BraPas88 ; BenGolLie05 ,

[TABLE]

The operator $\mathsf{S}_{\mathit{ff}}$ is the Schur-complement of $\mathsf{A}_{\mathit{uu}}$ in $\boldsymbol{\mathsf{A}}$. In anticipation of developments to come, we remark that this is a “first-level” Schur-complement, corresponding to elimination of one field from the original $3\times 3$ system. The second-level Schur-complement, stemming from elimination of two fields, will appear later in the discussion.

The guiding idea is to use the inverse of the lower-triangular factor as a template, and seek a preconditioning operator $\boldsymbol{\mathsf{M}}^{-1}\approx\boldsymbol{\mathsf{L}}^{-1}$, where

[TABLE]

This strategy is motivated by the observation that $\boldsymbol{\mathsf{L}}^{-1}\boldsymbol{\mathsf{A}}=\boldsymbol{\mathsf{U}}$. A property of triangular matrices is that their eigenvalues appear on the diagonal, from which we conclude that the matrix $\boldsymbol{\mathsf{U}}$ has a single eigenvalues $\lambda=1$ with multiplicity $n$. A direct application of $\boldsymbol{\mathsf{L}}^{-1}$ would therefore provide perfect conditioning. In a practical setting, however, we can only approximate $\boldsymbol{\mathsf{L}}^{-1}$ and so the construction of a useful scheme requires additional work.

The first challenge arises in the treatment of the Schur-complement operator. It may be expanded as

[TABLE]

This operator is dense and difficult to work with due to the presence of $\mathsf{A}_{\mathit{uu}}^{-1}$. As a first step, we therefore seek to replace the exact $\mathsf{S}_{\mathit{ff}}$ with a sparse approximation.

A first simplification comes from an examination of the block $\mathsf{A}_{\mathit{us}}$. This operator represents the impact of saturation changes on the momentum balance equations. Specifically, a change in saturation changes the mixture density appearing in equation (1b), which will impact the displacement field. In general, however, this coupling is quite weak, particularly if the saturation field evolves slowly in time or if phase density differences are small. For preconditioning purposes, we simply assume $\mathsf{A}_{\mathit{us}}\approx 0$ and eliminate the dense terms from the first block column of $\mathsf{S}_{\mathit{ff}}$.

For the remaining dense terms, we rely on a physically-motivated approach that has proven effective for single-phase poromechanics, the fixed-stress approximation. This approach was first introduced as an operator splitting method for sequential solution of poromechanical systems SetMou98 ; KimTchJua11a ; KimTchJua13 ; MikWhe13 . The sequential method was later reinterpreted in WhiCasTch16 as a block-preconditioned Richardson iteration involving a particular Schur-complement approximation. The authors showed that much better performance could be achieved by embedding this same preconditioning strategy within a Krylov iteration. Here, we consider the multiphase extension of this idea.

Let $\sigma_{v}=(K_{dr}\epsilon_{v}-bp)$ denote the mean total stress. Assuming the mean total stress remains fixed over an increment, i.e. $\Delta\sigma_{v}=0$, the volumetric strain can be expressed in terms of the pressure increment as

[TABLE]

Substituting expression (26) in the constitutive relationship (4) during the linearization process removes the explicit dependence of the porosity on the displacement vector, with the resulting discrete mass balance equations uncoupled from the linear momentum balance. Using this decoupling argument for computing an approximation to the exact Schur complement yields

[TABLE]

where the entry $(i,j)$ of matrices $\mathsf{D}_{sp}^{\textsc{(fs)}}$ and $\mathsf{D}_{pp}^{\textsc{(fs)}}$ are given by

[TABLE]

Noting the cell-wise constant interpolation used for both pressure and saturation, the matrices $\mathsf{D}_{sp}$ and $\mathsf{D}_{pp}$ are diagonal. Expanding the partial derivatives, the diagonal entries for an element $K_{i}\in\mathcal{T}^{h}$ with volume $V_{i}$ are

[TABLE]

The approximate Schur-complement defined in (27) is therefore sparse, with a storage pattern corresponding to the original finite-volume stencil of the mass balance equations. The fixed stress assumption only modifies the matrix diagonals associated with mass accumulation terms.

Remark 5.1*.*

For finite-volume based simulation of multiphase problems that account for a larger number of phases (or compositions), the matrix $\mathsf{D}_{\mathit{sp}}^{\textsc{(fs)}}$ will have a block-diagonal structure with blocks having dimensions $n_{s}^{K}\times 1$, where $n_{s}^{K}$ denotes the total number of non-pressure degrees-of-freedom per element.

Remark 5.2*.*

The relationship between volumetric strain and fluid pressure defined by equation (26) is an approximation based on the assumption that the loading path is primarily volumetric. The actual relationship will depend both on material properties and on the specific loading conditions for a given problem. While the modulus $K_{dr}$ will provide good convergence in all cases, it may not be the optimal value for any given problem. For example, under one-dimensional loading, a uniaxial modulus will more accurately approximate the coupling condition. See, for example, recent work in MikWhe13 ; BauRadKoc17 on using different fixed stress assumptions to improve convergence rates.

Remark 5.3*.*

An algebraic generalization of the fixed-stress approximation was proposed in Kle_etal16 based on a row-sum lumping (RSL) strategy. Using a probing vector $\mathsf{e}=\{1,\ldots,1\}^{T}\in\mathbb{R}^{n_{p}}$, the entries of the diagonal matrices defined in (28) can be computed as

[TABLE]

Here, $\texttt{diag}(\cdot)$ indicates construction of a diagonal matrix using an input vector. Computing matrices $\mathsf{D}_{\mathit{sp}}^{\textsc{(rsl)}}$ and $\mathsf{D}_{\mathit{pp}}^{\textsc{(rsl)}}$ in this manner requires two solves with $\mathsf{A}_{\mathit{uu}}$ up to a desired accuracy, and four matrix-vector multiplications. If a good preconditioner for $\mathsf{A}_{\mathit{uu}}$ is available, it can be used to approximate the action of $\mathsf{A}_{\mathit{uu}}^{-1}$ in the construction. Note that it is not necessary to update these vectors at each nonlinear iteration if the relevant Jacobian blocks do not change significantly. An alternative algebraic construction for the diagonal approximation $\mathsf{D}_{\mathit{pp}}^{\textsc{(fs)}}$ was presented and analyzed in CasWhiFer16 for single-phase poromechanics using an element-by-element approach.

Returning now to the template $\boldsymbol{\mathsf{M}}^{-1}\approx\boldsymbol{\mathsf{L}}^{-1}$, the top level preconditioning operator is defined as

[TABLE]

where two sub-preconditioners have been identified,

[TABLE]

corresponding to a mechanics preconditioner acting on the displacement components, and a flow preconditioner acting on the combined saturation and pressure components. We will describe particular forms for them in subsequent sections. In general, though, we remark that the sparse Schur complement system can be viewed as the discretization of a standard multiphase flow problem without mechanics, where the fixed stress assumption provides a particular form for the rock-compressibility term. Therefore, much of the existing knowledge on preconditioning multiphase flow problems may be brought to bear in designing $\mathsf{M}_{\mathit{ff}}^{-1}$. The matrix $\mathsf{A}_{\mathit{uu}}$ is an elasticity operator, and many existing tools also exist for constructing $\mathsf{M}_{\mathit{uu}}^{-1}$. Therefore, the block factorization strategy provides a convenient and modular framework for assembling a coupled multiphase poromechanics solver from highly-tuned component pieces.

5.2 Mechanics Preconditioner

The elasticity block $\mathsf{A}_{\mathit{uu}}$ is an elliptic operator that is well suited to multigrid techniques Baker:2010 . Here, we take a simple approach known as the separate displacement component approximation, introduced in AxeGus78 . If the displacement degrees of freedom are ordered by coordinate direction, then the resulting matrix has its own $3\times 3$ structure,

[TABLE]

Rather than building a preconditioner using the complete matrix, we consider the sparser approximation

[TABLE]

in which the off-diagonal blocks have been eliminated. Using Korn’s inequality, one can show that $\mathsf{A}^{\text{(SDC)}}_{\mathit{uu}}$ is spectrally equivalent to $\mathsf{A}_{\mathit{uu}}$ Bla94 ; GusLin98 . Note that this approximation breaks down in the incompressible elasticity limit, Poisson ratio $\nu\to 0.5$, though this is not a concern for most geologic materials. We then construct an algebraic multigrid preconditioner from the block-diagonal approximation,

[TABLE]

which can be thought of as three scalar AMG preconditioners, one for each coordinate direction. More specific details regarding the AMG implementation are summarized with the numerical examples below.

The mechanics preconditioner is the most expensive operator to construct during the solver setup phase. For an elastic model, however, the matrix $\mathsf{A}_{\mathit{uu}}$ is constant during the entire simulation. This AMG preconditioner therefore only needs to be initialized once, and its setup cost is quickly amortized over many solver calls. This is not the case for the flow preconditioner, which is typically re-initialized within every Newton iteration.

5.3 Flow Preconditioner

The preconditioning operator for the approximate Schur complement $\mathsf{S}^{\textsc{(fs)}}_{\mathit{ff}}$ is built following the Constrained Pressure Residual (CPR) approach Wal83 ; WalKenLit85 ; Cao_etal05 . This is a multistage strategy that combines preconditioners with different features in a multiplicative fashion. In general, given a matrix $\mathsf{B}$ in $\mathbb{R}^{n\times n}$, a multiplicative multistage preconditioner $\mathsf{M}_{B}^{-1}$ for $\mathsf{B}$ can be written as Cao_etal05

[TABLE]

with $n_{\text{st}}$ the number of stages, $\mathsf{M}_{i}^{-1}$ the $i$th preconditioner for $\mathsf{B}$, and $\mathsf{I}$ the identity matrix in $\mathbb{R}^{n\times n}$.

The design of the CPR multistage approach—the industry standard for fully-implicit simulation of multiphase flow—is guided by the underlying nature of the governing equations. Fields that exhibit near-elliptic behavior, such as pressure, are characterized by long-range (global) coupling. They require multilevel preconditioners to reduce error modes that are smooth across the computational domain. Conversely, fields with near-hyperbolic character, such as saturation, display short-range (local) coupling. Cheaper strategies can be employed to remove the associated error modes. Standard CPR is implemented as a two-stage method that combines a global (${\mathsf{M}}^{-1}_{G}$) preconditioner with a local preconditioner (${\mathsf{M}}^{-1}_{L}$). The overall preconditionioning operator $\mathsf{M}_{\mathit{ff}}^{-1}$ is formally

[TABLE]

In practice, $\mathsf{M}_{\mathit{ff}}^{-1}$ is applied through a sequence of vector-multiply operations.

The derivation of ${\mathsf{M}}_{G}^{-1}$ also relies on a block LU-factorization, this time of the approximate Schur-complement $\mathsf{S}_{\mathit{ff}}^{\textsc{(fs)}}$,

[TABLE]

Note that $\mathsf{S}_{\mathit{pp}}$ is the “second-level” Schur complement, obtained after recursive elimination of both displacements and saturations from the original $3\times 3$ system $\boldsymbol{\mathsf{A}}$. We again seek to sparsify the Schur-complement, this time by introducing diagonal approximations for $\mathsf{A}_{ss}$ and $\mathsf{A}_{ps}$. Two common strategies are True-IMPES reduction and Quasi-IMPES reduction.

True-IMPES reduction seeks to mimic the IMplicit-Pressure Explicit-Saturation (IMPES) time discretization scheme AziSet79 on the linear level. A column-sum lumping operation is used to lump off-diagonal entries—containing derivatives of the flux terms with respect to the saturation degrees of freedom—into the diagonal. This can be computed using the matrix transpose as

[TABLE]

with $\mathsf{e}=\{1,\ldots,1\}\in\mathbb{R}^{n_{s}}$. A simpler alternative is the Quasi-IMPES reduction, which extracts the diagonal entries from the respective matrices, without lumping in off-diagonal terms,

[TABLE]

To clarify notation, $\texttt{diag}(\cdot)$ indicates construction of a diagonal matrix using an input vector, while $\texttt{diagm}(\cdot)$ indicates construction of a diagonal matrix using an input matrix. In practice, the diagonal matrices are simply stored as a vector of their diagonal entries. Note that in a more general formulation where multiple saturation or compositional variables are present, the CPR reduction leads to block-diagonal matrices with small, dense blocks corresponding to degrees of freedom coupled in each cell. See Cao_etal05 for details.

With these diagonal approximations, the second-level Schur complement may be approximated as

[TABLE]

Note that this strategy again preserves the matrix sparsity pattern stemming from the original finite-volume stencil. The result is a near-elliptic pressure system that is well suited to AMG preconditioning. As before, we then use the inverse of the lower triangular factor appearing in (38) as a preconditioner template. This leads to the concrete definition of the global preconditioner

[TABLE]

In equation (42), we note that the matrix $A_{\mathit{ps}}$ could also be replaced here with its diagonal approximation. We have generally found better performance using the complete matrix, however, despite the larger cost for a matrix-vector multiply operation.

The global preconditioner $\mathsf{M}_{G}^{-1}$ corrects low-frequency errors associated with the pressure unknowns, but it only performs a simple (Jacobi) correction to the saturation unknowns using $\mathsf{D}_{ss}^{\textsc{(cpr)}}$. The second-stage preconditioner $\mathsf{M}_{L}^{-1}$ complements this first stage by smoothing out the remaining error modes. These residual errors are primarily associated with steep saturation gradients, and are highly localized. As a result, single-level and processor-local preconditioning strategies can be quite effective.

To construct the second stage, the matrix $\mathsf{S}_{\mathit{ff}}^{\textsc{(fs)}}$ is viewed as a monolithic system. Let $\mathsf{P}$ be a permutation matrix that reorders the partitioned vector of saturation and pressure unknowns into an interleaved ordering—that is, degrees of freedom are re-ordered as $\{\{s_{1},p_{1}\},\{s_{2},p_{2}\},...,\{s_{n},p_{n}\}\}$. This permutation creates a sparsity pattern in which dense $2\times 2$ blocks appear for each grid cell, with these cell blocks connected by the TPFA finite volume stencil. This hierarchical structure also appears in the more general case when multiple saturation or compositional variables exist. The size of the dense blocks corresponds to the total number of unknowns per cell. Block versions of relaxation or incomplete factorization preconditioners are therefore appealing, as dense multiplication and inversion operations can be applied to the small blocks.

Common strategies for $\mathsf{M}_{L}^{-1}$ are relaxation preconditioners (e.g. Jacobi, Gauss-Seidel) or incomplete factorizations (e.g. ILU(k), ILUT). In this work, we explore two options. Option 1 uses several sweeps of hybrid block Gauss-Seidel (hbgs), i.e.

[TABLE]

where the input matrix is the interleaved version of the partitioned flow system. Option 2 uses one sweep of processor-local, pointwise ILU(k) (iluk) on the interleaved system.

5.4 Final Algorithm

Using the developments above, the concrete implementation of the preconditioning strategy can be summarized as follows. In the simulator, linear solver calls are broken into two phases: a setup phase and an iteration phase. During the setup phase, auxiliary vectors and matrices associated with the fixed stress and CPR reductions are assembled. Then sub-preconditioners $\mathsf{M}_{\mathit{uu}}^{-1}$, $\mathsf{M}_{\mathit{pp}}^{-1}$, and $\mathsf{M}_{\mathit{L}}^{-1}$ can be constructed. The elasticity preconditioner $\mathsf{M}_{\mathit{uu}}^{-1}$ is the most expensive to build, but it only needs to be initialized once during the entire simulation. Once the necessary operators are assembled, a non-symmetric Krylov solver (e.g. GMRES) may be called with a preconditioning routine implementing the strategy proposed here. As an aid to the reader, Algorithm 1 summarizes the concrete steps necessary to apply the preconditioner $\boldsymbol{\mathsf{M}}^{-1}$ to an input vector. We observe that steps 1-5 implement a block-lower-triangular preconditioner involving two AMG sub-preconditioners and one Jacobi (diagonal) preconditioner. Steps 6-7 then implement a second stage correction to the saturation and pressure fields using the cheap local smoother.

Remark 5.4*.*

While this work focuses on “monolithic” solver strategy using a preconditioned Krylov iteration, we remark that a “sequential” solver strategy can be obtained easily by simply embedding the proposed preconditioner within a Richardson iteration. Step 1 of the preconditioner implements a correction to the displacement field, while steps 2-8 implement a subsequent update to the flow variables. One would then iterate back and forth between mechanics and flow in a sequential fashion within the fixed point iteration. See WhiCasTch16 for further details. In general, however, the Krylov approach will provide superior convergence properties.

6 Numerical Examples

We present two numerical examples to test the performance of the proposed preconditioner: (1) a weak scaling study using a simple, synthetic configuration; and (2) a strong scaling study using a more realistic, highly-heterogeneous reservoir. The code used for this study, Geocentric, relies heavily on the deal.ii Finite Element Library (BanHarKan07, ) for discretization functionality. Algebraic multigrid, relaxation, and incomplete factorization sub-preconditioners are provided by Hypre FalYan02 . All numerical results were obtained on a cluster of Intel Xeon E5-2695 processors. Each computational node contains two 18-core processors sharing 128 GiB of memory, with Intel Omni-Path interconnects between nodes. We use pure MPI-based parallelism.

In general, the performance of the proposed algorithm is heavily dependent on the specific design of the sub-preconditioners. The numerical results below are based on the following choices: For the separate-displacement-component elasticity operator, we use an AMG preconditioner with one level of aggressive coarsening, one sweep of hybrid forward $\ell_{1}$-Gauss-Seidel Baker:2011 for the down cycle and one sweep of hybrid backward $\ell_{1}$-Gauss-Seidel for the up cycle. The coarsest grid is solved directly with Gaussian elimination. For the first-stage flow system, we use quasi-IMPES reduction. The pressure system uses a standard scalar AMG approach with a Hybrid Modified Independent Set (HMIS) coarsening strategy DeSterck06 . Similar to the elasticity block, we use one sweep of hybrid forward $\ell_{1}$-Gauss-Seidel for the down cycle, one sweep of hybrid backward $\ell_{1}$-Gauss-Seidel for the up cycle, and a direct solve with Gaussian elimination for the coarsest grid. For the second stage correction, we explore two options. For the first benchmark in §6.1 we use three sweeps of hybrid block Gauss-Seidel. For the second benchmark in §6.2 we use one sweep of processor-local, pointwise ILU(0).

Details of the material and parameter values used in both examples are summarized in Table 1. In both examples, fluid phases follow the density model

[TABLE]

where $c_{l}$ is the fluid compressibility and $\rho_{l}^{o}$ is the phase density at a reference pressure $p^{o}$. This reference pressure is taken as initial formation pressure. For the rock, solid grains are intrinsically incompressible with Biot coefficient $b=1$, maximizing the coupling strength between the mechanical and flow problems. The relative permeability relationships are

[TABLE]

where $s_{wr}$ and $s_{nwr}$ are the wetting and non-wetting residual saturations. At $t=0$, the formations are assumed to be in hydrostatic and lithostatic equilibrium, with a uniform saturation $s(\mathbf{x},0)=s_{wr}$. In both examples, target well pressures are not applied instantaneously, but rather ramped up linearly over a fixed time period (1 day) to relax the stiffness of the nonlinear problem. We use a telescoping time-stepping scheme, where simulation timesteps grow exponentially from an initial value to a maximum value, after which they remain constant for the remainder of the simulation.

6.1 Staircase Benchmark

For the first benchmark, we consider the synthetic configuration illustrated in Figure 1. Injection and production wells are placed at opposite ends of a high-permeability channel that winds its way in a “staircase” spiral through a low permeability host rock. The reservoir channel has a 1000:1 permeability contrast with surrounding seal units, with isotropic permeability in both regions. All external boundaries are assumed to be no flow. For mechanics, the side and bottom boundaries have a zero normal-displacement constraint (roller boundary), while the top surface is allowed to deform freely. The water injector and fluid production wells completely perforate the channel. The injection (respectively, production) well reaches a target bottomhole pressure of 5 MPa above (below) initial formation pressure after 1 day of buildup (drawdown). An illustrative snapshot of the simulation results is presented in Figure 1. While the geometry is artificial, the problem is designed to exhibit strong coupling and provides a straightforward description for comparison with other codes.

The results of a weak-scaling study using this configuration are presented in Table 2. We begin with an 88k degree-of-freedom mesh run on one core. The mesh is then uniformly refined and run on eight cores, preserving (approximately) the same number of unknowns per core. This process is continued out to 512 cores and 42M degrees of freedom. We then repeat the process with 88k degrees-of-freedom on 2 cores and scale out to 1024 cores. One can also get an indication of strong scaling performance by comparing results in Tables 2a and 2b for the same problem size but different core counts.

Because the problem is nonlinear, we observe an increase in the average number of Newton iterations per timestep as the mesh is refined. As desired, however, the number of GMRES iterations per linear solve is relatively insensitive to the problem size, showing only minor growth with refinement. Using more complex smoothers, we can further drive down the iteration count but at the cost of overall timing performance. Table 2 also presents timing results. Both the setup and solver iteration phases exhibit some growth, though even at large processor counts the solve time is quite satisfactory. We note that the actual work per processor in this problem is quite small, and the timing growth is a symptom of communication costs in the AMG setup and solver iteration phases. In general, however, the overall framework provides a robust basis for scalable performance.

6.2 SPE10-based Benchmark

As a more realistic test case, we consider a water flood in a heterogeneous reservoir. The setup is based on Model 2 of the 10th SPE Comparative Solution Project ChrBlu01 but equipped with poroelastic mechanical behavior. In order to provide realistic mechanical boundary conditions, we add 300-m thick overburden and underburden layers to the original SPE10 reservoir. A sketch of the simulated domain is shown in Figure 2. The computational grid consists of $60\times 220\times 252$ cells. The total number of degrees of freedom is $n_{\text{total}}$ = 16,730,559 ($n_{u}$ = 10,077,759; $n_{s}$ = 3,326,400; and $n_{p}$ = 3,326,400). The original SPE10 anisotropic permeability distribution and relative permeability curves are used—details provided in ChrBlu01 —while the porosity field has been slightly modified so that the lowest porosity for a cell is thresholded at 1% to remove extremely low porosity values present in the original data files. An isotropic permeability value equal to $0.01$ mD is assigned to the overburden and underburden layers, and porosity is set to 1%. Homogeneous Young’s modulus $E$ = 5000 MPa, Poisson’s ratio $\nu$ = 0.25, and Biot’s coefficient $b$ = 1.0 are assumed everywhere. All boundaries are impervious to fluid flow and constrained to have zero normal displacement, except for the top which is allowed to deform freely. Five wells penetrate the entire reservoir as shown in Figure 2. The injector well, located in the middle of the domain, and the production wells, located at the corners, operate at a constant +5 and -5 MPa overpressure, respectively. Note that these well conditions differ from the original SPE10 specification, which include rate-controlled wells. Overall, 100 days of injection are simulated.

Table 3 reports results of a strong-scaling study in which the total problem size is kept fixed but allocated across an increasing number of cores. The GMRES iterations to convergence are well-controlled and mostly insensitive to the parallel partitioning. Good overall timing efficiency is observed with only a mild decrease from ideal performance for the 576 core case. The results suggest that one can use the proposed solver strategy and large processor counts to efficiently drive down long simulation times even for challenging heterogeneous problems.

7 Conclusion

In this work, we have presented a preconditioning scheme suitable for fully-implicit simulation of multiphase flow and transport with geomechanics. We have used physics-based arguments to break the inherently coupled problem into modular sub-problems for which “black-box” algebraic methods can be effective. This modularity also provides significant flexibility to test other approaches. For example, while we have used a CPR-like strategy for the reduced multiphase flow problem, research into effective reservoir simulation preconditioners remains an active area of research. If an alternative flow preconditioner has been identified, the partitioning framework here provides a simple strategy for adding poromechanical physics.

Regarding future work, a number of avenues would be interesting to explore. First, it is well known that the coupling strength between different PDEs evolves with both time and timestep size. Strong coupling is often observed at early time and small timesteps, while weaker coupling is observed at late time and large timesteps. One might consider an adaptive solver strategy where the coupling strength between various quantities is detected based on the magnitude of coupling contributions in the preconditioner. When sufficient “uncoupling” is detected between one or more fields via these coupling indicators, one might switch to a cheaper solver strategy.

Also, this work has focused almost exclusively on the solver strategy at the linearized equation level. There are likely significant performance improvements to be gained by embedding similar ideas at the nonlinear level. Finally, while the physics-based splitting strategies used here are effective, they are PDE-specific and cannot be immediately extended to other coupled systems. For example, the addition of thermal or fractured reservoir behavior would change the underlying PDEs and render the current solver strategy incomplete. There is therefore significant room for both new physics-based solvers as well as algebraic-coupling methods that can be (to the extent possible) agnostic as to the underlying physics.

Acknowledgements

Funding for JAW, NC, QMB and DO-K was provide by LLNL through Laboratory Directed Research and Development (LDRD) Project 18-ERD-027. Funding for SK and HAT was provided by Total S.A. through the FC-MAELSTROM Project. Portions of this work were performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07-NA27344.

Appendix A TPFA transmissibility computation

The transmissibility coefficient between two connected cells $K$ and $L$ in $\mathcal{T}^{h}$ sharing a face $f$ is defined by the following harmonic average

[TABLE]

where quantities $\Upsilon^{f,K}$ and $\Upsilon^{L,f}$, also known as half transmissibility coefficients, are given by TerVas13

[TABLE]

with $|f|$ the face area, and $\mathbf{x}_{f}$ a collocation point on $f$ that allows for enforcing point-wise pressure continuity across the face in case of incompressible single-phase flow. In our implementation $\mathbf{x}_{f}$ is selected as the intersection of the face and the line connecting the cell centroids $\mathbf{x}_{K}$ and $\mathbf{x}_{L}$ KarDur16 as shown in Fig. 1. For a boundary face $f\in\mathcal{F}_{K}\cap\Gamma_{f}^{D,h}$, transmissibility and half transmissibility coefficient coincide, i.e. $\Upsilon^{f}=\Upsilon^{f,K}$, and $\mathbf{x}_{f}$ is chosen as the orthogonal projection of $\mathbf{x}_{K}$ on $f$. For a non-planar face, a plane characterized by a mean normal $\mathbf{n}_{f}$ and a collocation point $\mathbf{x}_{f}$ are associated to $f$.

Appendix B Finite element and finite volume vectors and matrices

Repeated indices in the expressions given in this appendix do not imply summation.

B.1 Discrete residual vectors expressions

The residual vector $\mathsf{r}_{u}$ (13) is assembled as sum of element contributions. The residual vector $\mathsf{r}_{\ell}$ (13), $\ell=\{w,nw\}$, is assembled as sum of element ($\mathsf{r}^{K}_{\ell}$) and interface ($\mathsf{r}^{f}_{\ell}$) contributions, i.e. $\mathsf{r}_{\ell}=\mathsf{r}^{K}_{\ell}+\mathsf{r}^{f}_{\ell}$. Their expressions read:

[TABLE]

B.2 Jacobian sub-matrices expressions

The sub-matrices appearing in the global Jacobian matrix (15) are constructed as sum of element ($\mathsf{A}^{K}_{lm}$) and face ($\mathsf{A}^{f}_{lm}$) contributions, namely

[TABLE]

The expressions for sub-matrices $\mathsf{A}_{lm}^{K}$ read:

[TABLE]

The partial derivatives in the above equations are expanded as:

[TABLE]

The finite volume scheme provides a direct functional dependence for intercell fluxes $F^{f}_{\ell,n}$ on the algebraic vectors $\mathsf{s_{n}}$ and $\mathsf{p_{n}}$. The global expressions for sub-matrices $\mathsf{A}_{lk}^{f}$ are therefore

[TABLE]

The partial derivatives in the above equations are expanded as:

[TABLE]