A block preconditioner for non-isothermal flow in porous media

TL;DR

This paper introduces a novel block preconditioner tailored for non-isothermal flow in porous media, effectively incorporating thermal effects like heat diffusion into reservoir simulations.

Contribution

It develops a new block preconditioner with an efficient Schur complement approximation specifically designed for non-isothermal flow models.

Findings

The preconditioner is scalable with problem size and parallelization.

It effectively captures thermal effects such as heat diffusion.

The solver improves convergence in thermal reservoir simulations.

Abstract

In petroleum reservoir simulation, the industry standard preconditioner, the constrained pressure residual method (CPR), is a two-stage process which involves solving a restricted pressure system with Algebraic Multigrid (AMG). Initially designed for isothermal models, this approach is often used in the thermal case. However, it does not have a specific treatment of the additional energy conservation equation and temperature variable. We seek to develop preconditioners which better capture thermal effects such as heat diffusion. In order to study the effects of both pressure and temperature on fluid and heat flow, we consider a model of non-isothermal single phase flow through porous media. For this model, we develop a block preconditioner with an efficient Schur complement approximation. Both the pressure block and the approximate Schur complement are approximately inverted using an AMG V-cycle. The resulting solver is scalable with respect to problem size and parallelization.

Figures

Tables

Table 1. Table 1 : Parameters for the Bennison viscosity correlation

$A_1$	$A_2$	$A_3$	$A_4$
-0.8021	23.8765	0.31458	-9.21592

Table 2. Table 2 : Physical parameters for test cases

Initial pressure	4.1369$\times {10}^5$ Pa
Conductivity of oil	0.15 Wm-1K-1
Conductivity or rock	1.7295772056 Wm-1K-1
Specific heat of oil	2093.4 JK-1kg-1
Specific heat of rock	920 JK-1kg-1

Table 3. Table 3 : SPE10 test cases. Average linear iterations per nonlinear iteration.

method/case	W	H	W+H	h.p. W	h.p. W+H
Block	5.88	5.42	6.60	14.5	14.0
CPR	6.67	6.27	6.69	11.4	11.0

Table 4. Table 4 : Condition numbers (upper bounds) of the different matrices and Schur complement approximations for various cases

matrix/case	H	W	W+H	h.p. W	h.p. W+H
${\tilde{S}}_{\mathrm{diag}}^{-1}S$	1.061	20.75	3.323	8.703e7	2.191e7
${\tilde{S}}_{A_{\mathrm{TT}}}^{-1}S$	1.063	28.08	4.277	4117	2467
$\tilde{S}_T^{}{}_{}{}^{-1}S$	1.023	1.097	1.1717	5.969	5.939
$A_{TT}$	5.64e5	27.88	5.64e5	2324	2.143e5
$S$	5.479e5	2.862	5.717e5	20.60	4.976e5

Table 5. Table 5 : Case I: Heater case in isotropic medium. Average linear iterations per nonlinear iteration.

method/$N$	20	40	80	160	320
Block	2.57	3.23	2.86	3.44	3.71
CPR	3.4	5.38	9.09	16.3	30.7

Table 6. Table 6 : Case II: Well case in isotropic medium. Average linear iterations per nonlinear iteration.

method/$N$	20	40	80	160	320
Block	2.43	2.43	2.86	3.28	3.71
CPR	3.71	5.71	9.86	19.4	37.4
Block (${\tilde{S}}_{A_{TT}}$)	4.57	5	5.57	6.29	6.57
Block (${\tilde{S}}_{\mathrm{diag}}$)	4.14	4.43	5.29	5.86	6.43

Table 7. Table 7 : Case III: Higher injection well case in isotropic medium. Average linear iterations per nonlinear iteration.

method/$N$	20	40	80	160	320
Block	3.67	4.38	4.7	5.10	5.52
CPR	4.71	7.31	13.1	24.7	50.6

Table 8. Table 8 : Case IV: Heater case in anisotropic medium. Average linear iterations per nonlinear iteration.

method/$N$	20	40	80	160	320
Block	2.31	2.67	3.25	3.67	3.86
CPR	3.11	4.56	8.56	15.8	30.4

Table 9. Table 9 : Case V: Well case in anisotropic medium. Average linear iterations per nonlinear iteration.

method/$N$	20	40	80	160	320
Block	2.38	3.27	4.52	4.68	5.36
CPR	2.86	3.6	4.76	7.0	12.04
Block (${\tilde{S}}_{A_{TT}}$)	9	17.1	24.1	27.9	31.6
Block (${\tilde{S}}_{\mathrm{diag}}$)	$>200$	$>200$	$>200$	$>200$	$>200$

Table 10. Table 10 : Weak scaling: 3D Heating case. Average linear iterations per nonlinear iteration.

method/num. proc.	1	2	4	8	16
Block	7.5	7.9	8.25	8.75	9.29
CPR	15.75	22.3	29.9	38.5	45.6

Table 11. Table 11 : Weak scaling: 3D Well case. Average linear iterations per nonlinear iteration.

method/num. proc.	1	2	4	8	16
Block	4.29	4.43	4.71	5.25	5.89
CPR	6.57	10.0	12.3	16.3	18.4

Table 12. Table 12 : Strong scaling: 3D Heating case. Average linear iterations per nonlinear iteration.

method/num. proc.	1	2	4	8	16
Block	8.75	8.57	8.57	9.43	9.29
CPR	38.0	43.3	44.0	44.7	45.6

Table 13. Table 13 : Strong scaling: 3D Well case. Average linear iterations per nonlinear iteration.

method/num. proc.	1	2	4	8	16
Block	6.22	5.44	5.78	6.11	5.89
CPR	14.0	16.1	16.7	17.9	18.4

Equations

$$\varphi \frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=f\text{in }{\mathbb{R}}_{+}\times \Omega ,$$

$$\mathbf{u}=-\frac{\mathbf{K}}{\mu }(\nabla p-\rho \mathbf{g}),$$

$$\varphi \frac{\partial \rho }{\partial t}-\nabla \cdot \left(\rho \frac{\mathbf{K}}{\mu }(\nabla p-\rho \mathbf{g})\right)=f\text{in }{\mathbb{R}}_{+}\times \Omega .$$

$$-\frac{\mathbf{K}}{\mu }(\nabla p-\rho \mathbf{g})\cdot \mathbf{n}=g_N\text{ on }{\Gamma }_N,\text{and}p=g_D\text{ on }{\Gamma }_D,$$

$$\varphi \frac{\partial }{\partial t}(\rho c_{v}T)+(1-\varphi )\frac{\partial }{\partial t}({\rho }_r{c}_{r}T)+\nabla \cdot (\rho c_{v}T\mathbf{u})+\nabla \cdot \mathbf{q}=f_T\text{in }{\mathbb{R}}_{+}\times \Omega ,$$

$$\mathbf{q}=-k_T\nabla T,$$

$$k_T=\varphi k_{T,r}+(1-\varphi )k_{T,f},$$

