Algebraic Multigrid Preconditioners for Multiphase Flow in Porous Media
Abstract
Multiphase flow is a critical process in a wide range of applications, including carbon sequestration, contaminant remediation, and groundwater management. Typically, this process is modeled by a nonlinear system of partial differential equations derived by considering the mass conservation of each phase (e.g., oil, water), along with constitutive laws for the relationship of phase velocity to phase pressure. In this study, we develop and study efficient solution algorithms for solving the algebraic systems of equations derived from a fully coupled and timeimplicit treatment of models of multiphase flow. We explore the performance of several preconditioners based on algebraic multigrid (AMG) for solving the linearized problem, including “blackbox” AMG applied directly to the system, a new version of constrained pressure residual multigrid (CPRAMG) preconditioning, and a new preconditioner derived using an approximate Schur complement arising from the block factorization of the Jacobian. We show that the new methods are the most robust with respect to problem character as determined by varying effects of capillary pressures, and we show that the block factorization preconditioner is both efficient and scales optimally with problem size.
1 Introduction
Multiphase flow is a feature of many physical systems and models of it are used in applications such as reservoir simulation, carbon sequestration, ground water management and contaminant transport. Modeling multiphase flow in highly heterogeneous media with complex geometries is difficult, especially when realistic processes such as capillary pressure are included. The system describing multiphase flow consists of nonlinear partial differential equations, constitutive laws and constraints. In this paper, we focus on the iterative solution of linear systems arising in a fully implicit cellcentered finite volume discretization of single component isothermal twophase flow model with capillary pressure. This fully implicit timestepping scheme is among the most robust for simulation of subsurface flow. Moreover, it can serve as a basis for modeling more complex processes in which the physical quantities are tightly coupled. This additional complexity could include adding more components, miscibility between components, thermal effects, and phase transitions.
The fully implicit discretization gives rise to a nonlinear system of equations at each time step. We employ a variant of Newton’s method with an exact Jacobian of the discretized equations to solve this system. For the linear system, we use a preconditioned generalized minimal residual (GMRES) method. There is a vast literature on different approaches to precondition the Jacobian system. A very popular approach is to use incomplete LU factorization (ILU) for constructing the preconditioner. Though popular for its general applicability, ILUbased preconditioners are neither effective nor scalable in many cases. Another approach is to consider decoupled preconditioners for the coupled system
[6].
This methodology is based on a direct solution of the decoupled pressure system, followed by an iterative solution using ILU for the global system.
This formulation was refined in [31], where it was proposed to solve
the pressure system iteratively, giving rise to the decoupled implicit pressure explicit saturation (IMPES) preconditioner. The effect of the decoupling is to weaken the coupling between pressure and saturation. Thus, it is often used as a preprocessing step to produce a modified Jacobian system, for which new preconditioners can be developed [10, 30]. With recent development of algebraic multigrid (AMG) algorithms, the pressure block can be solved efficiently using AMG, resulting in the constrained pressure residual multigrid (CPRAMG) approach. In recent developments, AMG has also been applied to solve the coupled system with some success [10, 30], although developing a general AMG algorithm for these types of problems remains a topic of ongoing research [21]. Since the Jacobian matrix has a block structure, one can also consider a block LU decomposition with an approximate Schur complement, which has been successfully applied to other models of fluid dynamics [22, 33]. Besides AMGbased methods, geometric multigrid has also been applied successfully to solve these types of problems [4, 5]. Our focus in this study is on methods based on AMG because of its general applicability.
In this paper, we develop a new block preconditioner designed to respect the coupling inherent in models of multiphase flow, and we report our experience with the performance and scalability of four different preconditioning strategies: (1) a direct algebraic multigrid (AMG) preconditioner for the global system, (2) twostage CPRAMG method with correction for the pressure block, also known as the combinative twostage approach, (3) CPRAMG with corrections for both the pressure and saturation blocks, known as the twostage additive approach, and (4) the block factorization (BF) preconditioner. An outline of the paper is as follows. In section 2, we present the mathematical formulation for twophase flow in porous media and discretization schemes. In section 3, we describe the solution algorithms for the linearized system.
Numerical results for the algorithms are presented in section 4.
We conclude with some remarks and discussion of future research directions in section 5.
2 Problem Statement
We consider isothermal, immiscible twophase flow through a porous medium. For example, often in reservoir simulation, one phase is oil (the nonwetting phase) and the other is pure water (wetting phase); alternatively, in groundwater management, one may consider a system of contaminated water that infiltrates a domain saturated with air.
Conservation of mass of each of the phases leads to the following coupled PDEs:
(1)  
(2) 
in which are the saturation, are the densities, are the source terms of the wetting and nonwetting phases respectively, and is the porosity of the medium. We assume a common extension of Darcy’s law to multiphase flow and express the phase velocities as
(3) 
Here, is the absolute permeability tensor. The terms are the relative permeability, viscosity, and pressure of phase respectively, is the gravitational constant, and is the depth. We also define the phase mobility . To close the system, we also have the following constitutive law and constraint
(4)  
(5) 
From equations LABEL:mass_n and LABEL:mass_w, one can derive separate equations for pressure and saturation. The pressure equation is elliptic or parabolic (diffusion); the saturation equation is hyperbolic or convectiondominated. The pressure equation is solved implicitly, and depending on the time discretization strategies applied to the saturation equation, several methods have been developed. In the case where the saturation equation is discretized using an explicit method (e.g., forward Euler), it is referred as IMPES (implicit pressure explicit saturation)[3]; for an implicit time discretization of the saturation equation, the method is known as the sequential approach, which was first applied to the blackoil model by Watts in 1985 [32].
The appeal of these methods lies in the sequential decoupling between pressure and saturation variables. Each equation can be solved separately. In addition, knowing the characteristics of each equation facilitates the design of efficient preconditioners, which is critical to achieving high performance. Both of these methods have been successfully applied to many problems where the fully implicit method is difficult to implement or shown to be too costly. However, the solution obtained from these approaches may lose accuracy if pressure and saturation are strongly dependent, or if capillary pressure changes very quickly. The lack of accuracy of these methods can be even more pronounced if more complex processes such as miscibility, thermal, and phase transitions are included in the model. For a more complete summary of the advantages and disadvantages of these approaches, we refer to [19].
Substitution of LABEL:capillary_pressure and LABEL:darcys_law into LABEL:mass_n and LABEL:mass_w and using the constraint LABEL:constraint yields a system of two equations and two unknowns. Using one popular choice of primary variables, the pressure in the wetting phase and saturation in the nonwetting phase, , we obtain
(6)  
(7) 
In this paper, we consider solving the coupled system consisting of LABEL:balance2 and LABEL:balance1 fully implicitly. We use a cellcentered finite volume method for spatial discretization and the backward Euler method for time discretization, similar to an approach defined in [13]. This will serve as a base model for adding more complexity in the future. The finite volume method described below is known for its mass conservation property. In addition, it can deal with the case of discontinuous permeability coefficients, and it is relatively straightforward to implement. Under appropriate assumptions, this method also falls into the mixed finite element framework [23, 26]. For simplicity, we consider a uniform partitioning of the domain into equal sized cells , i.e., . Let denote the area of the face between cells and . For each cell , integration of the mass conservation equations and the divergence theorem gives
(8) 
where the storage and the flux terms are approximated using the midpoint rule which is secondorder accurate:
(9) 
The surface integrals are discretized using twopoint flux approximation (TPFA); dropping the phase subscript, this gives
(10)  
(11) 
The index signifies an appropriate averaging of properties at the interface between cell and . The coefficients are approximated by upwinding based on the direction of the velocity field, i.e.,
(12) 
and the absolute permeability tensor on the faces are computed using harmonic averaging,
(13) 
Discretization in time using the backward Euler method gives a fully discrete system of nonlinear equations,
(14) 
3 Solution Algorithms
The system of nonlinear equations (LABEL:fully_discrete_system) can be written generically as where . We solve the system using Newton’s method, which requires solution of a linear system at each iteration :
(15) 
In our case, the solution vector consists of all the pressure and saturation unknowns at all the cell centers. The Jacobian system resulting from the derivative is often very difficult to solve using iterative methods, and preconditioning is critical for rapid convergence of Krylov subspace methods such as GMRES [28]. Next, we discuss the linear system arising from the Newton’s method and give a detailed description of the solution algorithms we will use to solve this system.
3.1 Linear System
For the set of primary variables , assuming unknowns corresponding to physical variables are grouped together and unknowns associated with nodal points in the domain are ordered lexicographically, each nonlinear Newton iteration entails the solution of a discrete version of a block linear system of the form
(16) 
in which
(17)  
(18) 
All the coefficients in equation (LABEL:linearize_fc) are evaluated at the linearization point . In a more concise form, the Jacobian matrix of the system has block structure
(19) 
and the linear system is . The characteristics of the matrix have been discussed in numerous papers [5, 13, 17, 30]. We summarize important characteristics of the operators here:

