A block preconditioner for non-isothermal flow in porous media

# A block preconditioner for non-isothermal flow in porous media

Thomas Roy Tom B. Jönsthövel Christopher Lemon Andrew J. Wathen Mathematical Institute, University of Oxford, Oxford, United Kingdom Schlumberger Houston Production Technology Center, Houston, TX, United States of America Schlumberger Abingdon Technology Center, Abingdon, United Kingdom
###### 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.

###### keywords:
preconditioning, iterative solvers, porous media, thermal reservoir simulation
journal: 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  wallis1983incomplete (); wallis1985constrained () 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 lacroix2003iterative () with the use of Algebraic Multigrid (AMG) ruge1987algebraic () 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 clees2005amg (), often called System AMG (SAMG). While it seeks to replace CPR-AMG, the SAMG approach proposed in clees2010efficient (); gries2014preconditioning () 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.

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 li2017enhanced (). 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) booth2008miscible (). We present a block preconditioner for the solution of the resulting coupled pressure-temperature system. Additionally, this solver can be extended to the multiphase case as the first stage of a CPR-like two-stage preconditioner roy2019constrained ().

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

 ϕ∂ρ∂t+∇⋅(ρu)=fin% R+×Ω, (1)

where is the porosity field of the rock, is the density of the fluid, is the fluid velocity, is a source/sink term, and is the spatial domain in , . The source/sink term represents injection/production wells and is given in Section 2.1.5. We further assume that the velocity follows Darcy’s law darcy1856fontaines (), i.e.

 u=−Kμ(∇p−ρg), (2)

where is the pressure, is the permeability tensor field, is the viscosity, and is gravitational acceleration. The density and viscosity are functions of pressure and temperature given in Section 2.1.4. Then, (1) becomes

 ϕ∂ρ∂t−∇⋅(ρKμ(∇p−ρg))=fin R+×Ω. (3)

We also assume Neumann and Dirichlet boundary conditions

 −Kμ(∇p−ρg)⋅n=gN on% ΓN,andp=gD on ΓD, (4)

where is Neumann boundary data, is Dirichlet boundary data, is the unit outward normal vector on , and .

#### 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 debaun2005extensible (). Here, and are the specific heat of the fluid and rock, respectively, is the density of the rock, and is temperature. Here, represents the enthalpy of the fluid, and , 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:

 ϕ∂∂t(ρcvT)+(1−ϕ)∂∂t(ρrcrT)+∇⋅(ρcvTu)+∇⋅q=fTin R+×Ω, (5)

where is the heat flux, and 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.

 q=−kT∇T, (6)

where is the thermal conductivity field. It is given by

 kT=ϕkT,r+(1−ϕ)kT,f, (7)

where and are the conductivities of the rock and the fluid, respectively. Then, (5) becomes

 ϕ∂∂t(ρcvT)+(1−ϕ)∂∂t(ρrcrT)+∇⋅(ρcvTu)−∇⋅(kT∇T)=fTin R+×Ω. (8)

Then, assuming Darcy flow, we get

 ϕ∂∂t(ρcvT)+(1−ϕ)∂∂t(ρrcrT)−∇⋅(ρcvTKμ(∇p−ρg))−∇⋅(kT∇T)=fTin R+×Ω. (9)

We also assume Neumann and Dirichlet boundary conditions

 −(ρcvTKμ(∇p−ρg)+kT∇T)⋅n=gTN on ΓTN,andT=gTD on ΓTD, (10)

where is Neumann boundary data, is Dirichlet boundary data, , and .

#### 2.1.3 Coupled problem

We assume that and 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 , such that

 ϕ∂ρ∂t−∇⋅(ρKμ(∇p−ρg))=fin R+×Ω, (11)
 ϕ∂∂t(ρcvT)+(1−ϕ)∂∂t(ρrcrT)−∇⋅(ρcvTKμ(∇p−ρg))−∇⋅(kT∇T)=fTin R+×Ω, (12)
 −Kμ(∇p−ρg)⋅n=gN on% ΓN,andp=gD on ΓD, (13)
 −(ρcvTKμ(∇p−ρg)+kT∇T)⋅n=gTN on ΓTN,andT=gTD on ΓTD, (14)

where , , and initial conditions for and are prescribed.

#### 2.1.4 Nonlinear quantities

The density and viscosity 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 bennison1998prediction ():

 μ(TF)=10A1γAPI+A2TA3γAPI+A4F, (15)

which takes temperature in F and returns viscosity in cp (0.001 kg m s). The viscosity as a function of temperature (in Kelvin) is illustrated in Figure 1. The dimensionless parameters can be found in Table 1. The American Petroleum Institute (API) gravity 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:

 γAPI=141.5γSG−131.5. (16)