$$\varphi \frac{\partial }{\partial t}(\rho c_{v}T)+(1-\varphi )\frac{\partial }{\partial t}({\rho }_r{c}_{r}T)+\nabla \cdot (\rho c_{v}T\mathbf{u})-\nabla \cdot (k_T\nabla T)=f_T\text{in }{\mathbb{R}}_{+}\times \Omega .$$

$$\begin{gathered} \varphi \frac{\partial }{\partial t}(\rho c_{v}T)+(1-\varphi )\frac{\partial }{\partial t}({\rho }_r{c}_{r}T)-\nabla \cdot \left(\rho c_{v}T\frac{\mathbf{K}}{\mu }(\nabla p-\rho \mathbf{g})\right)-\nabla \cdot (k_T\nabla T)=f_T \\ \text{in }{\mathbb{R}}_{+}\times \Omega . \end{gathered}$$

$$-\left(\rho c_{v}T\frac{\mathbf{K}}{\mu }(\nabla p-\rho \mathbf{g})+k_T\nabla T\right)\cdot \mathbf{n}=g_N^T\text{ on }{\Gamma }_N^T,\text{and}T=g_D^T\text{ on }{\Gamma }_D^T,$$

$$\varphi \frac{\partial \rho }{\partial t}-\nabla \cdot \left(\rho \frac{\mathbf{K}}{\mu }(\nabla p-\rho \mathbf{g})\right)=f\text{in }{\mathbb{R}}_{+}\times \Omega ,$$

$$\begin{gathered} \varphi \frac{\partial }{\partial t}(\rho c_{v}T)+(1-\varphi )\frac{\partial }{\partial t}({\rho }_r{c}_{r}T)-\nabla \cdot \left(\rho c_{v}T\frac{\mathbf{K}}{\mu }(\nabla p-\rho \mathbf{g})\right)-\nabla \cdot (k_T\nabla T)=f_T \\ \text{in }{\mathbb{R}}_{+}\times \Omega , \end{gathered}$$

$$-\frac{\mathbf{K}}{\mu }(\nabla p-\rho \mathbf{g})\cdot \mathbf{n}=g_N\text{ on }{\Gamma }_N,\text{and}p=g_D\text{ on }{\Gamma }_D,$$

$$-\left(\rho c_{v}T\frac{\mathbf{K}}{\mu }(\nabla p-\rho \mathbf{g})+k_T\nabla T\right)\cdot \mathbf{n}=g_N^T\text{ on }{\Gamma }_N^T,\text{and}T=g_D^T\text{ on }{\Gamma }_D^T,$$

$$\mu (T_{\mathrm{F}})=10^{A_1{\gamma }_{\mathrm{API}}+A_2}{T}_{\mathrm{F}}^{A_3{\gamma }_{\mathrm{API}}+A_4},$$

$${\gamma }_{\mathrm{API}}=\frac{141.5}{{\gamma }_{\mathrm{SG}}}-\mathrm{131.5.}$$

$$\rho (p,T)={\rho }_0{e}^{c(p-p_0)}{e}^{\beta (T-T_0)},$$

$$f(\mathbf{x})=$$

$$-\sum _j{q}_{\mathrm{prod}}^j(p,T)\delta (\mathbf{x}-{\mathbf{x}}_{\mathrm{prod}}^j)\rho (p,T),$$

$$f_T(\mathbf{x})=$$

$$-\sum _j{q}_{\mathrm{prod}}^j(p,T)\delta (\mathbf{x}-{\mathbf{x}}_{\mathrm{prod}}^j)\rho (p,T)c_{v}T,$$

$$q=\frac{2\pi h{K}_e}{\mu \ln (r_e/r_w)}(p_{bh}-p),$$

$$r_e=\frac{0.14{\left((K_y/K_x{)}^{1/2}{D}_x^2+(K_x/K_y{)}^{1/2}{D}_y^2)\right)}^{1/2}}{0.5\left((K_y/K_x{)}^{1/4}+(K_x/K_y{)}^{1/4}\right)},$$

$$f_T=\sum _i{U}_{\mathrm{heater}}(p,T)\delta (\mathbf{x}-{\mathbf{x}}_{\mathrm{heater}}^i)(T_{\mathrm{heater}}-T),$$

$$\frac{\partial u}{\partial t}-{\nabla }^{2}u=0\text{ in }\Omega ,$$

$$u=f\text{ on }{\Gamma }_D,\nabla u\cdot \mathbf{n}=g\text{ on }{\Gamma }_N.$$

$${\int }_{E_i}\frac{\partial u}{\partial t}v\mathrm{d}x+{\int }_{E_i}\nabla u\cdot \nabla v\mathrm{d}x-{\int }_{\partial E_i}v\nabla u\cdot \mathbf{n}\mathrm{d}s=0\text{for all test functions }v,$$

$${\int }_{\partial E_i}v\nabla u\cdot \mathbf{n}\mathrm{d}s\approx \sum _{j\in \mathcal{N}(i)}{\int }_{e_{ij}}v{\mid }_{E_i}^{e_{ij}}\frac{u{\mid }_{E_j}^{e_{ij}}-u{\mid }_{E_i}^{e_{ij}}}{\Vert h_j-h_i\Vert }\mathrm{d}s.$$

$${\int }_{E_i}\frac{\partial u_i}{\partial t}{v}_i\mathrm{d}x-\sum _{j\in \mathcal{N}(i)}{\int }_{e_{ij}}{v}_i\frac{u_j-u_i}{\Vert h_j-h_i\Vert }\mathrm{d}s=0.$$

$$\frac{\partial u_i}{\partial t}|E_i|-\sum _{j\in \mathcal{N}(i)}\frac{u_j-u_i}{\Vert h_j-h_i\Vert }|e_{ij}|=0,$$

$${\int }_{\partial E_i}v\nabla u\cdot \mathbf{n}\mathrm{d}s=\sum _{j\in \mathcal{N}(i)}{\int }_{e_{ij}}v\nabla u\cdot {\mathbf{n}}_e+{\int }_{\partial E_i\cap {\Gamma }_D}v\nabla u\cdot \mathbf{n}\mathrm{d}s+{\int }_{\partial E_i\cap {\Gamma }_N}v\nabla u\cdot \mathbf{n}\mathrm{d}s,$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Models of fluid flow in porous media are used in the simulation of applications such as petroleum reservoirs, carbon storage, hydrogeology, and geothermal energy. In some cases, fluid flow must be coupled with heat flow in order to capture thermal effects. Petroleum reservoir simulation is used in optimizing oil recovery processes, which often involve heating and steam injection inside the reservoir in order to reduce the viscosity of the oil. This is especially important in the case of heavier hydrocarbons.

In the case of isothermal multiphase flow, a global pressure couples local concentration/saturation variables. The equations in the system are elliptic with respect to the pressure and hyperbolic with respect to the non-pressure variables. The industry standard constrained pressure residual (CPR) preconditioner introduced by Wallis  [1, 2] in the early 80s defines a discrete decoupling operator essentially splitting pressure and non-pressure variables, in order that each can be preconditioned separately. Indeed, the global nature of the pressure variable requires a more precise “global” preconditioning than the other variables for which “local” preconditioning is sufficient. In brief, the CPR preconditioner is a two-stage process in which pressure is solved first approximately, followed by solving approximately the full system.