is nonsymmetric and indefinite

The block has the structure of a discrete purely elliptic problem for pressure.

The coupling block has the structure of a discrete firstorder hyperbolic problem in the nonwetting phase saturation.

The coupling block has the structure of a discrete convectionfree parabolic problem in the wetting phase pressure.

The block has the structure of a discrete parabolic (convectiondiffusion) problem for saturation when capillary pressure is a nonconstant function of the saturation. When capillary pressure is zero or a constant, there is no diffusion term and the block has the form of a hyperbolic problem.

Under mild conditions, i.e. modest timestep size, the blocks are diagonally dominant.
In this paper, we present some numerical results that show how different models of capillary pressure affect the algebraic properties of the (2,2)block in particular and the global system in general, which consequently determines the success of AMG solution algorithms. Our emphasis is on the development and use of preconditioning operators denoted , for the purpose of solving preconditioned systems
(20) 
3.2 Algebraic Multigrid
Multigrid is a highly efficient and scalable method available for solving large sparse linear systems [34]. Geometric multigrid uses a hierarchy of nested grids, whose construction depends on the geometry of the problem and a priori knowledge of the grids. AMG methods such as those developed in [29] have the advantage of not requiring an explicit hierarchy of nested grids. AMG constructs coarse grids based on the matrix values only, which makes it suitable for solving a wide range of problems on complicated domains and unstructured grids. Despite its successful application to scalar problems, using AMG for coupled systems is still relatively limited. Some attempts to use AMG to solve fully coupled systems encountered in modeling multiphase flow for reservoir simulation include [10, 30]. In this work, we use BoomerAMG [16], part of the Hypre package [14, 15], as a blackbox AMG solver. We note that in order to use BoomerAMG for the coupled system in our case, the Jacobian matrix needs to be ordered by grid points, i.e.
(21) 
in which is the number of grid points, and are matrices representing the couplings between pressure and saturation at points and . This is called the “point” method in [30].
3.3 Twostage Preconditioning with AMG
Unlike AMG, which has not been popular in reservoir simulation until recently, twostage preconditioners are widely used [17]. This idea was developed and first appeared in the context of multiphase flow modeling in the work of Wallis [31]. Following [13], we refer to this method as the constrained pressure residual (CPR) approach. There are many variants of twostage preconditioners. We discuss two algorithms here: the twostage combinative preconditioner  CPRAMG(1), and the twostage additive preconditioner  CPRAMG(2) [2].
Algorithm 1. Twostage Combinative  CPRAMG(1)