For density, we use the following correlation:

 ρ(p,T)=ρ0ec(p−p0)eβ(T−T0), (17)

where , are reference pressure and temperature and is the density at those values, is a compressibility coefficient and is a thermal expansion coefficient. Values representative to those used in reservoir simulation are bar, = 288.7056 K (F), , and . Given a specific gravity, we have , where kg m 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

 f(x)= ∑iqiinj(p,T)δ(x−xiinj)ρ(p,Tinj) −∑jqjprod(p,T)δ(x−xjprod)ρ(p,T), (18)
 fT(x)= ∑iqiinj(p,T)δ(x−xiinj)ρ(p,Tinj)cvTinj −∑jqjprod(p,T)δ(x−xjprod)ρ(p,T)cvT, (19)

where and represent the location of injection and production wells, respectively, is the Dirac delta function, and are the wells’ injection and production rates, respectively.

The production rate is usually given by a constant target production rate. Similarly, the injection rate 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 peaceman1978interpretation (); chen2009well () for anisotropic media with as the permeability tensor field. In this case the rates are given by

 q=2πhKeμln(re/rw)(pbh−p), (20)

where is the height of well opening, is the equivalent permeability, is the bottom-hole pressure, is the well radius, and is the equivalent radius which can be calculated using

 re=0.14((Ky/Kx)1/2D2x+(Kx/Ky)1/2D2y))1/20.5((Ky/Kx)1/4+(Kx/Ky)1/4), (21)

where and 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 meters, and also choose meters and meters.

Oil recovery techniques for heavy oils can include electromagnetic heating sahni2000electromagnetic (). 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

 fT=∑iUheater(p,T)δ(x−xiheater)(Theater−T), (22)

where represent the location of heaters, is the heat transfer coefficient, and is the target heating temperature. For our simulations, we have a heating coefficient of JsK. For simplicity, we also choose to be the same as .

### 2.2 DG0 discretization

In reservoir simulation, Finite Volume methods are most commonly used leveque2002finite (). 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 riviere2008discontinuous () that is equivalent to a Finite Volume method used in reservoir simulation and is based on the description in riseth2015nonlinear (). The resulting weak formulation allows us to implement our problem in the open source Finite Element software Firedrake rathgeber2016firedrake ().

Let be a partition of into open element domains such that union of their closure is , where is a set of indices. Let the interior facet and let denote the union of all interior facets. Let connection set denote the set of indices such that . We begin by presenting a DG0 (piecewise constant) method for the heat equation

 ∂u∂t−∇2u=0 in Ω, (23)
 u=f on ΓD,∇u⋅n=g on ΓN. (24)

The variational problem for (23)-(24) on a single cell is: find such that

 ∫Ei∂u∂tvdx+∫Ei∇u⋅∇vdx−∫∂Eiv∇u⋅nds=0for all test functions v, (25)

where is the outward normal to . Let us first consider such that . Let denote the center point of cell . For the flux on the interior facets, we choose the following flux approximation

 ∫∂Eiv∇u⋅nds≈∑j∈N(i)∫eijv∣eijEiu∣eijEj−u∣eijEi∥hj−hi∥ds. (26)

Here, denotes the limit of in cell as it goes to the edge .

We consider a piecewise constant approximation of our solution, i.e. in the approximation space with basis . The DG0 approximation is . For this approximation, on , is constant and . Therefore, (25) becomes

 ∫Ei∂ui∂tvidx−∑j∈N(i)∫eijviuj−ui∥hj−hi∥ds=0. (27)

Note that this is equivalent to

 ∂ui∂t|Ei|−∑j∈N(i)uj−ui∥hj−hi∥|eij|=0, (28)

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 eymard2000finite (). 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 is a boundary element, then the boundary integral becomes

 ∫∂Eiv∇u⋅nds=∑j∈N(i)∫eijv∇u⋅ne+∫∂Ei∩ΓDv∇u⋅nds+∫∂Ei∩ΓNv∇u⋅nds, (29)

where is the unit outward pointing normal of a cell. We use the following flux approximation

 ∫∂Eiv∇u⋅nds\coloneqq ∑j∈N(i)∫eijviuj−ui∥hj−hi∥ds +∫∂Ei∩ΓDvi(f−ui)dΓD(hi)ds+∫∂Ei∩ΓNvigds, (30)

where is the shortest distance form to the boundary . For each , we have

 ∫Ei∂ui∂tvidx−∑j∈N(i)∫eijviuj−ui∥hj−hi∥ds−∫∂Ei∩ΓDvi(f−ui)dΓD(hi)ds −∫∂Ei∩ΓNvigds=0. (31)