A major improvement to CPR was introduced in [3] with the use of Algebraic Multigrid (AMG) [4] as a preconditioner for the pressure equation in the first stage. AMG is used as a solver for elliptic problems, usually as a preconditioner for a Krylov subspace method. Therefore, the elliptic-like nature of the pressure equation makes it an ideal candidate for the use of AMG. This improved preconditioner, often denoted CPR-AMG, is widely used in modern reservoir simulators.

The non-isothermal case adds a conservation of energy equation and a temperature (or enthalpy) variable to the system of PDEs. In the standard preconditioning approach, the energy conservation equation and the temperature unknowns are treated similarly to the secondary equations and unknowns. This means that the thermal effects are only treated in the second stage of CPR, usually an Incomplete LU factorization (ILU) method. More dense incomplete factors are often needed in the thermal case. While this results in a lower iteration count, it is not ideal in terms of computational time, memory requirements, and parallelization.

Alternatives to the usual approach were recently proposed. The Fraunhofer Institute for Algorithms and Scientific Computing (SCAI) focuses on AMG for systems of PDEs based on [5], often called System AMG (SAMG). While it seeks to replace CPR-AMG, the SAMG approach proposed in [6, 7] is quite similar in the isothermal case. Indeed, AMG is applied to the whole system, but all non-pressure variables remain on the fine level. In the thermal case, however, SAMG allows the consideration of both pressure and temperature for the hierarchy. A proper comparison with CPR-AMG for thermal simulation cases has yet to be done. Other AMG methods for systems of PDEs include BoomerAMG [8] and multigrid reduction (MGR) in the hypre library [9], as well as Smoothed Aggregation in the ML package [10]. In particular, BoomerAMG has been shown to be effective in diffusion-dominated two-phase flow problems [11], and MGR has also had some success with multiphase flow problems [12, 13].

The inclusion of temperature in the AMG hierarchy is still not well understood. The temperature is not always descriptive of the flow everywhere in the reservoir (for example in regions of faster flows). This justifies an adaptive method where only variables which are descriptive be included in the first stage of CPR. Retaining the CPR structure, Enhanced CPR (ECPR) constructs a “strong” subsystem for the first stage of CPR by looking at the coupling in the system matrix [14]. The resulting subsystem has no real physical interpretation and it is unclear if it is possible to solve it via AMG.

In the context of this paper, we consider single phase non-isothermal flow. This single phase case is relevant for geothermal models and simple reservoir simulation examples, but can also be applied to miscible displacement problems (where a concentration plays a similar role to temperature) [15]. We present a block preconditioner for the solution of the resulting coupled pressure-temperature system.

In Section 2, we present the mathematical model for non-isothermal flow in porous media and the discretization. In Section 3, we describe the preconditioning approaches for the linearized system. Numerical results for the preconditioners are presented in Section 4. We conclude in Section 5 with a discussion on the future direction of this research.

2 Problem statement

In this section, we describe a coupled PDE system and its discretization.

2.1 Single phase thermal flow in porous media

We describe the equations for single phase flow in porous media coupled with thermal effects.

2.1.1 Conservation of mass

We start with the continuity equation which states that the rate at which mass enters the system is equal to the rate of mass which leaves the system plus the accumulation of mass within the system. Additionally, we include a source/sink term which accounts for mass which is added or removed from the system. We have

[TABLE]

where $\phi$ is the porosity field of the rock, $\rho$ is the density of the fluid, $\mathbf{u}$ is the fluid velocity, $f$ is a source/sink term, and $\Omega$ is the spatial domain in $\mathbb{R}^{d}$, $d=2,3$. The source/sink term $f$ represents injection/production wells and is given in Section 2.1.5. We further assume that the velocity follows Darcy’s law [16], i.e.

[TABLE]

where $p$ is the pressure, $\mathbf{K}$ is the permeability tensor field, $\mu$ is the viscosity, and $\mathbf{g}$ is gravitational acceleration. The density and viscosity are functions of pressure and temperature given in Section 2.1.4. Then, (1) becomes

[TABLE]

We also assume Neumann and Dirichlet boundary conditions

[TABLE]

where $g_{N}$ is Neumann boundary data, $g_{D}$ is Dirichlet boundary data, $\mathbf{n}$ is the unit outward normal vector on $\partial\Omega=\Gamma_{N}\cup\Gamma_{D}$, and $\Gamma_{D}\cap\Gamma_{N}=\varnothing$.

2.1.2 Conservation of energy

Similarly, we have a conservation of energy equation for the heat energy. Note that formulations where enthalpy is an independent variable are common, but here we consider temperature as an independent variable as in a reference commercial reservoir simulator [17]. Here, $c_{v}$ and $c_{r}$ are the specific heat of the fluid and rock, respectively, $\rho_{r}$ is the density of the rock, and $T$ is temperature. Here, $c_{v}T$ represents the enthalpy of the fluid, and $\rho c_{v}$, its energy density. Heat energy is not only transported by a heat flux, but also by the fluid flux. We get the following advection-diffusion equation:

[TABLE]

where $\mathbf{q}$ is the heat flux, and $f_{T}$ is a source/sink term representing wells or heaters and given in Section 2.1.5. Furthermore, we assume that the heat flux follows Fourier’s law, i.e.

[TABLE]

where $k_{T}$ is the thermal conductivity field. It is given by

[TABLE]

where $k_{T,r}$ and $k_{T,f}$ are the conductivities of the rock and the fluid, respectively. Then, (5) becomes

[TABLE]

Then, assuming Darcy flow, we get

[TABLE]

We also assume Neumann and Dirichlet boundary conditions

[TABLE]

where $g_{N}^{T}$ is Neumann boundary data, $g_{D}^{T}$ is Dirichlet boundary data, $\partial\Omega=\Gamma^{T}_{N}\cup\Gamma^{T}_{D}$, and $\Gamma^{T}_{D}\cap\Gamma^{T}_{N}=\varnothing$.

2.1.3 Coupled problem

We assume that $\rho$ and $\mu$ are empirically determined functions of pressure and temperature. Our choices are given in Section 2.1.4.

We are interested in solving the following boundary value problem:

find $p$, $T$ such that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\partial\Omega=\Gamma_{N}\cup\Gamma_{D}=\Gamma_{N}^{T}\cup\Gamma_{D}^{T}$, $\Gamma_{D}\cap\Gamma_{N}=\Gamma_{D}^{T}\cap\Gamma_{N}^{T}=\varnothing$, and initial conditions for $p$ and $T$ are prescribed.

2.1.4 Nonlinear quantities

The density $\rho$ and viscosity $\mu$ are empirically determined functions of temperature and pressure. For the examples in this paper, we will consider the flow of heavy oil in porous media and thus use the following empirical laws.

For viscosity, we choose the following correlation [18]:

[TABLE]