At each iteration let the residual be .

Solve .

Update the residual .

Solve for the pressure correction

Update the solution
Algorithm 2. Twostage Additive  CPRAMG(2)

At each iteration let the residual be .

Solve .

Update the residual .

Solve for the pressure correction

Solve for the saturation correction

Update the solution
The matrices denote the restriction of the global unknown vector to those associated with pressure and saturation respectively. That is, and for
(22) 
Then, in matrix form, the action of the twostage preconditioners can be expressed as
(23)  
(24) 
For the preconditioner in step 2 of both algorithms, we use the incomplete factorization with no fill ILU(0) method. For the correction solve, we apply AMG with one Vcycle iteration. The combinative approach with AMG was presented in [18, 20]. However, this method does not work well in the presence of fast changing capillary pressure. We confirm this observation in the next section. To deal with fast changing capillary pressure, we employ an additive CPRAMG approach, which involves one more AMG solve for the correction of the saturation block. The intuition is that when the absolute value of the derivative of capillary pressure is large, the block becomes diffusion dominated, and AMG can handle it efficiently.
3.4 Block Factorization Preconditioners
Consider the following decomposition of the Jacobian,
where is the Schur complement
(25) 
We could choose
(26) 
as an uppertriangular block preconditioner; this incorporates the effects of the coupling block , which contains the time derivative and gravity terms (LABEL:linearize_fc). We use an approximation of the Schur complement in which is replaced by its diagonal values:
(27) 
The purpose of this is to keep the Schur complement sparse so that the action of its inverse can be applied efficiently. This idea is the basis of the SIMPLE method used in other models of fluid dynamics [22]. A similar approach has also been applied to problems in single phase flow coupled with geomechanics in [33].
Algorithm 3. Block factorization preconditioner

At each iteration let the residual be .

Solve for the saturation using AMG.

Compute the residual for pressure .

Solve for the pressure using AMG.
An important advantage of this algorithm is that it does not rely on an ILU factorization. In matrix form,
(28) 
4 Numerical Results
In this section, we perform numerical experiments for the four aforementioned preconditioners. All of them are implemented in Amanzi, a parallel opensource multiphysics C++ code developed as a part of the ASCEM project [1]. Although Amanzi was first designed for simulation of subsurface flow and reactive transport, its modular framework and concept of process kernels [11] allow new physics to be added relatively easily for other applications. The twophase flow simulator employed in this work is one such example. Amanzi works on a variety of platforms, from laptops to supercomputers. It also leverages several popular packages for mesh infrastructure and solvers through a unified input file. Here, all of our experiments use a classical AMG solver through BoomerAMG in Hypre. The ILU(0) method is from Euclid, also a part of Hypre. GMRES is provided within Amanzi. For simplicity, we employ structured Cartesian grids for the test cases, but we can also use unstructured Korthogonal grids. This section has three parts. In the first part, we show the results for a twodimensional oilwater model problem. Although the problem is small, it is difficult to solve due to the heterogeneity of the permeability field. In the second part, we report the results for a threedimensional example. In the last part, we examine the scalability of the three preconditioning strategies. Unless specified otherwise, we use the benchmark problem SPE10 [8] for permeability data and porosity.
4.1 Twodimensional oilwater problem
The domain is a rectangle with dimensions meters. The mesh is , which means that the problem is truly twodimensional in the plane. We inject pure water into the domain through the boundary at the lower left corner, and oil and water exit the domain through the top right corner. These correspond to the at the inlet, and at the outlet. The simulation is run for 200 days with time step .
For capillary pressure models, we employ a simple linear model and the BrooksCorey [7] model:
(29) 
in which is the effective saturation, is the entry pressure, and is related to the poresize distribution. For the BrooksCorey model, the typical range of is [4, 12]. In general, is greater than 2 for narrow distributions of pore sizes, and is less than 2 for wide distributions. For example, sandpacks with broader distributions of particle sizes have ranging from to [24]. The BrooksCorey capillary pressure curves for various values of are plotted in Figure LABEL:brooks_corey. Other parameters are listed in Table LABEL:params_q5 and example 1 of Table LABEL:pc_params.
For all of the simulations presented here, the convergence tolerance for Newton’s method is , and the linear tolerance for GMRES is , which is the default in Amanzi. BoomerAMG is used as a preconditioner. The number of Vcycle steps is set to 1. The coarsening strategy is the parallel ClearyLubyJonesPlassman (CLJP) coarsening [9]. The interpolation method is the classical interpolation defined in [25], and the smoother is the forward hybrid GaussSeidel / SOR scheme.
Initial wetting phase pressure  