For a given ordering of the indices in , we denote by and the limit value of for two cells sharing an edge. Now, summing over all , and noting that each interior facet is visited twice, we obtain

 ∫Ω∂u∂tvdx+∫Γint(v+−v−)u+−u−∥h+−h−∥dS−∫ΓDvf−udΓD(h)ds−∫ΓNvgds=0. (32)

We define the jump of as . We then get the following problem: find such that

 ∫Ω∂u∂tvdx+∫Γint[v][u]∥h+−h−∥dS−∫ΓDvf−udΓD(h)ds−∫ΓNvgds=0, (33)

for all .

#### 2.2.1 Upwinding

We now consider an upwind Godunov method godunov1959difference () for the advection equation

 ∂u∂t+∇⋅(uw)=0 on Ω, (34)
 u=f on ΓD,uw⋅n=g on ΓN, (35)

where is a given vector field. For an interior , the upwind scheme is given by

 ∫Ei∂u∂tvdx+∫∂Eivuupw⋅nds=0, (36)

where is the upwind value of , which, for a facet shared by and and pointing from to , is given by

 uup={u∣eE1if w⋅ne≥0,u∣eE2if w⋅ne<0. (37)

For the full discretized problem we have: find such that

 ∫Ω∂u∂tvdx+∫Γint[v]uupw⋅nds−∫ΓDvf−udΓD(h)ds−∫ΓNvgds=0, (38)

for all .

#### 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 such that

 ϕ∫Ω∂ρ∂tqdx+∫Γint[q]({{K}}ρupμup([p]∥h+−h−∥−{ρ}g⋅ne))dS −∫Ωfqdx=0, (39)
 ϕ∫Ωcv∂ρT∂trdx+(1−ϕ)∫Ωρrcr∂T∂trdx+∫Γint[r]{{K}}cvρupμupTup([p]∥h+−h−∥−{ρ}g⋅ne)dS+∫Γint[r]{{kT}}[T]∥h+−h−∥dS−∫ΩfTrdx=0, (40)

for all . 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 eymard2000finite (). The upwind quantities are given by

 (u)up=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩u∣eE1if [p]∥h+−h−∥−{ρ}g⋅ne≥0,u∣eE2if [p]∥h+−h−∥−{ρ}g⋅ne<0. (41)

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

 δ(x)={1/|Ei|if x∈Ei,0otherwise. (42)

#### 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:

 (43)
 Fe(pn+1,Tn+1;r)\coloneqqϕ∫Ωcvρn+1Tn+1−ρnTnΔtrdx+(1−ϕ)∫ΩρrcrTn+1−TnΔtrdx+∫Γint[r]{{K}}(ρn+1)up(μn+1)up(Tn+1)up([pn+1]∥h+−h−∥−{ρn+1}g⋅ne)dS+∫Γint[r]{{kT}}[Tn+1]∥h+−h−∥dS−∫Ωfn+1Trdx. (44)

Let , which is linear in both and , but nonlinear in and . At each time-step, given the previous solution , we search for such that

 F(pn+1,Tn+1;q,r)=0for all (q,r)∈P0DG×P0DG. (45)

## 3 Solution algorithms

The system of nonlinear equations (45) can be written as a system of nonlinear equations for the real coefficients and of the DG0 functions and , respectively. Let be the vector of these coefficients and the function such that is equivalent to (45). By linearizing this equation with Newton’s method, we must solve at each iteration

 ∂G∂x∣x=xk(xk+1−xk)=−G(xk). (46)

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

 Aδx=[AppApTATpATT][δpδT]=[bpbT]=b, (47)

where is the Newton increment. The different blocks are the discrete versions of Jacobian terms as follows

 App∼ϕ1Δt(ρ)p+∇⋅(ρu)p−(f)p, (48)
 ApT∼ϕ1Δt(ρ)T+∇⋅(ρu)T−(f)T, (49)
 ATp∼ϕ1Δt(ρ)pcvT+∇⋅(cvT(ρu)p)−(fT)p, (50)
 ATT∼ϕcv(ρ+(ρ)TT)Δt+(1−ϕ)ρrcrΔt+∇⋅(cvρu)+∇⋅(cvT(ρu)T)−∇⋅(kT∇)−(fT)T, (51)

where

 (ρu)p=−Kμ(ρ(∇−(ρ)pg)+(ρ)p(∇p−ρg)), (52)

and

 (ρu)T=−K(ρμ)T(∇p−ρg). (53)

All coefficients in (48)-(53) are evaluated at the previous Newton iterate , and and denote the partial derivatives with respect to and , 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 saad2003iterative (). 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 by constructing a sequence of Krylov subspaces, . The generalized minimal residual method (GMRES) saad1986gmres () is a Krylov subspace method suitable for general linear systems. The approximate solution is formed by minimizing the Euclidean norm of the residual over the subspace .

Incomplete LU factorization (ILU) saad2003iterative (); meijerink1977iterative () is a general preconditioning technique in which sparse triangular factors are used to approximate the system matrix . 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 . This is known as ILU(0). More generally, choosing the sparsity pattern of is called ILU().

Multigrid methods brandt1977multi (); briggs2000multigrid (); trottenberg2000multigrid () 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) ruge1987algebraic (); stuben2001introduction () 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 and be two preconditioners for the linear system for which we have the action of their (generally approximate) inverse , and . Applying a multiplicative two-stage preconditioner can be done as follows:

1. Precondition using : ;

2. Compute the new residual: ;

3. Precondition using and correct: .

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

 M−1=M−12(I−AM−11)+M−11. (54)

In the case of multiphase flow in porous media, the standard preconditioner is the Constrained Pressure Residual method (CPR) wallis1983incomplete (). 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 is given by

 M−11≈[A−1pp000], (55)

where is approximated using an AMG V-cycle. The second preconditioner is chosen such that , 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 is performed in the first stage of CPR where the saturation coupling is ignored. A decoupling operator is a left preconditioner applied a priori to (47) of the form

 [I−D0I]. (56)

The most often-used approximations for multiphase flow are Quasi-IMPES (QI) and True-IMPES (TI)  lacroix2003iterative (); lacroix2000iterative (); scheichl2003decoupling (). The approximations are , . Here, is a diagonal matrix with entries the sums of the entries in the columns of , 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 instead of the original pressure block. However, the properties of the resulting need to be amenable to the application of AMG (for example M-matrix properties). While this is nearly guaranteed in the black-oil case gries2015system (), 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

 A=[I0ATpA−1ppI][App00ST][IA−1ppApT0I], (57)

where is the Schur complement. The inverse of the Jacobian is given by

 A−1=[I−A−1ppApT0I][A−1pp00S−1T][I0−ATpA−1ppI]. (58)

Even if is sparse, the Schur complement 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 elman2014finite ().

### 3.3 Schur complement approximation

Common sparse approximations for the Schur complement are and . 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 mardal2011preconditioning () 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 , such that

 ∇⋅(ρu)=0in Ω, (59)
 ∇⋅(ρcvTu)−∇⋅(kT∇T)=0in Ω, (60)

where 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:

 App=∇⋅(ρu)p,ApT=∇⋅(ρu)T, (61)
 ATp=∇⋅(cvT(ρu)p)=cv[∇T⋅(ρu)p+T∇⋅(ρu)p], (62)
 ATT =cv[∇⋅(ρu)+∇⋅(T(ρu)T)]−∇⋅(kT∇) (63)

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

 ATpA−1ppApT =cv[∇T⋅(ρu)p+T∇⋅(ρu)p](∇⋅(ρu)p)−1∇⋅(ρu)T =cvT∇⋅(ρu)T+cv∇T⋅(ρu)p(∇⋅(ρu)p)−1∇⋅(ρu)T. (64)

We notice that in , the terms cancel. We are left with

 ST=cv∇⋅(ρu)+cv∇T⋅(ρu)T−∇⋅(kT∇)+cv∇T⋅(ρu)T−cv∇T⋅(ρu)p(∇⋅(ρu)p)−1∇⋅(ρu)T. (65)

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 , which is close to the operator . This holds if is close to , i.e. if is close to being constant with respect to .

Assuming that the operator is applied to a sufficiently smooth vector field , we can use Helmholtz decomposition to decompose this field into the sum of its curl-free and divergence-free part , where is a scalar potential and a vector potential. Since removes the divergence-free part of a field, applying to , we obtain

 sF=−∇(∇⋅∇)−1∇⋅∇Φ=−∇Φ, (66)

assuming that satisfies the same boundary conditions as the operator . Hence 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 . In (65), we see that this operator is applied to which is of the form , where is a scalar field. In order for to be a curl-free vector field, we need , i.e. we need and to be parallel vectors. In the discretized case, our grid satisfies an orthogonality property as mentioned in Section 2.2, and the gradients of and 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 by the identity and obtain the following Schur complement approximation

 ~ST=cv∇⋅(ρu)−∇⋅(kT∇). (67)

Similar heuristic arguments for replacing by the identity operator in the case of the Stokes problem can be found, for example, in maday1993analysis (), and for the Navier-Stokes equations, in kay2002preconditioner ().

#### 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 , while injection wells satisfy . Thus,

 ST= ∇⋅(cvρu)−∇⋅(kT∇)−cvfprod+cv∇T⋅(ρu)T +cvT∇⋅(ρu)T−cvT(fprod)T−cvTinj(finj)T −(cv∇T⋅(ρu)p+cvT∇⋅(ρu)p−cvT(fprod)p−cvTinj(finj)p) (∇⋅(ρu)p−(fprod)p−(finj)