which takes temperature $T_{\mathrm{F}}$ in ∘F and returns viscosity in cp (0.001 kg m*-1* s*-1*). The viscosity as a function of temperature (in Kelvin) is illustrated in Figure 1. The dimensionless parameters $A_{i}$ can be found in Table 1. The American Petroleum Institute (API) gravity $\gamma_{\mathrm{API}}$ is a measure of how heavy or light a petroleum liquid is compared to water: if its API is greater than 10, than it is lighter and floats on water; if less than 10, it is heavier and sinks. We can calculate API gravity from specific gravity (SG) (ratio of the density of the petroleum liquid to the density of water, at 60∘ F) using the following formula:

[TABLE]

For density, we use the following correlation:

[TABLE]

where $p_{0}$, $T_{0}$ are reference pressure and temperature and $\rho_{0}$ is the density at those values, $c$ is a compressibility coefficient and $\beta$ is a thermal expansion coefficient. Values representative to those used in reservoir simulation are $p_{0}=1.01325$ bar, $T_{0}$= 288.7056 K ($60\;^{\circ}$F), $c=5.5\times 10^{-5}\;\mathrm{bar}^{-1}$, and $\beta=2.5\times 10^{-4}\;\mathrm{K}^{-1}$. Given a specific gravity, we have $\rho_{0}=\gamma_{\mathrm{SG}}\;\rho_{w}$, where $\rho_{w}=999$ kg m*-3* is the density of water at the reference temperature.

2.1.5 Source/sink terms

We first consider source/sink terms representing injection and production wells. A simple way to model these is by using point sources/sinks

[TABLE]

[TABLE]

where $\mathbf{x}_{\mathrm{inj}}$ and $\mathbf{x}_{\mathrm{prod}}$ represent the location of injection and production wells, respectively, $\delta(\mathbf{x})$ is the Dirac delta function, $q_{\mathrm{inj}}^{i}$ and $q_{\mathrm{prod}}^{j}$ are the wells’ injection and production rates, respectively.

The production rate $q_{\mathrm{prod}}$ is usually given by a constant target production rate. Similarly, the injection rate $q_{\mathrm{inj}}$ is given by a target injection rate. These rates can only be maintained if the pressure at the production well does not drop below a minimum pressure, and the pressure at the injection well does not go above a maximum pressure. In those cases, a well model is required. We consider the commonly used Peaceman well model [19, 20] for anisotropic media with $\mathbf{K}=\mathrm{diag}(K_{x},K_{y},K_{z})$ as the permeability tensor field. In this case the rates are given by

[TABLE]

where $h$ is the height of well opening, $K_{e}=\sqrt{K_{x}K_{y}}$ is the equivalent permeability, $p_{bh}$ is the bottom-hole pressure, $r_{w}$ is the well radius, and $r_{e}$ is the equivalent radius which can be calculated using

[TABLE]

where $D_{x}$ and $D_{y}$ are the horizontal lengths of the grid cell. Since we want to allow mesh refinements, we do not want the model to change as we vary the grid size. Therefore, we arbitrarily fix $D_{x}=D_{y}=5$ meters, and also choose $h=5$ meters and $r_{w}=0.1$ meters.

Oil recovery techniques for heavy oils can include electromagnetic heating [21]. These can be expressed as source terms for the energy equation. For simplicity, we do not use an electromagnetic model and choose the simple function

[TABLE]

where $\mathbf{x}_{\mathrm{heater}}$ represent the location of heaters, $U_{\mathrm{heater}}$ is the heat transfer coefficient, and $T_{\mathrm{heater}}$ is the target heating temperature. For our simulations, we have a heating coefficient of $5.44409\times 10^{-6}$ Js*-1K-1*. For simplicity, we also choose $T_{\mathrm{heater}}$ to be the same as $T_{\mathrm{inj}}$.

2.2 DG0 discretization

In reservoir simulation, Finite Volume methods are most commonly used [22]. Since the flux entering a given volume is identical to that leaving an adjacent one, these methods are conservative. Additionally, upwind schemes introduce substantial numerical diffusion, which helps with stability. In this section, we present a discontinuous Galerkin (DG) method [23] that is equivalent to a Finite Volume method used in reservoir simulation and is based on the description in [24]. The resulting weak formulation allows us to implement our problem in the open source Finite Element software Firedrake [25].

Let $\mathcal{T}=\{E_{i},i\in\mathcal{I}\}$ be a partition of $\Omega$ into open element domains $E_{i}$ such that union of their closure is $\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu$, where $\mathcal{I}$ is a set of indices. Let the interior facet $e_{ij}=\mkern 1.5mu\overline{\mkern-1.5muE\mkern-1.5mu}\mkern 1.5mu_{i}\cap\mkern 1.5mu\overline{\mkern-1.5muE\mkern-1.5mu}\mkern 1.5mu_{j}$ and let $\Gamma_{\mathrm{int}}$ denote the union of all interior facets. Let connection set $\mathcal{N}(i)$ denote the set of indices $j$ such that $|e_{ij}|>0$. We begin by presenting a DG0 (piecewise constant) method for the heat equation

[TABLE]

[TABLE]

The variational problem for (23)-(24) on a single cell $E_{i}$ is: find $u$ such that

[TABLE]

where $\mathbf{n}$ is the outward normal to $E_{i}$. Let us first consider $E_{i}$ such that $\partial E_{i}\in\Gamma_{\mathrm{int}}$. Let $h_{i}$ denote the center point of cell $E_{i}$. For the flux on the interior facets, we choose the following flux approximation

[TABLE]

Here, $u\mid_{E_{i}}^{e_{ij}}$ denotes the limit of $u$ in cell $E_{i}$ as it goes to the edge $e_{ij}$.

We consider a piecewise constant approximation of our solution, i.e. in the approximation space $\mathcal{V}_{h}=\mathbb{P}_{\mathrm{DG}}^{0}$ with basis $\{\phi_{i}=\mathbf{1}_{E_{i}}\mid i\in\mathcal{I}\}$. The DG0 approximation is $u_{h}=\sum_{i\in\mathcal{I}}u_{i}\phi_{i}$. For this approximation, on $E_{i}$, $v_{h}\in\mathcal{V}_{h}$ is constant and $\nabla v_{h}=0$. Therefore, (25) becomes

[TABLE]

Note that this is equivalent to

[TABLE]

which is a Finite Volume approximation of the heat equation. In reservoir simulation, this way of approximating the interior facet integrals is known as a “two-point flux” (TPFA) approximation. In order for such a Finite Volume method to converge, the grid must satisfy a certain orthogonality property [26]. In brief, in each cell, there exists a point called the center of the cell such that for any adjacent cell, the straight line between the two centers is orthogonal to the boundary between the cells. For the examples in this paper, we choose quadrilateral meshes, which easily satisfy this condition.

If instead $E_{i}$ is a boundary element, then the boundary integral becomes

[TABLE]

where $\mathbf{n}_{e}$ is the unit outward pointing normal of a cell. We use the following flux approximation

[TABLE]

