Reservoir simulation linear equation solver

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

....................(1)

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

....................(2)

The well variables are then obtained by back substitution as

....................(3)

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:

• Orthomin
• GMRES

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

Preconditioning

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 r₀ = b- Ax₀. Set p₀ = M^-1r₀.
2. For j = 0,1,2, …, until convergence Do :
3. 
4. x_j+1 = x_j + α_jp_j
5. r_j+1 = r_j - α_jAp_j
6. 
7. 
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). 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. 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, 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 or CPR 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. 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.

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

• Reverse Cuthill-McKee (RCMK)
• Minimum Degree Fill (MDF).

Nomenclature

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

Subscripts

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

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

PEH:Reservoir_Simulation