Initial nonwetting phase saturation  0.8 
Residual wetting phase saturation  0.0 
Nonwetting phase density  700 
Wetting phase density  1000 
Nonwetting phase viscosity  10.0 
Wetting phase viscosity  1 
Porosity  0.2 
Parameters  Ex 1  Ex 2  Ex 3  Ex 4 
Linear entry pressure  
BrooksCorey entry pressure  
BrooksCorey  2.5  0.8  2.5  0.8 
In order to explore the effects of different models for capillary pressure on solver performance, we use the four sets of parameters listed in Table LABEL:pc_params. In Example 1, the parameters are chosen such that the norm of the derivative of capillary pressure is large, leading to a diffusiondominated case (see equation (LABEL:linearize_fc)). In Example 2, the parameters are tuned to reduce the norm of , leading to an advectiondominated case. Example 3 is a more extreme case of example 2, in which is further decreased, leading to a strongly advectiondominated case. We also note the difference between the linear model and the BrooksCorey model for capillary pressure. The derivative for the linear model is a constant value, which means that the character of the problem, i.e. diffusiondominated or advectiondominated, is the same everywhere for the whole domain. In the BrooksCorey model, depends on the saturation of the wetting phase, and the problem can be diffusiondominated in one part of the domain, and advectiondominated in another part. This can cause further difficulties for AMGbased solvers, whose optimal performance is sensitive to the characteristics of the problem.
Methods/Models  Linear  Brooks Corey  

NI  LI  LI/NI  Time  NI  LI  LI/NI  Time  
AMG  32  368  11.5  27.2  36  470  13.1  37.24 
CPRAMG(1)  32  3695  115.5  324.15  36  5831  162  567.7 
CPRAMG(2)  32  899  28.1  103.94  36  1102  30.6  134.6 
BF  32  524  16.4  33.17  36  599  16.6  46.2 
The performance of the three strategies is summarized in Tables LABEL:performance_q5_1, LABEL:performance_q5_2, and LABEL:performance_q5_3. NI denotes the number of nonlinear iterations, LI the number of linear iterations, LI/NI the average number of linear iterations per nonlinear iterations, and Time the total time in seconds of the whole simulation. For the diffusiondominated problem for which the results are shown in Table LABEL:performance_q5_1, AMG is the most efficient method, about 25% more efficient than the block preconditioner in terms of both iteration counts (linear iterations per Newton step) and total run time. Note that for this example, the diffusion terms in the (2,2)block () are large and the performance of AMG tends to be very strong for problems of this type. For the linear model, the block factorization approach still takes about 8 times fewer linear iterations, and it is about 10 times faster in total run time than CPRAMG(1). The reason for this discrepancy is that CPRAMG(1) is a twostage preconditioner, and it requires an extra global solve using ILU. The block factorization preconditioner does not rely on ILU, which helps improve the run time significantly. CPRAMG(2) also performs well in this case. Although it requires one more AMG solve per Newton iteration than CPRAMG(1), it still outperforms CPRAMG(1) in terms of both the number of linear iterations per Newton step and the total run time. The same conclusion can be made for the BrooksCorey model.
Methods/Models  Linear  Brooks Corey  