where $d_{\Gamma_{D}}(h_{i})$ is the shortest distance form $h_{i}$ to the boundary $\Gamma_{D}$. For each $i\in\mathcal{I}$, we have

[TABLE]

For a given ordering of the indices in $\mathcal{I}$, we denote by $u^{+}$ and $u^{-}$ the limit value of $u$ for two cells sharing an edge. Now, summing over all $i\in\mathcal{I}$, and noting that each interior facet is visited twice, we obtain

[TABLE]

We define the jump of $v$ as $[v]=v^{+}-v^{-}$. We then get the following problem: find $u\in\mathbb{P}_{\mathrm{DG}}^{0}$ such that

[TABLE]

for all $v\in\mathbb{P}_{\mathrm{DG}}^{0}$.

2.2.1 Upwinding

We now consider an upwind Godunov method [27] for the advection equation

[TABLE]

[TABLE]

where $\mathbf{w}$ is a given vector field. For an interior $E_{i}$, the upwind scheme is given by

[TABLE]

where $u^{\mathrm{up}}$ is the upwind value of $u$, which, for a facet $e$ shared by $E_{1}$ and $E_{2}$ and $\mathbf{n}_{e}$ pointing from $E_{1}$ to $E_{2}$, is given by

[TABLE]

For the full discretized problem we have: find $u\in\mathbb{P}_{\mathrm{DG}}^{0}$ such that

[TABLE]

for all $v\in\mathbb{P}_{\mathrm{DG}}^{0}$.

2.2.2 Semidiscrete problem

We now discretize (11)-(14) in space using the semidiscrete DG0 formulation described above. Assuming homogeneous Neumann boundary conditions, the variational problem is: find the approximation $(p,T)\in\mathbb{P}_{\mathrm{DG}}^{0}\times\mathbb{P}_{\mathrm{DG}}^{0}$ such that

[TABLE]

[TABLE]

for all $(q,r)\in\mathbb{P}_{\mathrm{DG}}^{0}\times\mathbb{P}_{\mathrm{DG}}^{0}$. The brackets $\{\}$ denote the average across the facets, and the double brackets $\{\!\!\{\}\!\!\}$ denote the harmonic average across the facets. The use of the harmonic average is standard for two-point flux approximation, and is obtained by considering piecewise constant permeabilities [26]. The upwind quantities are given by

[TABLE]

For the delta functions in the source/sink terms, we choose the simple approximation:

[TABLE]

2.2.3 Fully discretized problem

For time discretization, we use the backward Euler method. We define the two following forms, which are linear with respect with their last argument:

[TABLE]

[TABLE]

Let $F(p,T;q,r)\coloneqq F_{m}(p,T;q)+F_{e}(p,T;r)$, which is linear in both $q$ and $r$, but nonlinear in $p$ and $T$. At each time-step, given the previous solution $(p^{n},T^{n})$, we search for $(p^{n+1},T^{n+1})\in\mathbb{P}_{\mathrm{DG}}^{0}\times\mathbb{P}_{\mathrm{DG}}^{0}$ such that

[TABLE]

3 Solution algorithms

The system of nonlinear equations (45) can be written as a system of nonlinear equations for the real coefficients $p_{i}$ and $T_{i}$ of the DG0 functions $p^{n+1}$ and $T^{n+1}$, respectively. Let $x$ be the vector of these coefficients and $G$ the function such that $G(x)=0$ is equivalent to (45). By linearizing this equation with Newton’s method, we must solve at each iteration

[TABLE]

The resulting linearized systems can be written as a block system of the form

[TABLE]

where $\delta x=x_{k+1}-x_{k}$ is the Newton increment. The different blocks are the discrete versions of Jacobian terms as follows

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

and

[TABLE]

All coefficients in (48)-(53) are evaluated at the previous Newton iterate $(p_{k},T_{k})$, and $(.)_{p}$ and $(.)_{T}$ denote the partial derivatives with respect to $p$ and $T$, respectively.

The linearized systems are often very difficult to solve using iterative methods. Indeed, efficient preconditioning is required in order to achieve rapid convergence with linear solvers [28]. In this section, we will detail different preconditioning techniques used to solve (47). We first mention some methods which are important ingredients of the preconditioning techniques.

Krylov subspace methods are used to approximate the solution of $Ax=b$ by constructing a sequence of Krylov subspaces, $\mathcal{K}_{n}=\left\{b,Ab,A^{2}b,\dots,A^{n-1}b\right\}$. The generalized minimal residual method (GMRES) [29] is a Krylov subspace method suitable for general linear systems. The approximate solution $x_{n}$ is formed by minimizing the Euclidean norm of the residual $r_{n}=Ax_{n}-b$ over the subspace $\mathcal{K}_{n}$.

Incomplete LU factorization (ILU) [28, 30] is a general preconditioning technique in which sparse triangular factors are used to approximate the system matrix $A$. This preconditioner requires assembling the factors and then solving two triangular systems. A popular way to determine the sparsity pattern of the factors is to simply choose the relevant triangular parts of the sparsity pattern of $A$. This is known as ILU(0). More generally, choosing the sparsity pattern of $A^{k+1}$ is called ILU($k$).

Multigrid methods [31, 32, 33] use hierarchies of coarse grid approximations in order to solve differential equations. Smoothing operations (such as a Jacobi or Gauss-Seidel iteration) are combined with coarse grid corrections on increasingly coarser grids. For positive definite elliptic PDEs, it is known that multigrid methods can provide optimal solvers (in the sense of linear scalability with the dimension of the discretized problem).

Algebraic Multigrid (AMG) [4, 34] uses information from the entries of the system matrix rather than that of the geometric grid. This makes AMG an ideal black-box solver for elliptic problems. Although it can be used to solve simpler problems, AMG is often used as a preconditioner for Krylov subspace methods in problems which are essentially elliptic. Relative to preconditioners such as ILU, parallel variants of multigrid methods retain more effectiveness.

3.1 Two-stage preconditioning: CPR

Let $M_{1}$ and $M_{2}$ be two preconditioners for the linear system $Ax=b$ for which we have the action of their (generally approximate) inverse $M_{1}^{-1}$, and $M_{2}^{-1}$. Applying a multiplicative two-stage preconditioner can be done as follows:

Precondition using $M_{1}$: $x_{1}=M_{1}^{-1}b$; 2. 2.

Compute the new residual: $b_{1}=b-Ax_{1}$; 3. 3.

Precondition using $M_{2}$ and correct: $x=M_{2}^{-1}b_{1}+x_{1}$.

The action of the two-stage preconditioner can be written as

[TABLE]

In the case of multiphase flow in porous media, the standard preconditioner is the Constrained Pressure Residual method (CPR) [1]. In the multiphase case, the linear systems are like in (47) except that the temperature blocks are replaced (or combined in the thermal case) with saturation blocks.

In the case of CPR, the first stage preconditioner $M_{1}$ is given by

[TABLE]

where $A_{pp}^{-1}$ is approximated using an AMG V-cycle. The second preconditioner is chosen such that $M_{2}^{-1}\approx A^{-1}$, usually with an incomplete LU factorization method (ILU).

