You must log in to edit PetroWiki. Help with editing

Content of PetroWiki is intended for personal use only and to supplement, not replace, engineering judgment. SPE disclaims any and all liability for your use of such content. More information

Reservoir simulation linear equation solver

Jump to navigation Jump to search

The linear equation solver is an important component in a reservoir simulator. It is used in the Newton step to solve the discretized nonlinear partial differential equations. These equations describe mass balances on the individual components treated in the model.

Linear solver

For nonisothermal problems, an energy balance is added to the system. The matrix problem involves solving Ax=b, where A is typically a large sparse matrix, b is the right-side vector, and x is the vector of unknowns. In the Implicit Pressure Explicit Saturation (IMPES) formulation, there is a single unknown per cell pressure. In the fully implicit formulation, there is a fixed number n of unknowns per cell where n ≥ 2. In the adaptive implicit formulation, there is a variable number of unknowns per cell. In most formulations, pressure is an unknown for each cell. The matrix A typically has associated well constraint equations and well variables and may be partitioned in block 2 × 2 form as


where xw is the well variable-solution vector and xR is the reservoir variable-solution vector. The matrix Aww is often diagonal. In this case, the well variables may be directly eliminated, and the iterative solution is on the implicitly defined matrix system


The well variables are then obtained by back substitution as


If A is large, solution of the matrix equations is impractical using direct methods such as Gaussian elimination because of computer storage or CPU time requirements. Iterative solution based on projection onto Krylov subspaces is typically used. These Krylov subspaces are spaces spanned by vectors of the form p(A)v, where p is a polynomial. Basically, these techniques approximate A-1b by p(A)b . The commonly used methods for constructing p(A)b are:

Both methods minimize the residual norm over all vectors in span{b, Ab, A2b, …, Am-1b} at iteration m. They should yield identical results. From a practical standpoint, it does not matter which is used.


A technique known as preconditioning can improve both the efficiency (speed in a typical problem) and robustness (ability to solve a wide range of problems at least reasonably well) of Orthomin or GMRES. Preconditioning involves transforming the original matrix system into one with the same solution that is easier to solve. As a rule, the robustness of the iterative scheme is far more dependent on the preconditioning than on the specific Krylov subspace accelerator used. The preconditioner M is a matrix that approximates A, and has the property that the linear systems of the form Mx = b are easily and inexpensively solved. For most linear solvers the following preconditioned system is solved:


The preconditioned Orthomin(k) algorithm, which retains the last k A-orthogonal direction vectors, is given by

  1. Compute r0 = b- Ax0. Set p0 = M-1r0.
  2. For j = 0,1,2, …, until convergence Do :
  4. xj+1 = xj + αjpj
  5. rj+1 = rj - αjApj
  8. End Do

Because of the nature of the reservoir simulation equations, only certain preconditioners are effective in solving them. Reservoirs are typically shaped like pancakes, being much broader than they are thick. This geometry leads to strong vertical connectivity. Some preconditioners exploit this property. The most commonly used such preconditioner is nested factorization (NF).[3] The convergence of NF is sensitive to the cell ordering. Best results are usually obtained by ordering the cells first along the direction of highest transmissibility and then successively along directions of decreasing transmissibility. This nearly always means that NF should be ordered first in the vertical direction.

The other commonly used preconditioners are incomplete lower triangular/upper triangular factorizations of the matrix, or incomplete LU factorization (ILU(n)), where n is the level of infill that is retained during the elimination process. Performance of these can be improved by using a red-black checkerboard ordering (also called D4 ordering) of the nodes.[4] Red-black ordering on a five-point (in 2D) or seven-point (in 3D) grid leads to direct elimination of the unknowns at the red cells, leaving a system containing only the unknowns at the black cells. The result is a halving of the number of unknowns. An ILU preconditioner using red-black ordering with zero infill on the reduced system is referred to as RBILU(0) and is the most frequently-used form of ILU.

In IMPES models with either no faults or vertical faults only, RBILU(0) or ILU(1) combined with z-line additive corrections,[5][6] typically converges very rapidly. NF can also be used effectively in cases involving vertical faults and pinchouts because the matrix retains the structure required by NF. Nonvertical faults interrupt the matrix structure that makes red-black orderings attractive. In models containing them, ILU is the method of choice. Nine-point discretizations also cause problems for red-black orderings, but cause no difficulty for NF.

Combinative or CPR method

NF and RBILU(0) are commonly used in implicit models. Another very effective approach exploits the fact that pressure is the "stiff" variable. The Combinative[7] or CPR[6] method is a two-step preconditioning that extracts a pressure equation from the implicit matrix. It iteratively solves for a pressure correction at each iteration, uses the pressure correction to form a new residual, applies an inexpensive implicit preconditioning such as diagonal scaling or line Gauss-Seidel to the new residual, and then uses the sum of the two steps as the approximate solution. In compositional models, this two-step method can be much faster than one-step methods.