NI  LI  LI/NI  Time  NI  LI  LI/NI  Time  
AMG  37  2575  69.6  138.8         
CPRAMG(1)  37  1919  51.9  175.5  55  4851  88.2  605.7 
CPRAMG(2)  37  1222  33.0  157.1  55  3701  67.3  506.8 
BF  37  684  18.5  51.7  55  1633  29.7  131.1 
Methods/Models  Linear  Brooks Corey  
NI  LI  LI/NI  Time  NI  LI  LI/NI  Time  
AMG                 
CPRAMG(1)  43  1079  25.1  122.8  48  2173  45.3  247.6 
CPRAMG(2)  43  1442  35.5  169.8  48  4805  100.1  560.5 
BF  43  1002  23.3  69.8  48  1829  38.1  121.8 
The results reported in Table LABEL:performance_q5_2 reveal the lack of robustness of AMG when applied to the coupled system. In contrast to the diffusiondominated case, for the linear model of capillary pressure, AMG requires the highest number of linear iterations per Newton step for the advectiondominated case, and it even diverges for the BrooksCorey model. The block factorization preconditioner still shows good performance, taking about half the number of iterations and running four times faster than the next best method, which is CPRAMG(2). CPRAMG(1) is still the least effective method in this case for both capillary pressure models.
For the strongly advectiondominated problem with parameters in example 3, AMG diverges for both the linear and BrooksCorey capillary pressure models. The performance of CPRAMG(2) is also affected in this case, trailing that of CPRAMG(1). CPRAMG(2) is still more robust than direct application of AMG, however, since unlike AMG, this method still converges. The block factorization preconditioner is again the most effective method, requiring fewer number of iterations and about half the run time of CPRAMG(1). This suggests that when the diffusion term in the block gets small, the coupling block which has the structure of a discrete firstorder hyperbolic problem for the saturation, becomes important and needs to be taken into account.
4.2 Twodimensional problem with gravity
In this example, we compare the performance of the different strategies for a problem in which gravity plays a dominant role. The domain is a square box of size meters. The absolute permeability is a homogeneous field of millidarcy. Water is injected into the domain through the boundary at the top left corner, and the outlet is at the top right corner. The rate of injection is 5 /day. For spatial discretization, we use uniform grids of size , , , and respectively. The initial conditions are the same as the heterogeneous twodimensional example above. The time steps are 10, 4, and 1 days, and the final times are 20, 8, and 2 days respectively.
Methods/Mesh sizes  

AMG  7  7  7  7 
CPRAMG(1)  15.1  25.9  49.4  95.2 
CPRAMG(2)  22.0  30.9  38.1  40.7 
BF  19.9  21.0  21.1  21.1 
Methods/Mesh sizes  

AMG        17.3 
CPRAMG(1)  17.9  17.8  16.1  25.2 
CPRAMG(2)  30.3  29.8  22.7  30.2 
BF  13.0  17.0  21.1  23.9 
Methods/Mesh sizes  

AMG         
CPRAMG(1)  18.6  19.6  19.7  18.8 
CPRAMG(2)  31.5  34.6  36.6  36.3 
BF  13.1  9.3  12.0  16.4 
The diffusiondominated case, shown in Table LABEL:2d_grav_diff, exhibits the same pattern as in the previous example: the AMG preconditioner is the most efficient method, followed by the block factorization method, CPRAMG(2), and CPRAMG(1). AMG and the block factorization method exhibit optimal performance with respect to problem size. The number of iterations for CPRAMG(2) also seems to reach a plateau as the mesh size is refined. In contrast, the performance of CPRAMG(1) does not scale well with respect to mesh size for this case, taking about twice the number of iterations for each level of mesh refinement.
The results for the advectiondominated case are shown in Table LABEL:2d_grav_adv. The AMG method is not robust and only converges for the largest mesh size (for which it takes the fewest iterations). The block factorization preconditioner is highly robust and also appears to require iteration counts tending to a constant as the mesh is refined. The performance of CPRAMG(2) is consistent except for the mesh. Although it requires more iterations than CPRAMG(1), this method shows promising scaling property, similar to the previous example, since the number of iterations does not grow as the mesh is refined. CPRAMG(1) performs quite well for this case, but it still exhibits poor scalability as the number of iterations grow quickly between and .
In the strongly advectiondominated case, AMG diverges for all mesh sizes. The new block factorization is the most efficient method in this case, requiring the smallest number of iterations across all mesh sizes. Here, CPRAMG(1) is more efficient than CPRAMG(2), requiring about half the number of iterations. Both CPRAMG(1) and CPRAMG(2) show good scaling property in this case. The scaling result for the block factorization method is not as clear as in the diffusion and advection dominated cases, but we suspect that the mesh is not fine enough for a consistent pattern to emerge.
Besides varying the mesh size, we also experimented with changing the time step size for a fixed mesh of for the same problem. The results are reported in Table LABEL:varying_time_steps. The final time for the simulation is 8 days. The results are reported in Table LABEL:varying_time_steps.
Methods/Time steps  = 1 day  = 2 days  = 4 days  = 8 days  