In addition to the two-stage preconditioner, decoupling operators are often used to reduce the coupling between the pressure equation and the saturation variables. Indeed, an approximation of the pressure equation $A_{pp}\delta p+A_{pT}\delta T=b_{T}$ is performed in the first stage of CPR where the saturation coupling $A_{pT}$ is ignored. A decoupling operator is a left preconditioner applied a priori to (47) of the form

[TABLE]

The most often-used approximations for multiphase flow are Quasi-IMPES (QI) and True-IMPES (TI)  [3, 35, 36]. The approximations are $D_{QI}=\mathrm{diag}(A_{pT})\mathrm{diag}(A_{TT})^{-1}$, $D_{TI}=\mathrm{colsum}(A_{pT})\mathrm{colsum}(A_{TT})^{-1}$. Here, $\mathrm{colsum}(A)$ is a diagonal matrix with entries the sums of the entries in the columns of $A$, which is equivalent to the mass accumulation terms when discretized with the two-point flux approximation as outlined in this paper (as the fluxes sum up to zero in a given column).

By performing this decoupling operation on the system (47) before CPR, the first stage now consists in solving a subsystem for the approximate Schur complement $S_{p}=A_{pp}-DA_{Tp}$ instead of the original pressure block. However, the properties of the resulting $S_{p}$ need to be amenable to the application of AMG (for example M-matrix properties). While this is nearly guaranteed in the black-oil case [37], it does not necessarily follow for compositional flow or thermal flow. For the single phase test cases detailed in Section 4, we observe that CPR performs best without decoupling operators (results not shown here).

3.2 Block factorization preconditioner

Consider the following decomposition of the Jacobian

[TABLE]

where $S_{T}=A_{TT}-A_{Tp}A_{pp}^{-1}A_{pT}$ is the Schur complement. The inverse of the Jacobian is given by

[TABLE]

Even if $A$ is sparse, the Schur complement $S_{T}$ is generally dense. A common preconditioning technique is to use the blocks of the factorization (57) combined with a sparse approximation of the Schur complement [38]. Given an appropriate Schur complement approximation $\tilde{S}_{T}$, applying the block preconditioner can be done as follows:

Solve the pressure subsystem: $A_{pp}x_{p}=b_{p}$; 2. 2.

Compute the new energy equation residual: $\tilde{b}_{T}=b_{T}-A_{Tp}x_{p}$; 3. 3.

Solve the Schur complement subsystem: $\tilde{S}_{T}\delta_{T}=\tilde{b}_{T}$; 4. 4.

Compute the new mass equation residual: $\tilde{b}_{p}=x_{p}-A_{pT}\delta_{T}$; 5. 5.

Solve the pressure subsystem: $A_{pp}\delta_{p}=\tilde{b}_{p}$.

In our case, $A_{pp}^{-1}$ and $\tilde{S}_{T}^{-1}$ are both approximated using an AMG V-cycle.

3.3 Schur complement approximation

Common sparse approximations for the Schur complement are ${\tilde{S}_{A_{TT}}}=A_{TT}$ and $\tilde{S}_{\mathrm{diag}}=A_{TT}-A_{Tp}\mathrm{diag}\left(A_{pp}\right)^{-1}A_{pT}$. Here we present a Schur complement approximation which performs significantly better than such simple approximations.

For the derivation of our Schur complement approximation, we consider the linearized problem before discretization. This approach results in an approximation which holds as we refine the mesh. See [39] for a theoretical framework in using the infinite-dimensional setting to find mesh-independent preconditioners for self-adjoint problems.

3.3.1 Steady-state case

We first consider a steady-state single phase thermal problem: find $p$, $T$ such that

[TABLE]

[TABLE]

where $\mathbf{u}$ is given by (2), and we have no-flux boundary conditions for the fluid and heat. Here we will consider the linearized system in a continuous setting. Applying a Newton method to (59)-(60), we obtain a block systems of the form (47) where the blocks are:

[TABLE]

[TABLE]

[TABLE]

where we have used the product rule for the divergence operator in (62) and (63). Then the second term of the Schur complement (which corresponds in the continuous setting to the Poincaré-Steklov operator) becomes

[TABLE]

We notice that in $A_{TT}-A_{Tp}A_{pp}^{-1}A_{pT}$, the terms $c_{v}T\nabla\cdot(\rho\mathbf{u})_{T}$ cancel. We are left with

[TABLE]

One of the nonlinear terms has canceled, and so we consider if it is possible that the last two terms also cancel. Consider the operator $(\rho\mathbf{u})_{p}\left(\nabla\cdot(\rho\mathbf{u})_{p}\right)^{-1}\nabla\cdot$, which is close to the operator $\nabla(\nabla\cdot\nabla)^{-1}\nabla\cdot\eqqcolon s$. This holds if $(\rho\mathbf{u})_{p}$ is close to $\nabla$, i.e. if $\rho$ is close to being constant with respect to $p$. While this approximation holds for liquid water and hydrocarbons, it may be less applicable in the case of gases. Extending this to multiphase flow is straightforward and is part of ongoing work.

Assuming that the operator $s$ is applied to a sufficiently smooth vector field $\mathbf{F}$, we can use Helmholtz decomposition to decompose this field into the sum of its curl-free and divergence-free part $\mathbf{F}=-\nabla\Phi+\nabla\times\mathbf{A}$, where $\Phi$ is a scalar potential and $\mathbf{A}$ a vector potential. Since $s$ removes the divergence-free part of a field, applying $s$ to $\mathbf{F}$, we obtain

[TABLE]

assuming that $\Phi$ satisfies the same boundary conditions as the operator $(\nabla\cdot\nabla)^{-1}$. Hence $s$ is a projection to the curl-free subspace, and it acts like the identity operator when applied to curl-free vector fields.

We assume that this also holds for $(\rho\mathbf{u})_{p}\left(\nabla\cdot(\rho\mathbf{u})_{p}\right)^{-1}\nabla\cdot$. In (65), we see that this operator is applied to $(\rho\mathbf{u})_{T}$ which is of the form $\gamma\nabla p$, where $\gamma$ is a scalar field. In order for $\gamma\nabla p$ to be a curl-free vector field, we need $\nabla\times(\gamma\nabla p)=\nabla\gamma\times\nabla p=0$, i.e. we need $\nabla\gamma$ and $\nabla p$ to be parallel vectors. In the discretized case, our grid satisfies an orthogonality property as mentioned in Section 2.2, and the gradients of $p$ and $\gamma$ are approximated using a two-point flux approximation. In this case, the gradients are always orthogonal to the facets, and thus parallel. Accordingly, we replace the operator $(\rho\mathbf{u})_{p}\left(\nabla\cdot(\rho\mathbf{u})_{p}\right)^{-1}\nabla\cdot$ by the identity and obtain the following Schur complement approximation

[TABLE]