Parallel iteration

Many models include well-constraint equations that add well pressures to the set of unknowns. A simple but effective way of dealing with these equations is to order them first in the global matrix and then eliminate the well pressures from the set of unknowns.[8] In this approximate elimination, any infill terms are either column-summed or row-summed into the main diagonal of the reservoir matrix, which is then factored using NF. The reservoir matrix problem is then solved iteratively and the well variables are obtained by back substitution.

Parallel iterative solution typically uses a domain decomposition approach in which the grid is partitioned into domains that contain approximately the same number of cells. The partitioning should be done such that coupling is strongest between cells within a domain. As a result, domains normally are groups of columns of cells. One way to accelerate the iteration is to color the domains in red-black fashion and apply an NF-type procedure in which the outer level of nesting is the coupling between domains and either NF or an ILU variation is used to factor the individual domains.[9][10]

The solution of matrices arising from unstructured grids typically involves some variant of ILU with cell orderings such as:

  • Reverse Cuthill-McKee (RCMK)[11]
  • Minimum Degree Fill (MDF).[12][13]


A = a matrix
b = right-side vector
x = an unknown vector or concentration of a component in the liquid phase
α = mole fraction vapor phase


I = component number, index counter, or initial condition
j = index counter
J = phase number
k = index counter
r = residual
R = reservoir
o = oil phase
w = water phase or well


  1. Vinesome, P.K.W. 1976. Orthomin, an Iterative Method for Solving Sparse Banded Sets of Simultaneous Linear Equations. Presented at the SPE Symposium on Numerical Simulation of Reservoir Performance, Los Angeles, 19–20 February. SPE-5729-MS.
  2. Saad, Y. and Schultz, M.H. 1986. GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing 7 (3): 856-869.
  3. Appleyard, J.R. 1983. Nested Factorization. Presented at the SPE Reservoir Simulation Symposium, San Francisco, California, USA, 15-18 November. SPE-12264-MS.
  4. Tan, T.B.S. and Letkeman, J.P. 1982. Application of D4 Ordering And Minimization In An Effective Partial Matrix Inverse Iterative Method. Presented at the SPE Reservoir Simulation Symposium, New Orleans, Louisiana, 31 January-3 Feburary 1982. SPE-10493-MS.
  5. Watts, J.W. 1971. An Iterative Matrix Solution Method Suitable for Anisotropic Problems. Society of Petroleum Engineers Journal 11 (1): 47-51. SPE-2802-PA.
  6. 6.0 6.1 Wallis, J.R. 1983. Incomplete Gaussian Elimination as a Preconditioning for Generalized Conjugate Gradient Acceleration. Presented at the SPE Reservoir Simulation Symposium, San Francisco, California, USA, 15-18 November. SPE-12265-MS.
  7. Behie, A. and Vinsome, P.K.W. 1982. Block Iterative Methods for Fully Implicit Reservoir Simulation. Society of Petroleum Engineers Journal 22 (5): 658-668. SPE-9303-PA.
  8. Holmes, J.A. 1983. Enhancements to the Strongly Coupled, Fully Implicit Well Model: Wellbore Crossflow Modeling and Collective Well Control. Presented at the SPE Reservoir Simulation Symposium, San Francisco, California, 15-18 November 1983. SPE-12259-MS.
  9. Burrows, R., Ponting, D., and Wood, L. 1996. Parallel Reservoir Simulation with Nested Factorization. Presented at the 1996 European Conference on the Mathematics of Oil Recovery, Leoben, Austria, 3–6 September.
  10. Killough, J.E. et al. 1996. Parallelization of a General-Purpose Reservoir Simulator. Presented at the 1996 European Conference on the Mathematics of Oil Recovery, Leoben, Austria, 3–6 September.
  11. Cuthill, E. and McKee, J. 1969. Reducing the Bandwidth of Sparse Symmetric Matrices. Proc., 24th National Conference of the Assoc. for Computing Machinery, Brandon Press, Princeton, New Jersey, 157–172.
  12. D'Azevedo, E., Forsyth, P.A., and Tang, W.-P. 1991. An Automatic Ordering Method for Incomplete Factorization Iterative Solvers. Presented at the SPE Symposium on Reservoir Simulation, Anaheim, California, 17-20 February 1991. SPE-21226-MS.
  13. D’Azevedo, E.F., Forsyth, P.A., and Tang, W.-P. 1992. Ordering Methods for Preconditioned Conjugate Gradient Methods Applied to Unstructured Grid Problems. SIAM Journal of Matrix Analysis and Applications 13 (3): 944-961.

Noteworthy papers in OnePetro

Use this section to list papers in OnePetro that a reader who wants to learn more should definitely read

External links

Use this section to provide links to relevant material on websites other than PetroWiki and OnePetro

See also

Reservoir simulation

Gridding in reservoir simulation

Upscaling of grid properties in reservoir simulation

High performance computing and reservoir simulation

Reservoir simulation applications