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
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.
Contents
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 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
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^{-1}b by p(A)b . The commonly used methods for constructing p(A)b are:
- Orthomin^{[1]}
- GMRES^{[2]}
Both methods minimize the residual norm over all vectors in span{b, Ab, A2b, …, A^{m-1}b} 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
- Compute r_{0} = b- Ax_{0}. Set p_{0} = M^{-1}r_{0}.
- For j = 0,1,2, …, until convergence Do :
- x_{j+1} = x_{j} + α_{j}p_{j}
- r_{j+1} = r_{j} - α_{j}Ap_{j}
- 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]}
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 |
References
- ↑ 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. http://dx.doi.org/10.2118/5729-MS
- ↑ 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. http://dx.doi.org/10.1137/0907058
- ↑ Appleyard, J.R. 1983. Nested Factorization. Presented at the SPE Reservoir Simulation Symposium, San Francisco, California, USA, 15-18 November. SPE-12264-MS. http://dx.doi.org/10.2118/12264-MS
- ↑ 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. http://dx.doi.org/10.2118/10493-MS
- ↑ 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. http://dx.doi.org/10.2118/2802-PA
- ↑ ^{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. http://dx.doi.org/10.2118/12265-MS
- ↑ 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. http://dx.doi.org/10.2118/9303-PA
- ↑ 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. http://dx.doi.org/10.2118/12259-MS
- ↑ 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.
- ↑ 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.
- ↑ 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.
- ↑ 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. http://dx.doi.org/10.2118/21226-MS
- ↑ 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. http://dx.doi.org/10.1137/0613057
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
Gridding in reservoir simulation
Upscaling of grid properties in reservoir simulation
High performance computing and reservoir simulation