NI/TS  LI/NI  NI/TS  LI/NI  NI/TS  LI/NI  NI/TS  LI/NI  
CPRAMG(1)  12.4  16.9  17  15.9  23  16.1  28  19.4 
CPRAMG(2)  12.4  29.0  17  24.0  23  22.7  28  27.2 
BF  12.4  16.7  17  19.0  23  21.1  28  23.0 
NI/TS is the number of Newton iteration per time step.
Since AMG does not converge in this experiment, we exclude it from the results. From Table LABEL:varying_time_steps, it is clear that as the time step gets larger, Newton’s method takes more iterations to converge. For 8 days, there is only one time step and it is the most difficult case. CPRAMG(1) number of iterations is not significantly affected by the time step except for the largest time step size of 8 days. Meanwhile, CPRAMG(2) number of iterations decreases as the time step gets larger, but goes up again at days. The block factorization method shows consistent increase in the number of iterations for larger time steps. Overall, there is not much of a difference in terms of iteration counts for these three methods, but it is worth noting that the block factorization method is much faster than the others in terms of run time, as it does not require a global ILU solve.
4.3 Behavior of Eigenvalues
It is often possible to obtain insight into the properties of preconditioning operators from the eigenvalues of the preconditioned matrix . In particular, recall a standard analysis of the convergence behavior of GMRES for solving the preconditioned system (LABEL:preconsystem) [27]. Assume the preconditioned matrix is diagonalizable, , where is a diagonal matrix containing the eigenvalues of the preconditioned matrix and the columns of are the corresponding eigenvectors. If are the iterates obtained at the th step of GMRES iteration, with residual , then
(30) 
where the minimum in (LABEL:GMRESbound) is over all polynomials of degree at most that have the value at the origin, is the set of eigenvalues of , and the norm is the vector Euclidian norm. Thus, a good preconditioner tends to produce a preconditioned operator with a compressed spectrum whose entries are not near the origin. In this section, we explore the behavior of the eigenvalues of the preconditioned matrix with an eye toward understanding the effects of features of the discrete problem such as discretization mesh size and qualitative features of the model such as the relative weights of diffusion and advection and the degree of coupling between the components.
Figure LABEL:evalues_plot gives a representative depiction of the eigenvalues of preconditioned operators for three of the preconditioners considered. These results are for benchmark problems for which performance is considered in section LABEL:2Doilwater, the twodimensional linear oilwater model discretized on a grid. The plots on the left side of the figure show eigenvalues for the diffusiondominated case (for which solution performance is shown in Table LABEL:performance_q5_1), and those on the right show eigenvalues for the advectiondominated case (performance in Table LABEL:performance_q5_2). ^{1}^{1}1These computations were done using the eig function in Matlab, and they use Matlab backslash to perform the actions of the inverses of , and the modified Schur complement. This contrasts with the solution algorithms tested, which approximate these operations using one AMG Vcycle.
These displays indicate that the spectra for the preconditioned systems for the CPRAMG(2) and BF preconditioners are bounded away from the origin, whereas for the CPRAMG(1) preconditioner there are many small eigenvalues. Performance of CPRAMG(1) improves in the advectiondominated case, and the smallest associated eigenvalues are somewhat further from the origin. In contrast, the latter two preconditioners are largely unchanged in the advectiondominated case, where they are still effective, and the associated eigenvalues are also contained in similarly structured regions far from the origin. We believe the superior performance of the BF preconditioner comes from its greater emphasis on the coupling between pressure and saturation, derived from use of the approximate Schur complement (LABEL:SIMPLESchur).
4.4 Threedimensional Problem
We use a homogeneous permeability field of millidarcy, and the grid is stretched to induce anisotropy. The model dimensions are meters and the cell size is meter. Thus, the mesh is , and the problem has 1.2 million unknowns in total. Water is injected into the domain at one bottom corner and the outlet is at the opposite corner. The injection rate is 0.75 . The parameters for the capillary pressure model is from example 1 of Table LABEL:pc_params. The simulation is run for 100 days with time step days.
Methods/Models  Linear  Brooks Corey  

NI  LI  LI/NI  Time  NI  LI  LI/NI  Time  
AMG  16  282  17.6  103.1  20  452  22.6  144.7 
CPRAMG(1)  16  2698  168.6  803.2  20  6069  303.45  1940.8 
CPRAMG(2)  16  712  44.5  299.5  20  1900  95.0  741.1 
BF  16  355  22.2  133.6  20  752  37.6  231.1 
Table LABEL:performance_3d_1 shows the performance results of the diffusiondominated case for this 3D example, which are consistent with those of the previous twodimensional example. AMG preconditioner shows the best results for both the iteration counts per Newton step and the time it takes to complete the simulation for both capillary pressure models. CPRAMG(2) does not perform quite as well as AMG, but it is much more efficient than CPRAMG(1) for both performance measures and capillary pressure models. As in the twodimensional case, the new block factorization method performs well, requiring fewer than half the iterations than CPRAMG(2) for both the linear and BrooksCorey models, and running in about one third the CPU time.
We also tested the threedimensional SPE10 problem with the linear model of capillary pressure for the different preconditioning strategies. Here, AMG diverges even for the diffusiondominated case, even though it was the most efficient method for the twodimensional example. The block factorization method is about four times faster than CPRAMG(2) and five times faster than CPRAMG(1) in the diffusion dominated case. CPRAMG(2) still outperforms CPRAMG(1) both in terms of iteration counts and run time, but the margin is smaller than for the twodimensional problem. In the advectiondominated case, unlike in the twodimensional example, CPRAMG(1) is more efficient than CPRAMG(2), requiring about 23% fewer number iterations and 45% run time. The block factorization approach is still the most efficient method, taking fewer than half the number of iterations and less than half the run time of CPRAMG(1). We also note that the number of iterations for the block factorization method is very consistent with respect to the characteristics of the problem, i.e. it does not change significantly whether the problem is diffusiondominated, advectiondominated, or strongly advectiondominated.
Methods/Models  Linear  