Similar heuristic arguments for replacing $\nabla\cdot(\nabla\cdot\nabla)^{-1}\nabla$ by the identity operator in the case of the Stokes problem can be found, for example, in [40], and for the Navier-Stokes equations, in [41].

3.3.2 Source terms

Similarly, we consider the steady-state case with the addition of source/sink terms. In this case, production wells satisfy $f^{\mathrm{prod}}_{T}=c_{v}Tf_{\mathrm{prod}}$, while injection wells satisfy $f^{\mathrm{inj}}_{T}=c_{v}T_{\mathrm{inj}}f_{\mathrm{inj}}$. Thus,

[TABLE]

Since the injection term is weighted by $T_{\mathrm{inj}}$, we cannot directly cancel the $c_{v}T$ terms as in (65). However, after a certain amount of injection, $T$ tends to $T_{\mathrm{inj}}$ where the injection well is located. Furthermore, in the infinite-dimensional setting, this effect will be instantaneous since the well terms are defined using a delta function in (19). Using this argument, we get

[TABLE]

Further assuming that the mass source/sink terms are almost constant in $p$ and $T$, i.e. $(\rho)_{T}$ and $(\rho)_{p}$ are small and the injection/production rates are independent of pressure and temperature (which is the case when operating at a target rate), we ignore the derivatives of the source/sink terms. Then, using the same argument as for the steady-state case, we obtain the Schur complement approximation

[TABLE]

In the case where the the source terms are heaters, we have $f=0$, and $f_{T}=U(T_{\mathrm{heater}}-T)D_{\mathrm{heaters}}$, where $D_{\mathrm{heaters}}$ is the sum of delta functions for the location of heaters. The Schur complement is given by

[TABLE]

We see that heaters do not affect the right-hand side term. Using the same argument as above, we get the approximation

[TABLE]

3.3.3 Time-dependent case

We now generalize our analysis to the time-dependent problem. We first consider the case without source/sink terms. The blocks are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The second term of the Schur complement is given by

[TABLE]

and thus the Schur complement is

[TABLE]

To justify further simplification, we need to assume that either $\rho$ is almost constant in $p$ and $T$, or that the time-step is very large. We get the following Schur complement approximation:

[TABLE]

In the case where we have source/sink terms, the Schur complement is given by

[TABLE]

Using the same argument for the injection temperature $T_{\mathrm{inj}}$ as in Section 3.3.2, we get

[TABLE]

Then, again assuming that the mass source/sink terms are independent of $p$ and $T$, we obtain the following Schur complement approximation:

[TABLE]

We can obtain the discretized version of this operator from (44) by removing the terms depending on the previous time-step, and evaluating the nonlinear terms at the previous Newton iteration. We get the following bilinear operator:

[TABLE]

4 Numerical results

In this section, we perform numerical experiments for our block preconditioner and CPR. These are implemented on the open source Finite Element software Firedrake [25]. The linear algebra backend is the PETSc library [42], allowing efficient and parallel computations. The CPR preconditioner (without decoupling) is implemented by providing PETSc options. Our custom block preconditioner is implemented through Firedrake’s Python interface. Recent work from [43] allows us to easily assemble our Schur complement approximation preconditioner by providing a weak form with the bilinear operator (82). We modified the custom preconditioner class from [43] to allow the use of matrix formats other than matfree. For example, the default aij matrix format allows for faster computations for lower order methods such as the one described in Section 2.2. Our implementation is available on GitHub111https://github.com/tlroy/thermalporous.

For the block preconditioner, we use our Schur complement approximation (82), unless stated otherwise. Both the pressure block $A_{pp}$ and the approximate Schur complement are inverted using a V-cycle of AMG. We use BoomerAMG [8] from the hypre library [9] with default parameters, i.e. a symmetric-SOR/Jacobi relaxation scheme (one sweep up, one sweep down), Falgout coarsening, classical Ruge-Stüben interpolation, and Gaussian Elimination as the coarse grid solver. This implementation has a very efficient parallel version of AMG.

For the second stage of CPR, we use ILU(0) as provided from PETSc. In parallel, we use block Jacobi with ILU(0) for each block (the partition is assigned when Firedrake does the discretization).

The nonlinear solver is Newton’s method with line search, and the linear solver is right-preconditioned GMRES [29], restarted after 30 iterations. The convergence tolerance of Newton’s methods is $10^{-8}$ for the relative function norm and relative step size norm. The convergence tolerance for GMRES is $10^{-10}$ for the relative residual norm for the tests in Section 4.1, and $10^{-5}$ for the tests in Sections 4.3 and 4.4.

For all cases, we consider a heavy oil with density and viscosity as described in Section 2.1.4. The other physical parameters are shown in Table 2. These parameters are representative with those used in commercial reservoir simulators.

For all cases, we evaluate the performance of the methods by comparing the number of linear iterations per nonlinear iteration. We note that, for our proof-of-concept implementation, the cost of applying the block preconditioner is around two times more computationally expensive (in serial) than CPR. The difference may not be as significant in an optimized implementation.

4.1 SPE10 test cases

The domain is a square with dimensions 365.76$\times$365.76 meters, and the mesh is 60$\times$120. For permeability and porosity fields, we use the benchmark problem SPE10 [44]. This problem has a highly heterogeneous permeability field. We consider a 60$\times$120 slice in the xy direction. The permeability, which is isotropic in the xy plane, is illustrated in Figure 2. We do not include gravity for the 2D simulations.

For the well case (W), we have one production well and one injection well. These are located in the upper half of the domain in the regions of high permeability. For the injection and production rates, we use the Peaceman well model. The bottom-hole pressure for the injection well is fixed at $6.895\times 10^{7}$ Pa, and $2.7579\times 10^{7}$ Pa for the production well. The maximum rate is set to $q=1.8\times 10^{-3}$m3s*-1*, although this is only achieved for the high permeability cases. The initial temperature in the reservoir is 288.706 K and the injection temperature is 422.039 K. For the heater case (H), heater placement is the same as for the well case, and so are the initial and heating temperatures. For the well and heater case (W+H), we combine both wells and heaters. For the high permeability cases (h.p.), we increase the permeability by a factor of 1,000. While the resulting permeability values are not representative of physical ones, they give a simple example of advection-dominated heat flow.

For each case, we simulate injection and production for 1000 days where the time steps are chosen adaptively such that Newton’s method converges in around 4 iterations. The average linear iterations per nonlinear iteration are shown in Table 3.

We observe that for the first three cases in Table 3, the block preconditioner performs better than CPR in terms of the number of GMRES iterations, but that CPR performs best for the high permeability cases. The heat flow for the first three cases is diffusion-dominated, especially when the oil is not yet heated. For the high permeability cases, advection dominates. This change in performance appears later in the simulation when temperature has increased everywhere between the two wells. This indicates that CPR can still be a good choice if temperature is simply transported by the fluid flow. However, we will see in the next sections that the block preconditioner is a more scalable method.

4.2 Numerical justification of the Schur complement approximation

We now perform a numerical comparison of the action of the inverses of the different Schur complement approximations. We use the cases given in Section 4.1. In Table 4, we compare the different Schur complement approximations by looking at the condition number of their inverse applied to the full Schur complement. While this condition number does not directly inform us about how well the preconditioner performs, it is a good indication of the quality of the approximations. For the cases, H and W stand for heaters and wells, respectively, and h.p. stands for high permeability (increased by a factor 1,000). We observe that $\tilde{S}_{T}$ is a good Schur complement approximation even for the high permeability cases where the other approximations struggle.

In terms of the performance of the solver, $\tilde{S}_{T}$ always results in fewer GMRES iterations (results not shown here). For harder cases (for example high permeability), this difference is significant; the linear solver can even fail to converge before the prescribed maximum number of iterations. In the next section, we will see that the other Schur complement approximations struggle in anisotropic medium.

4.3 Problem size scaling

We now investigate the performance of CPR and our block preconditioner as we refine a mesh. For two cases, we will also consider the Schur complement approximations $\tilde{S}_{A_{TT}}$ and $\tilde{S}_{\mathrm{diag}}$. To this end, we test a case with homogeneous permeability and porosity fields. The domain is a square with dimensions 20$\times$ 20 meters and uniform porosity $\phi=0.2$. We test both isotropic and anisotropic permeability fields. We refine the mesh from a $20\times 20$ grid to $320\times 320$.

We begin with an isotropic permeability of $3\times 10^{-13}$ m2. For all cases, the injection/heating temperature is 422.039K. For all cases except Case III, the initial temperature is 288.706K. For each case, we take two time steps and calculate the average number of linear iterations per nonlinear iteration. For Case I-IV, the time step is 10 days, and for Case V, 12 hours.

For Case I, we have 6 heaters in the domain. In Table 5, we observe that the number iteration increases by 9 times for CPR, while it increases by less than 50% for the block preconditioner.

For Case II, we have injection wells and 3 production wells. The wells operate at constant injection and production rates of $q=5\times 10^{-8}m^{3}s^{-1}$. In Table 6, we observe that the number of iterations for CPR increases by 10 times while it only increases by less than 50 % for the block preconditioner with the Schur complement approximation $\tilde{S}_{T}$. We observe a similar increase in iterations for the block preconditioner with the Schur complement approximations $\tilde{S}_{A_{TT}}$ and $\tilde{S}_{\mathrm{diag}}$.

For Case III, we also have 3 injection wells and 3 production wells. To allow higher rates and faster flow, we increase the initial temperature to 320K. The wells operate at injection and production rates $q=10^{-6}m^{3}s^{-1}$. In Table 6, we observe that the number of iterations for CPR increases by more than 10 times while it only increases by around 50 % for the block preconditioner.

For case IV and V, we increase the permeability in the $x$-direction to $3\times 10^{-11}$m2. For Case IV, we have 6 heaters and observe the same trend as the previous case in Table 8.

For Case V, we have 3 injection wells and 3 production wells. The wells operate at constant injection and production rate $q=1\times 10^{-6}m^{3}s^{-1}$. In this case, the flow is much faster and thus the time step size is reduced to half a day for the convergence of Newton’s method. In Table 9, we observe that the number of iterations is doubled for the block preconditioner with $\tilde{S}_{T}$, increased by 4 times for CPR, and slightly less for the block preconditioner with $\tilde{S}_{A_{TT}}$. Additionally, the block preconditioner with $\tilde{S}_{\mathrm{diag}}$ fails to converge within 200 GMRES iterations.

In summary, as we refine the mesh, the number of iterations has a very small increase for the block preconditioner, but a large increase for CPR. The heat diffusion is much more noticeable on fine meshes, which CPR does not treat appropriately. However, coarser meshes are more common in commercial reservoir simulators.

Note that the success of the block preconditioner is also due the linear scalability of AMG for elliptic problems. By removing the need for ILU, we get a nearly mesh-independent preconditioner.

4.4 Parallel scaling

We now compare the performance of the two methods in parallel. We look at both weak and strong scaling.

4.4.1 Weak scaling

For weak scaling, we compare the parallel performance of the methods as we increase the number of processors and problem size. The domain is $50\times 50\times 50$ meters with an $N\times N\times N$ grid. Since this is a 3D case, we include gravity. The permeability is $3\times 10^{-13}$m2 and the porosity is 0.2. We seek to have around 100,000 degrees of freedom per processor. Thus, for the number of processors 1, 2, 4, 8, and 16, we have $N$ = 36, 46, 58, 73, 92.

For the heating case, we have two sets of 21 heaters near the top and bottom of the domain. We take two time steps of 100 days and illustrate the results in Table 10. We observe that the number of iterations increases by around 20 % for the block preconditioner and triples for CPR.

For the well case, we have 21 injection wells near the top of the domain, and 21 production wells near the bottom. All wells operate at a constant injection/production rate $q=10^{-7}m^{3}s^{-1}$. We take two time steps of 10 days and illustrate the performance of the methods in Table 11. We observe that the number of iterations for the block preconditioner increases by around 40% while the number of iterations for CPR nearly triples.

4.4.2 Strong scaling

We use the same problem as the previous section on the finest mesh. We keep the problem size fixed while increasing the number of processors. For reservoir simulation, strong scaling is usually more relevant than weak scaling. Indeed, reservoir models often come with a (usually rather coarse) fixed grid. As observed in Section 4.3, CPR does not behave as well on a fine mesh. Thus, keeping the mesh size constant is a good way of isolating the parallel performance of the methods.

In Tables 12 and 13, we illustrate the strong scaling results for the heating and well cases, respectively. For the block preconditioner, we observe that the number of iteration is essentially independent of the number of processors used. This is thanks to the parallel capability of BoomerAMG. On the other hand, the number of iterations for CPR exhibit a small but progressive increase. This is because the second stage of CPR uses Block ILU, which becomes a weaker preconditioner as the number of blocks increases. Therefore, this trend will continue as the number of processors increases.

5 Conclusion

In this work, we have implemented a fully implicit parallel non-isothermal porous media flow simulator including two preconditioning strategies, CPR and a block preconditioner with our own Schur complement approximation. We have tested the performance of these methods as preconditioners for GMRES. Our Schur complement approximation performs better than simple one, especially in cases with heterogeneous or anisotropic permeability. While the block preconditioner performs well for diffusion-dominated cases, CPR is still the method of choice for advection-dominated (manufactured) cases, at least in serial. However, the block preconditioner scales optimally with problem size while CPR does not do well under mesh refinement. Additionally, the block preconditioner remains efficient in parallel, while the CPR iteration count increases gradually as we increase the number of processors.

This research demonstrates that a preconditioning strategy which considers the diffusive effect of temperature is important for diffusion-dominated cases. In non-isothermal multiphase flow, the energy equation is treated in CPR like a hyperbolic equation. A coupled solution of pressure and temperature using multigrid is key to methods for multiphase flow currently being developed.

6 Acknowledgments

This publication is based on work partially supported by the EPSRC Centre For Doctoral Training in Industrially Focused Mathematical Modelling (EP/L015803/1) in collaboration with Schlumberger.