NI  LI  LI/NI  Time (s)  
AMG         
CPRAMG(1)  17  2410  141.8  614.14 
CPRAMG(2)  17  1661  97.7  448.21 
BF  17  490  28.8  121.71 
Methods/Models  Linear  

NI  LI  LI/NI  Time (s)  
AMG         
CPRAMG(1)  18  1122  62.3  354.38 
CPRAMG(2)  18  1554  86.3  657.12 
BF  18  474  26.3  157.24 
4.5 Scaling Results
To perform a scalability study, we run a test problem on a box of dimensions meters. The initial mesh is and is repeatedly refined in the zdirection. The domain has constant material properties. The parameters for the water retention models are listed in example 4 of Table LABEL:pc_params. Note that this set of parameters corresponds to a diffusiondominated problem.
The results shown in Figure LABEL:weak_scaling. indicate that the performance of the block factorization, CPRAMG(2), and AMG methods is independent of the mesh size. The number of linear iterations per Newton step does not grow as the mesh is refined which is optimal multigrid performance. The block factorization method’s performance is nearly identical to that of AMG for the linear model, and still quite close for the BrooksCorey model, compared to CPRAMG(2). CPRAMG(1), however, does not scale as well as the other two methods. The linear iteration counts for CPRAMG(1) grows linearly as the mesh is refined.
5 Conclusions
In this work, we have implemented a fully implicit parallel isothermal twophase flow simulator along with four different preconditioning strategies to solve the linear systems resulting from linearization of the coupled equations, and we have tested the performance of these methods as preconditioners for GMRES. We have also developed a new block factorization preconditioner whose performance is robust and efficient across all benchmark problems studied. In contrast, although AMG preconditioning applied to the coupled systems is the most efficient choice in some cases (both twodimensional and threedimensional diffusiondominated examples), it exhibits slow convergence and sometimes diverges for advectiondominated cases. The new block factorization preconditioner achieves consistently low iteration counts across all the tests and varying examples of capillary pressure, and it scales optimally with problem size. The combinative CPRAMG(1), though robust across all the tests, is the least efficient method, with the exception of the near hyperbolic case where it is faster than CPRAMG(2). The additive CPRAMG(2) method performs well in most cases except the strongly advectiondominated case. It also scales optimally with problem size for both advectiondominated and diffusiondominated case.
Acknowledgement
This work was supported by the US Department of Energy under grant DESC0009301, the U.S. National Science Foundation under grant DMS1418754, and by the Department of Energy at Los Alamos National Laboratory under contract DEAC5206NA25396 through the DOE Advanced Simulation Capability for Environmental Management (ASCEM) program, and the LANL Summer Student Program at the Center for Nonlinear Studies (CNLS).
References
 [1] ASCEM. http://esd.lbl.gov/research/projects/ascem/thrusts/hpc/, 2009.
 [2] O. Axelsson. Iterative Solution Methods. Cambridge University Press, 1994.
 [3] K. Aziz and A. Settari. Petroleum Reservoir Simulation. Applied Science Publishers, 1979.
 [4] P. Bastian. Numerical computation of multiphase flow in porous media. Habilitationsschrift, 1999.
 [5] P. Bastian and R. Helmig. Efficient fullycoupled solution techniques for twophase flow in porous media: Parallel multigrid solution and large scale computations. Advances in Water Resources, 23(3):199 – 216, 1999.
 [6] A. Behie and P.K.W. Vinsome. Block iterative methods for fully implicit reservoir simulation. Society of Petroleum Engineers Journal, 22(05):658–668, oct 1982.
 [7] R. Brooks and A. Corey. Hydraulic Properties of Porous Media. Hydrology Papers, Colorado State Univeristy, 1964.
 [8] M.A. Christie and M.J. Blunt. Tenth SPE comparative solution project: A comparison of upscaling techniques. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers (SPE), 2001.
 [9] A.J. Cleary, R.D. Falgout, V.E. Henson, and J.E. Jones. Coarsegrid selection for parallel algebraic multigrid. In Alfonso Ferreira, JosÃ© Rolim, Horst Simon, and ShangHua Teng, editors, Solving Irregularly Structured Problems in Parallel, volume 1457 of Lecture Notes in Computer Science, pages 104–115. Springer Berlin Heidelberg, 1998.
 [10] T. Clees and L. Ganzer. An efficient algebraic multigrid solver strategy for adaptive implicit methods in oilreservoir simulation. SPE Journal, 15(03):670–681, 2010.
 [11] E. Coon, J. D. Moulton, and S. Painter. Managing complexity in simulations of land surface and nearsurface processes. Technical Report LAUR 1425386, Applied Mathematics and Plasma Physics Group, Los Alamos National Laboratory, 2014.
 [12] A. T. Corey. Mechanics of Immiscible Fluids in Porous Media. Water Resources Pubns, 1995.
 [13] C.N. Dawson, H. Klie, M.F. Wheeler, and C.S. Woodward. A parallel, implicit, cellcentered method for twophase flow with a preconditioned NewtonâKrylov solver. Computational Geosciences, 1(34):215–249, 1997.
 [14] R.D. Falgout, J.E. Jones, and U.M. Yang. The design and implementation of hypre, a library of parallel high performance preconditioners. In AreMagnus Bruaset and Aslak Tveito, editors, Numerical Solution of Partial Differential Equations on Parallel Computers, volume 51 of Lecture Notes in Computational Science and Engineering, pages 267–294. Springer Berlin Heidelberg, 2006.
 [15] R.D. Falgout and U.M. Yang. HYPRE: a library of high performance preconditioners. In Preconditioners,â Lecture Notes in Computer Science, pages 632–641, 2002.
 [16] V.E. Henson and U.M. Yang. BoomerAMG: a parallel algebraic multigrid solver and preconditioner. Applied Numerical Mathematics, 41:155–177, 2000.
 [17] H. Klie, M. Rame, and M.F. Wheeler. Twostage preconditioners for inexact Newton methods in multiphase reservoir simulation. Technical report, Center for Research on Parallel Computation, Rice University, 1996.
 [18] S. Lacroix, Yu. Vassilevski, J. Wheeler, and M.F. Wheeler. Iterative solution methods for modeling multiphase flow in porous media fully implicitly. SIAM J. Sci. Comput., 25(3):905–926, jan 2003.
 [19] B. Lu. Iteratively Coupled Reservoir Simulation for Multiphase Flow in Porous Media. PhD thesis, Center for Subsurface Modeling, ICES, The University of Texas  Austin, 2008.
 [20] R. Masson, R. Scheichl, P. Quandalle, C. Guichard, R. Eymard, R. Herbin, and K. Brenner et al. CPRAMG preconditioner for multiphase porous media flows: a few experiments. La page numérique du LATP  UMR 7353 Marseille, 2012.
 [21] D. OseiKuffuor, L. Wang, R. D. Falgout, and I. D. Mishev. An algebraic multigrid solver for fullyimplicit solution methods in reservoir simulation. SIAM Geosciences, Society for Industrial and Applied Mathematics (SIAM), 2015.
 [22] S.V Patankar and D.B Spalding. A calculation procedure for heat, mass and momentum transfer in threedimensional parabolic flows. International Journal of Heat and Mass Transfer, 15(10):1787 – 1806, 1972.
 [23] D.W. Peaceman. Fundamentals of Numerical Reservoir Simulation. Elsevier, Amsterdam, 1977.
 [24] PetroWiki. Relative permeability models. http://petrowiki.org/Relative_permeability_models.
 [25] J. W. Ruge and K. Stueben. 4. Algebraic Multigrid, chapter 4, pages 73–130. 1986.
 [26] T. F. Russell and M.F. Wheeler. Finite element and finite difference methods for continuous flows in porous media. In The Mathematics of Reservoir Simulation, chapter 2, pages 35–106. Society for Industrial and Applied Mathematics, 1983.
 [27] Y. Saad. Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics, second edition, 2003.
 [28] Y. Saad and M.H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput., 7(3):856–869, July 1986.
 [29] K. Stueben. A review of algebraic multigrid. Journal of Computational and Applied Mathematics, 128(1â2):281 – 309, 2001. Numerical Analysis 2000. Vol. VII: Partial Differential Equations.
 [30] K. Stueben, T. Clees, H. Klie, B. Lu, and M.F. Wheeler. Algebraic multigrid methods (AMG) for the efficient solution of fully implicit formulations in reservoir simulation. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers (SPE), 2007.
 [31] J.R. Wallis, R.P. Kendall, and L.E. Little. Constrained residual acceleration of conjugate residual methods. 1985.
 [32] J.W. Watts. A compositional formulation of the pressure and saturation equations. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers (SPE), SPE 12244, 1985.
 [33] Joshua A. White and Ronaldo I. Borja. Blockpreconditioned Newton–Krylov solvers for fully coupled flow and geomechanics. Computational Geosciences, 15(4):647–659, 2011.
 [34] U.M. Yang. Parallel algebraic multigrid methods â high performance preconditioners. In Numerical Solution of Partial Differential Equations on Parallel Computers, volume 51 of Lecture Notes in Computational Science and Engineering, pages 209–236. Springer Berlin Heidelberg, 2006.