Solving systems of linear equations in simulators

Quick refresher

Every implicit step is a system $\mathbf{A}\mathbf{x} = \mathbf{b}$ , with one equation per cell (or per phase of a cell). A row holds the diagonal and one number for each connection: seven numbers in a three-dimensional Cartesian grid, and the rest are zeros.
The matrix of the pressure equation is symmetric (the number in row $i$ and column $j$ equals the number in row $j$ and column $i$ ). It is diagonally dominant thanks to the storage term $S/\Delta t$ (the diagonal is larger than the sum of the other numbers in the row). A uniform shift of pressure is an eigenvector of this matrix: a pattern that the matrix only stretches, by the factor $S/\Delta t$ , which is its smallest eigenvalue. So the conditioning (how sensitive the solution is to errors) gets worse as the step grows. The Jacobian matrix of the fully implicit method is non-symmetric and is made of blocks.
Gaussian elimination removes the unknowns one by one, which leaves a triangular system (each equation has one unknown more than the one below it). Then it finds the unknowns by back substitution (from the last equation upwards). The elimination fills the band of the matrix. The work is of the order of $nb^2$ and the memory is $nb$ , where $n$ is the number of unknowns and $b$ is the band width. In a three-dimensional grid the width of the band is a whole layer of cells, so banded methods are unsuitable for large models.
Jacobi, Gauss–Seidel and over-relaxation improve each unknown from its own equation. They converge if the spectral radius of the iteration matrix (the largest factor by which one iteration can multiply an error pattern) is smaller than one. But on large grids they converge slowly: with $N$ cells along a side, Jacobi and Gauss–Seidel need $\propto N^2$ iterations, and over-relaxation needs $\propto N$ .
Krylov methods (conjugate gradients, GMRES) build the best approximation from the vectors $\mathbf{b}$ , $\mathbf{A}\mathbf{b}$ , $\mathbf{A}^2\mathbf{b}$ , … (each new vector costs only one multiplication by the matrix). For conjugate gradients the number of iterations grows as $\sqrt\kappa$ , where $\kappa$ is the condition number. A preconditioner, a rough but cheap copy of $\mathbf{A}$ , reduces the number of iterations several times over, especially in heterogeneous models.
Iterations are stopped by the relative residual. The residual of the pressure equation is the imbalance of a cell in m³/d: an incompletely solved system “creates” fluid, and its error goes straight into the material balance.

Key relationships

\left(\frac{S_i}{\Delta t} + \sum_j a_{ij}\right)\delta_i - \sum_j a_{ij}\,\delta_j = -q_i

a row of the pressure system of an implicit step (the pressures at the start of the step are equal)

W_{\text{band}} \approx n\,b^2,\qquad W_{\text{iter}} \approx \nu\,c\,n

work of banded elimination and of an iterative method

\nu \approx \frac{\ln(1/\varepsilon)}{-\ln\rho(\mathbf{G})}

iterations needed to reduce the error by a factor of $1/\varepsilon$

\frac{\lVert\mathbf{x} - \mathbf{x}^{(\nu)}\rVert}{\lVert\mathbf{x}\rVert} \le \kappa\,\frac{\lVert\mathbf{r}^{(\nu)}\rVert}{\lVert\mathbf{b}\rVert}

a small residual is not yet a small error if $\kappa$ is large

Common pitfalls

Inverting the matrix or the Jacobian: the inverse of a sparse system is almost full.
Believing that a small residual guarantees accurate pressures: the error is bounded by the residual multiplied by the condition number.
Loosening the solver tolerance for speed: an incompletely solved system creates a material balance error and extra Newton iterations.
Comparing iterative methods by their first iterations: the speed in the limit is set by the spectral radius.

After this lesson you can

describe the structure of the matrix of a simulation model: sparsity, the band, symmetry, diagonal dominance, and the dependence of the condition number on the time step
solve a small system by Gaussian elimination, explain the LU decomposition, fill-in and the role of the cell numbering, and estimate the work of a banded method
carry out iterations of the Jacobi, Gauss–Seidel and successive over-relaxation methods, and explain what their convergence depends on
explain the idea of Krylov subspace methods and of preconditioning (incomplete LU decomposition, the two-stage CPR, multigrid methods), and why heterogeneous models need them
choose a stopping criterion, link the residual to the material balance, and explain why the solver decides the run time

Where the systems of equations come from and how they are built

An implicit scheme (explicit and implicit schemes) turns every time step into a system of equations. The flows are calculated from the pressures at the end of the step, which are not yet known, and the equations of all cells are solved together. In the IMPES method this is one linear system for the pressures, with one equation per cell. In the fully implicit method the nonlinear equations are solved by Newton's method, and at every iteration we solve a linear system with the Jacobian matrix (a table of derivatives), with one equation per phase of a cell (residuals and the Jacobian matrix). A forecast run takes hundreds of steps with several iterations each, which makes thousands of systems with hundreds of thousands or millions of unknowns.

We write the pressure system of one implicit step for a single-phase, slightly compressible fluid. The unknowns are the changes of pressure over the step, $\delta_i = p_i^{n+1} - p_i^n$ :

A row of the pressure system: the cell equation at an implicit step

\left(\frac{S_i}{\Delta t} + \sum_j a_{ij}\right)\delta_i - \sum_j a_{ij}\,\delta_j = -q_i + \sum_j a_{ij}\left(p_j^{n} - p_i^{n}\right)

where:

$\delta_i,\ \delta_j$ changes of the pressure of the cell and its neighbours over the step, the unknowns, MPa
$S_i$ storage capacity of the cell $V_{p,i}c_t$ , m³/MPa
$a_{ij}$ link coefficient $0.0853\,T_{ij}/\mu$ , m³/(d·MPa)
$q_i$ withdrawal by a well at reservoir conditions (injection with a minus sign), m³/d
$p_i^{n},\ p_j^{n}$ pressures at the start of the step, MPa
$\Delta t$ time step, d

Units: metric oilfield units

Holds when

A linear pressure equation: the properties of the fluid and the rock are constant during the step (a nonlinear case is solved by iterations, and each one has such a system).
The sums run over all connections of the cell: neighbours across faces and non-neighbour connections. A well with a given bottomhole pressure adds its own coefficient to the diagonal.

This is the system of the lesson on the implicit scheme, written for the pressure changes. The link coefficients come from the transmissibilities (Interblock flow and transmissibility). In a row only the diagonal and the numbers of the neighbours are non-zero: the matrix is sparse.

Sparsity (sparse matrix). Nearly all the numbers of the matrix are zeros. A row has as many non-zero numbers as the cell has connections, plus the diagonal: three in a one-dimensional strip, five in a two-dimensional grid and seven in a three-dimensional Cartesian grid. Non-neighbour connections and wells add only a few.
The band. With natural numbering, first along the $x$ axis, then by rows and by layers, the neighbours of a cell are 1, $N_x$ and $N_xN_y$ places away from it in the numbering. The non-zero numbers lie on seven diagonals, which form a strip around the main diagonal, and this strip is the band. The width of the band is a whole layer of cells.
Symmetry. The link coefficient is the same for the pairs $(i, j)$ and $(j, i)$ , so the pressure matrix is symmetric. The Jacobian matrix of the fully implicit method is non-symmetric, because the mobility is taken from the upstream cell (what the Jacobian matrix is made of), and each cell corresponds to a block of 2 × 2 or 3 × 3 in it.
Diagonal dominance. The diagonal is larger than the sum of the absolute values of the other numbers in the row by $S_i/\Delta t$ (and by the coefficient of a well with a given bottomhole pressure). A symmetric matrix with such dominance and a positive diagonal is positive definite (all its eigenvalues are positive): the solution is unique, and many methods are certain to converge.

The sum of the numbers in each row equals the storage term $S_i/\Delta t$ : each link coefficient enters it twice, with a plus sign on the diagonal and with a minus sign off it. So, if the cells have the same capacity, a uniform change of pressure in the whole formation is an eigenvector of the matrix:

A uniform shift of pressure is an eigenvector of the matrix

\mathbf{A} = \frac{S}{\Delta t}\,\mathbf{I} + \mathbf{C},\qquad \mathbf{C}\,\mathbf{1} = 0,\qquad \lambda_{\min}(\mathbf{A}) = \frac{S}{\Delta t},\qquad \kappa = 1 + \frac{\Delta t}{S}\,\lambda_{\max}(\mathbf{C})

where:

$\mathbf{A}$ matrix of the pressure system, m³/(d·MPa)
$\mathbf{C}$ matrix of links: $\sum_j a_{ij}$ on the diagonal, $-a_{ij}$ off it, m³/(d·MPa)
$\mathbf{I},\ \mathbf{1}$ identity matrix and the vector of ones, an equal change of pressure in all cells, —
$S$ storage capacity of a cell (the same for all cells), m³/MPa
$\lambda_{\min},\ \lambda_{\max}$ smallest and largest eigenvalues, m³/(d·MPa)
$\kappa$ condition number $\lambda_{\max}(\mathbf{A})/\lambda_{\min}(\mathbf{A})$ , the ratio of the largest to the smallest eigenvalue; it shows how sensitive the solution is to errors, —

Units: metric oilfield units; $\kappa$ is dimensionless

Holds when

The cells have the same capacity. The outer boundaries are closed, and the wells work at given rates.
A well with a given bottomhole pressure, or a constant-pressure boundary, adds to the diagonal and increases $\lambda_{\min}$ .

The links do not “see” a uniform shift of pressure; only the capacity determines it. The larger the step, the smaller $S/\Delta t$ and the worse the conditioning of the system. In the limit of an incompressible fluid with given rates, the system is degenerate (it has no single solution): the pressure is known only up to a constant, because adding the same number to every pressure gives another solution. It is exactly such smooth components of the error, shared by a large part of the model, that are hardest to remove by simple iterations.

SchematicThe block of four cells of the examples (a sketch): the permeabilities of the cells, the transmissibilities of the connections and the wells, and the structure of the matrix of the system, with the non-zero elements filled in. After

\delta_1

is eliminated, new non-zero elements appear at the places (2, 3) and (3, 2). This is fill-in.

Worked exampleThe system of equations for a block of four cells

Given (assumed data)

Block: 2 × 2 cells of 300 × 300 × 12 m: 1 is south-west, 2 is south-east, 3 is north-west, 4 is north-east; the outer faces are closed
Permeability of cells 1–4: 120; 40; 200; 80 mD
Rock and liquid: $\phi$ = 0.20; $c_t = 1.5\cdot10^{-3}$ MPa⁻¹; viscosity 1.0 mPa·s
Wells and step: injection of 400 m³/d into cell 1, withdrawal of 400 m³/d from cell 4 (reservoir conditions); implicit step $\Delta t$ = 10 d from a uniform pressure

Solution

The transmissibilities are two half-cells in series (transmissibility of a connection). For cells of the same size, with a face of $A = 300\cdot12 = 3600$ m² and half a cell of 150 m, $T = 3600/(150/k_i + 150/k_j)$ : $T_{12} = 720$ , $T_{34} = 1371.43$ , $T_{13} = 1800$ , $T_{24} = 640$ mD·m.
Link coefficients $a = 0.0853\,T/1.0$ : $a_{12} = 61.416$ , $a_{34} = 116.983$ , $a_{13} = 153.540$ , $a_{24} = 54.592$ m³/(d·MPa). The storage capacity of a cell is $S = 300\cdot300\cdot12\cdot0.20\cdot1.5\cdot10^{-3} = 324$ m³/MPa, and $S/\Delta t = 32.4$ .
The diagonal is $S/\Delta t$ plus the link coefficients of the cell: 247.356, 148.408, 302.923 and 203.975.
The cell equations (the pressures at the start of the step are equal, so only the wells stand on the right): 1: $247.356\,\delta_1 - 61.416\,\delta_2 - 153.540\,\delta_3 = 400$ ; 2: $-61.416\,\delta_1 + 148.408\,\delta_2 - 54.592\,\delta_4 = 0$ ; 3: $-153.540\,\delta_1 + 302.923\,\delta_3 - 116.983\,\delta_4 = 0$ ; 4: $-54.592\,\delta_2 - 116.983\,\delta_3 + 203.975\,\delta_4 = -400$ .
The eigenvalues of the matrix are 32.4, 176.4, 221.0 and 472.9. The smallest equals $S/\Delta t$ . The condition number is $472.9/32.4 = 14.6$ .

Answer: The 4 × 4 matrix is symmetric. In each row the diagonal is larger than the sum of the absolute values of the other numbers by 32.4. Four elements out of sixteen are zeros, because cells 1 and 4, and cells 2 and 3, do not share a face.

What it means. The zeros at the places (1, 4) and (2, 3) are a special case of the main property: a row holds only its own connections. In a model of a million cells a row has the same seven numbers, and there are almost a million zeros. The smallest eigenvalue equals $S/\Delta t$ , which is the uniform shift of pressure. With a step of 100 d it would be ten times smaller, and the condition number would grow almost ten times (see the table below). The data in this example are assumed for illustration.

Direct methods: Gaussian elimination and LU decomposition

Gaussian elimination solves a system in two passes. The forward pass eliminates the unknowns one by one. From all the equations below the first we subtract the first one, multiplied by numbers chosen so that the first unknown disappears. Then we do the same with the second equation, and so on to the last. The matrix becomes upper triangular (all numbers below the diagonal are zero). The backward pass is a substitution from the bottom up: the last equation contains one unknown, and each earlier one contains one new unknown.

A step of the forward pass of Gaussian elimination

m_{ki} = \frac{A_{ki}}{A_{ii}},\qquad A_{kl} \leftarrow A_{kl} - m_{ki}\,A_{il},\qquad b_k \leftarrow b_k - m_{ki}\,b_i\qquad (k > i)

where:

$A_{kl},\ b_k$ elements of the matrix and of the right-hand side that change in the elimination, in the units of the system
$i$ number of the unknown being eliminated (the pivot row: the row that we subtract from the rows below it), —
$m_{ki}$ multiplier of row $k$ ; the multipliers form the lower triangular matrix $\mathbf{L}$ of the decomposition, —

Units: in the units of the system

Holds when

The pivot element $A_{ii}$ (the diagonal number that we divide by) is not zero. With diagonal dominance this is guaranteed, and no row swaps are needed.
Rows in which $A_{ki} = 0$ do not change.

After all the unknowns are eliminated, the matrix is upper triangular ( $\mathbf{U}$ ), and the multipliers form a lower triangular matrix $\mathbf{L}$ with ones on the diagonal: $\mathbf{A} = \mathbf{L}\mathbf{U}$ . This is the LU decomposition. We do it once, and then for any right-hand side we solve two triangular systems. For a symmetric positive definite matrix we use the Cholesky decomposition $\mathbf{A} = \mathbf{G}\mathbf{G}^{\mathsf T}$ , which is twice as cheap.

Back substitution

x_n = \frac{\tilde b_n}{U_{nn}},\qquad x_i = \frac{1}{U_{ii}}\left(\tilde b_i - \sum_{l > i} U_{il}\,x_l\right),\quad i = n - 1,\ \ldots,\ 1

where:

$U_{il}$ elements of the upper triangular matrix after the forward pass, in the units of the system
$\tilde b_i$ right-hand side after the forward pass, in the units of the system
$x_i$ unknowns, in the units of the system

Units: in the units of the system

Holds when

An upper triangular system with a non-zero diagonal.

For a tridiagonal matrix the two passes make up the Thomas algorithm: the number of operations is proportional to the number of unknowns.

Worked exampleGaussian elimination for the block of four cells

Given (assumed data)

System: from the previous example; the coefficients have three decimals, the multipliers five

Solution

We eliminate $\delta_1$ from row 2. The multiplier is $61.416/247.356 = 0.24829$ . We add row 1 multiplied by it to row 2. The diagonal becomes $148.408 - 61.416\cdot0.24829 = 133.159$ , the right-hand side becomes $400\cdot0.24829 = 99.316$ , and at $\delta_3$ the zero is replaced by $-153.540\cdot0.24829 = -38.122$ . This is the fill-in.
Row 3: the multiplier is $153.540/247.356 = 0.62072$ . The diagonal becomes $302.923 - 153.540\cdot0.62072 = 207.618$ , at $\delta_2$ a new number appears, $-61.416\cdot0.62072 = -38.122$ , and the right-hand side becomes $400\cdot0.62072 = 248.288$ . The elimination of $\delta_1$ does not change row 4, because it has no $\delta_1$ .
We eliminate $\delta_2$ with the help of the new row 2. Row 3: the multiplier is $38.122/133.159 = 0.28629$ . The diagonal becomes $207.618 - 38.122\cdot0.28629 = 196.704$ , the coefficient at $\delta_4$ becomes $-(116.983 + 54.592\cdot0.28629) = -132.612$ , and the right-hand side becomes $248.288 + 99.316\cdot0.28629 = 276.721$ .
Row 4: the multiplier is $54.592/133.159 = 0.40998$ . The diagonal becomes $203.975 - 54.592\cdot0.40998 = 181.593$ , the coefficient at $\delta_3$ becomes $-(116.983 + 38.122\cdot0.40998) = -132.612$ , and the right-hand side becomes $-400 + 99.316\cdot0.40998 = -359.282$ .
We eliminate $\delta_3$ from row 4: the multiplier is $132.612/196.704 = 0.67417$ , the diagonal becomes $181.593 - 132.612\cdot0.67417 = 92.190$ , and the right-hand side becomes $-359.282 + 276.721\cdot0.67417 = -172.725$ . The matrix has become upper triangular.
Back substitution: $\delta_4 = -172.725/92.190 = -1.87358$ ; $\delta_3 = (276.721 + 132.612\,\delta_4)/196.704 = 0.14368$ ; $\delta_2 = (99.316 + 38.122\,\delta_3 + 54.592\,\delta_4)/133.159 = 0.01886$ ; $\delta_1 = (400 + 61.416\,\delta_2 + 153.540\,\delta_3)/247.356 = 1.71097$ MPa.
Check: the solution without rounding is 1.71099, 0.01887, 0.14370 and −1.87356 MPa. The rounding in the chain changed the answer by not more than 0.00002 MPa.

Answer: Over the step the pressure in the cell with injection grew by 1.71 MPa, and in the cell with withdrawal it fell by 1.87 MPa. In cells 2 and 3 it changed by 0.019 and 0.14 MPa.

What it means. By eliminating cell 1 we linked its neighbours, cells 2 and 3, which do not share a face: a new non-zero element appeared in the matrix. This always happens, and in a large grid the fill-in occupies the whole band (see the figure below). Note $\delta_2$ : it is a small difference of large numbers, so a calculation with few digits would lose accuracy. A computer keeps about 16 digits, but a badly conditioned system loses more of them the larger $\kappa$ is. The data in this example are assumed for illustration.

CalculatedThe matrix of the pressure equation of a 4 × 3 grid with natural numbering (along rows of 4 cells), and the multipliers of its LU decomposition. The non-zero elements of the matrix lie on five diagonal lines, and the elimination fills almost everything between the outer ones: a band as wide as a row of the grid (36 new elements). With numbering along the short side (3 cells in a row) the band is narrower and the fill-in is smaller, 24 elements.

Fill-in and numbering. The elimination of an unknown links all its neighbours that have not been eliminated yet to each other. So in a grid with natural numbering almost the whole band between the outer non-zero diagonals fills in. The work and the memory grow with its width, and the order of numbering matters: we number first along the shortest direction. For a grid of 10 × 5, numbering along the long side gives 648 new elements, and numbering along the short side gives 288. The inverse of a sparse system is almost completely full, so we never build it (we do not invert the Jacobian matrix).

Work and memory of banded elimination

W \approx n\,b^2,\qquad M \approx n\,b;\qquad b = N_x\ \ (\text{2D}),\qquad b = N_xN_y\ \ (\text{3D})

where:

$W$ number of multiply–add operations, —
$M$ number of stored numbers of the Cholesky decomposition (for the LU decomposition it is twice as many: $\mathbf{L}$ and $\mathbf{U}$ ), —
$n$ number of unknowns, —
$b$ half-width of the band: the largest difference between the numbers of linked unknowns, —
$N_x,\ N_y$ numbers of cells along the axes by which the numbering goes (fastest along $x$ ), —

Units: dimensionless counts

Holds when

The leading term of the estimate for $b \ll n$ . Elimination without row swaps.
For the Cholesky decomposition of a symmetric matrix the work is about half as large. For a dense matrix the work is $n^3/3$ .

The work grows as the square of the band width, so in a three-dimensional model we number first along the vertical, which is usually the shortest direction. Even then the band width is in the thousands, and for models of hundreds of thousands of cells banded elimination is too expensive (an example is in the last section).

For sparse matrices there are also more economical direct methods. The unknowns are renumbered so that the fill-in is smaller (the minimum degree and nested dissection methods), and only the non-zero numbers are stored. In two-dimensional problems they work well. In three dimensions their work grows about as the square of the number of unknowns, and for models of millions of cells, which are solved thousands of times in a run, they are too expensive. In simulators direct methods remain for small systems (two-dimensional cross-sections, well equations, the coarse levels of multigrid methods) and as the basis of preconditioners.

Simple iterative methods: Jacobi, Gauss–Seidel and over-relaxation

An iterative method starts with a first guess, for example zero changes of pressure (the pressures at the start of the step), and improves it again and again. The simplest idea is this. We solve each equation for “its own” unknown, put in the values of the neighbours from the previous guess, and repeat until the changes become small.

The Jacobi, Gauss–Seidel and successive over-relaxation methods

\begin{aligned} &\text{Jacobi:} && x_i^{(\nu+1)} = \frac{1}{A_{ii}}\Bigl(b_i - \sum_{l \ne i} A_{il}\,x_l^{(\nu)}\Bigr)\\ &\text{Gauss–Seidel:} && x_i^{(\nu+1)} = \frac{1}{A_{ii}}\Bigl(b_i - \sum_{l < i} A_{il}\,x_l^{(\nu+1)} - \sum_{l > i} A_{il}\,x_l^{(\nu)}\Bigr)\\ &\text{SOR:} && x_i^{(\nu+1)} = x_i^{(\nu)} + \omega\,\bigl(\hat x_i - x_i^{(\nu)}\bigr)\end{aligned}

where:

$x_i^{(\nu)}$ approximation of the unknown $i$ at iteration $\nu$ , in the units of the unknowns
$A_{il},\ b_i$ elements of the matrix and of the right-hand side, in the units of the system
$\hat x_i$ the value that the Gauss–Seidel formula would give, in the units of the unknowns
$\omega$ over-relaxation parameter, $1 < \omega < 2$ , —

Units: in the units of the system

Holds when

$A_{ii} \ne 0$ ; Gauss–Seidel and SOR go through the equations in the order of their numbers and use the new values at once;
convergence is guaranteed with strict diagonal dominance (Jacobi and Gauss–Seidel), and for a symmetric positive definite matrix for Gauss–Seidel and for SOR with $0 < \omega < 2$ .

Jacobi recalculates all the unknowns from the values of the previous iteration. Gauss–Seidel uses the values that are already improved. Successive over-relaxation (SOR) deliberately “overshoots”: it takes a step $\omega$ times as long as the Gauss–Seidel correction. An iteration costs about one multiplication of the matrix by a vector.

Worked exampleIterations for a block of four cells

Given (assumed data)

System: the 2 × 2 block of the first example, a step of 10 d
First guess: $\delta = 0$ , the pressures at the start of the step
Stopping criterion: relative residual $\lVert\mathbf{b} - \mathbf{A}\boldsymbol\delta\rVert/\lVert\mathbf{b}\rVert < 10^{-6}$

Solution

Jacobi, iteration 1: in equations 1 and 4 the neighbours are still zero, so $\delta_1 = 400/247.356 = 1.6171$ and $\delta_4 = -400/203.975 = -1.9610$ . In equations 2 and 3 the right-hand sides are zero, and $\delta_2 = \delta_3 = 0$ .
Jacobi, iteration 2: $\delta_2 = (61.416\cdot1.6171 + 54.592\cdot(-1.9610))/148.408$ , that is $(99.316 -107.055)/148.408 = -0.0521$ ; $\delta_3 = (248.290 -229.404)/302.923 = 0.0623$ . The values $\delta_1$ and $\delta_4$ did not change, because their neighbours at iteration 1 were zero.
Gauss–Seidel, iteration 1: $\delta_1 = 1.6171$ , as in Jacobi. Then the new $\delta_1$ goes at once into equations 2 and 3: $\delta_2 = 99.316/148.408 = 0.6692$ , $\delta_3 = 248.290/302.923 = 0.8196$ . Equation 4 uses the new $\delta_2$ and $\delta_3$ : $\delta_4 = (-400 + 36.533 + 95.879)/203.975 = -1.3119$ .
The spectral radius of the Jacobi iteration matrix is $\rho_J$ = 0.858. For Gauss–Seidel it is $\rho_J^2$ = 0.736. The optimal parameter $\omega = 2/(1 + \sqrt{1 - \rho_J^2}) = 1.32$ gives $\rho$ = 0.32.
Iterations to a residual of $10^{-6}$ : Jacobi needs 65, Gauss–Seidel 39, over-relaxation with $\omega$ = 1.32 needs 15, and the conjugate gradient method (next section) needs 3.

Answer: Jacobi needs 65 iterations, Gauss–Seidel 39, over-relaxation 15 and conjugate gradients 3.

What it means. On the first iteration Gauss–Seidel “overshot”: $\delta_2$ and $\delta_3$ received 0.6692 and 0.8196 MPa, much more than the solution. The reason is that in equations 2 and 3 it already saw the new $\delta_1$ but did not yet see the negative $\delta_4$ . We do not compare the methods by their first iterations: the speed in the long run is set by the spectral radius, and Gauss–Seidel then converges twice as fast as Jacobi. Conjugate gradients finished in 3 iterations: for such a small system a Krylov method is exact after a number of steps that is not larger than the number of unknowns. The data in this example are assumed for illustration.

Why do the iterations settle on the answer (we say that they converge), and how fast? At every iteration the error of the approximation, $\mathbf{e}^{(\nu)} = \mathbf{x} - \mathbf{x}^{(\nu)}$ , is multiplied by the same iteration matrix $\mathbf{G}$ (see the advanced section). The iterations converge from any first guess if the spectral radius $\rho(\mathbf{G})$ is smaller than one. The spectral radius is the largest absolute value of an eigenvalue of the iteration matrix, that is, the largest factor by which the matrix can stretch a pattern of the error. In the long run the error shrinks by a factor of $1/\rho$ at every iteration:

The speed of convergence of simple iterations

\nu \approx \frac{\ln(1/\varepsilon)}{-\ln\rho(\mathbf{G})},\qquad \rho_{GS} = \rho_J^2,\qquad \omega_{opt} = \frac{2}{1 + \sqrt{1 - \rho_J^2}},\qquad \rho_{SOR} = \omega_{opt} - 1

where:

$\nu$ number of iterations needed to reduce the error by a factor of $1/\varepsilon$ , —
$\rho(\mathbf{G})$ spectral radius of the iteration matrix of the method, —
$\rho_J,\ \rho_{GS},\ \rho_{SOR}$ spectral radii of the Jacobi, Gauss–Seidel and over-relaxation methods, —
$\omega_{opt}$ optimal over-relaxation parameter, —

Units: dimensionless

Holds when

an estimate for a large number of iterations; at the first iterations the speed depends on the initial error;
the relations for Gauss–Seidel and SOR hold for “consistently ordered” matrices, a technical condition on the order of the unknowns. The five-point matrices of Cartesian grids with natural numbering satisfy it.

Gauss–Seidel converges twice as fast as Jacobi, and over-relaxation with the optimal $\omega$ is several times faster than Gauss–Seidel. But $\rho_J$ tends to one when the grid is refined and when the time step grows. The errors that are smooth over the scale of many cells die out the slowest.

The 2 × 2 block: time step, conditioning and number of iterations

Step $\Delta t$ , d: 1
$\rho_J$: 0.387
$\kappa$: 2.36
Jacobi; Gauss–Seidel; SOR: 12; 8; 6
Conjugate gradients: 3

Step $\Delta t$ , d: 10
$\rho_J$: 0.858
$\kappa$: 14.6
Jacobi; Gauss–Seidel; SOR: 65; 39; 15
Conjugate gradients: 3

Step $\Delta t$ , d: 100
$\rho_J$: 0.984
$\kappa$: 137
Jacobi; Gauss–Seidel; SOR: 464; 287; 46
Conjugate gradients: 3

Step $\Delta t$ , d	$\rho_J$	$\kappa$	Jacobi; Gauss–Seidel; SOR	Conjugate gradients
1	0.387	2.36	12; 8; 6	3
10	0.858	14.6	65; 39; 15	3
100	0.984	137	464; 287; 46	3

Iterations from the zero first guess to a relative residual of $10^{-6}$ ; SOR with the optimal $\omega$ (1.04; 1.32; 1.69). The quantity $\kappa - 1$ grows exactly in proportion to the step, as the formula for the uniform shift of pressure predicts, and the number of simple iterations grows with it.

On a large grid the simple methods have an even harder time. In one Jacobi iteration a correction travels from a cell only to its neighbours. To fix a smooth error, for example the common pressure level in a large block of the formation, the information must pass across the whole grid. So the number of iterations of Jacobi and Gauss–Seidel grows as the square of the number of cells along a side, and that of over-relaxation grows in direct proportion to it (the table in the next section). Simulators used block variants for a long time, for example line successive over-relaxation (line SOR, LSOR): a whole row of cells is solved at once by the Thomas algorithm. Today simple iterations work mostly inside stronger methods, as smoothers (steps that remove the rough, cell-to-cell part of the error) in multigrid methods and in preconditioners.

Krylov subspace methods and preconditioning

Simple iterations use only the latest approximation. Krylov subspace methods remember the history. They look for the approximation at iteration $\nu$ as the best combination of the vectors $\mathbf{b}$ , $\mathbf{A}\mathbf{b}$ , …, $\mathbf{A}^{\nu - 1}\mathbf{b}$ , which we get by multiplying the matrix by a vector. For a symmetric positive definite matrix, such as the pressure equation, this is the conjugate gradient method (CG). It minimises the error in the “energy” norm (a way of measuring the error that takes the matrix into account) and, in exact arithmetic with no rounding, finishes in no more than $n$ iterations. In the example of the 2 × 2 block it takes three: the right-hand side has no component of the uniform shift (injection equals withdrawal), and three directions are left. For non-symmetric Jacobian matrices we use GMRES and ORTHOMIN. They keep the history of directions, and its length is limited by dropping old directions or by starting again from time to time. We can also use BiCGStab, which manages with a short history but converges less smoothly.

The convergence of the conjugate gradient method

\lVert\mathbf{e}^{(\nu)}\rVert_{\mathbf{A}} \le 2\left(\frac{\sqrt\kappa - 1}{\sqrt\kappa + 1}\right)^{\nu}\lVert\mathbf{e}^{(0)}\rVert_{\mathbf{A}},\qquad \nu \lesssim \frac{\sqrt\kappa}{2}\,\ln\frac{2}{\varepsilon}

where:

$\mathbf{e}^{(\nu)}$ error of the approximation after $\nu$ iterations, in the units of the unknowns
$\lVert\mathbf{e}\rVert_{\mathbf{A}}$ energy norm $\sqrt{\mathbf{e}^{\mathsf T}\mathbf{A}\mathbf{e}}$ , —
$\kappa$ condition number of the matrix, —
$\varepsilon$ how many times the error must be reduced (the equation uses the reciprocal of that number), —

Units: dimensionless

Holds when

a symmetric positive definite matrix; exact arithmetic;
an upper bound: in practice convergence is often faster, especially if the eigenvalues are gathered in groups.

The number of iterations grows as the square root of the condition number. For the model problem it grows in proportion to the number of cells along a side, as with over-relaxation, but we do not have to choose a parameter. This is the point of preconditioning: to reduce $\kappa$ .

Preconditioning (preconditioner). Instead of $\mathbf{A}\mathbf{x} = \mathbf{b}$ we solve an equivalent system in which the matrix is closer to the identity matrix (the matrix with ones on the diagonal and zeros everywhere else):

The preconditioned system

\mathbf{M}^{-1}\mathbf{A}\,\mathbf{x} = \mathbf{M}^{-1}\mathbf{b},\qquad \mathbf{M} \approx \mathbf{A},\qquad \kappa\left(\mathbf{M}^{-1}\mathbf{A}\right) \ll \kappa(\mathbf{A})

where:

$\mathbf{M}$ the preconditioner: a rough copy of $\mathbf{A}$ , with which a system is cheap to solve, in the units of $\mathbf{A}$
$\kappa$ condition number, —

Units: in the units of the system

Holds when

the matrix $\mathbf{M}^{-1}$ is not built: at every iteration we solve the system $\mathbf{M}\mathbf{z} = \mathbf{r}$ ;
for the conjugate gradient method $\mathbf{M}$ must be symmetric positive definite.

A good preconditioner is a compromise: the closer $\mathbf{M}$ is to $\mathbf{A}$ , the fewer iterations we need, but the more each one costs. The extreme cases are $\mathbf{M} = \mathbf{A}$ (one iteration, but this is the direct solution) and the diagonal of $\mathbf{A}$ (cheap, but it helps little).

Incomplete LU decomposition (ILU): Gaussian elimination in which the fill-in is thrown away. In ILU(0) the multipliers are stored only at the places of the non-zero elements of $\mathbf{A}$ , and in ILU(k) part of the fill-in is kept as well. It is cheap and especially useful for heterogeneous models.
Nested factorization: a preconditioner invented for reservoir simulation models. It uses a nested order of the links along the axes of the grid (inside a line of cells, between lines, between planes), and it is widely used in simulators with Cartesian and corner-point grids.
Two-stage preconditioning CPR (constrained pressure residual) for the fully implicit scheme. First we build an equation for the pressure from the equations of a cell and solve it separately. This part of the problem is “elliptic”: every cell feels all the others, and the pressure error is smooth and shared by the whole model. Then the rest of the error, which is linked mostly to the saturations and is local, is removed by a cheap ILU for the whole system.
Multigrid methods (algebraic multigrid, AMG, builds its coarse levels straight from the matrix): a smooth error, which simple iterations remove slowly, is corrected on coarse grids, where it stops being smooth. The number of iterations hardly depends on the size of the grid, so AMG often solves the pressure equation at the first stage of CPR.

The model problem: number of iterations to a relative residual of

10^{-6}

Method: Jacobi
10 × 10: 258
20 × 20: 986
40 × 40: 3729
20 × 20, heterogeneous: 2575

Method: Gauss–Seidel
10 × 10: 130
20 × 20: 494
40 × 40: 1866
20 × 20, heterogeneous: 1297

Method: SOR with the optimal $\omega$
10 × 10: 28
20 × 20: 53
40 × 40: 102
20 × 20, heterogeneous: 92

Method: Conjugate gradients
10 × 10: 13
20 × 20: 33
40 × 40: 66
20 × 20, heterogeneous: 299

Method: Conjugate gradients + ILU(0)
10 × 10: 10
20 × 20: 16
40 × 40: 30
20 × 20, heterogeneous: 24

Method	10 × 10	20 × 20	40 × 40	20 × 20, heterogeneous
Jacobi	258	986	3729	2575
Gauss–Seidel	130	494	1866	1297
SOR with the optimal $\omega$	28	53	102	92
Conjugate gradients	13	33	66	299
Conjugate gradients + ILU(0)	10	16	30	24

A square block of the formation surrounded by an aquifer at constant pressure; an incompressible fluid, injection in the centre; iterations from the zero first guess. In the heterogeneous variant the permeability is log-normal (its logarithm follows a bell-shaped distribution) with a standard deviation of $\ln k$ equal to 2: it changes from cell to cell by more than 4 orders of magnitude. The optimal $\omega$ for SOR is 1.53; 1.73; 1.85; 1.82. The data are assumed for illustration.

The table shows three things. Each time the number of cells along a side is doubled, the number of iterations of Jacobi and Gauss–Seidel grows almost fourfold (Jacobi: 258 → 986 → 3729), and that of over-relaxation and conjugate gradients roughly twofold. Heterogeneity slows the simple methods 1.7–2.6 times, and conjugate gradients without preconditioning roughly 9 times. Their speed is set by the condition number, and contrasts in transmissibility between cells worsen it by orders of magnitude. Incomplete LU decomposition almost removes this difficulty: 24 iterations against 299. In the homogeneous problem ILU(0) cuts the number of iterations only 1.3–2.2 times and does not change its growth with the size of the grid (10 → 16 → 30). Multigrid methods remove that growth.

CalculatedThe relative residual of the 20 × 20 model problem (homogeneous formation) over the iterations: Jacobi, Gauss–Seidel, over-relaxation (SOR) with the optimal

\omega

, and conjugate gradients (CG) without preconditioning and with ILU(0). The vertical axis is the base-10 logarithm of the residual: −6 corresponds to

10^{-6}

. Jacobi needs about a thousand iterations to reach

10^{-6}

, and Gauss–Seidel about five hundred.

When to stop the iterations: residual, error and balance

The stopping criterion and the link between the residual and the error

\frac{\lVert\mathbf{r}^{(\nu)}\rVert}{\lVert\mathbf{b}\rVert} \le \varepsilon,\qquad \mathbf{r}^{(\nu)} = \mathbf{b} - \mathbf{A}\,\mathbf{x}^{(\nu)},\qquad \frac{\lVert\mathbf{x} - \mathbf{x}^{(\nu)}\rVert}{\lVert\mathbf{x}\rVert} \le \kappa\,\frac{\lVert\mathbf{r}^{(\nu)}\rVert}{\lVert\mathbf{b}\rVert}

where:

$\mathbf{r}^{(\nu)}$ residual of the approximation $\mathbf{x}^{(\nu)}$ (by how much the equations are still not satisfied), in the units of the right-hand side
$\varepsilon$ tolerance: the relative residual at which the iterations are stopped, —
$\kappa$ condition number, —
$\lVert\cdot\rVert$ norm (size) of a vector, for example the square root of the sum of the squares, —

Units: dimensionless ratios

Holds when

we can calculate the residual but not the error, because the exact solution is unknown; the inequality for the error is an upper bound.

A small relative residual guarantees a small error only in a well-conditioned system. In a badly conditioned system the pressures may be noticeably inaccurate with the same residual. Such systems come from a large step, large contrasts in transmissibility, or small cells next to large ones.

The residual of the pressure equation has a physical meaning. Each equation is the balance of a cell over the step in cubic metres per day, and the residual is the volume “lost” or “created” in the cell per day. The sum of the residuals over all cells is the error of the material balance of the step: the flows between cells cancel out when we add them (why the grid does not lose mass). So simulators write their convergence criteria through the normalised residuals of the cells and the total balance (convergence criteria).

Worked exampleWhat a loose tolerance costs: the balance of a block of four cells

Given (assumed data)

System: the 2 × 2 block, a step of 10 d; the Gauss–Seidel method from the zero first guess
Tolerances: relative residual $10^{-2}$ and $10^{-4}$

Solution

To reach a residual of $10^{-2}$ we need 9 iterations (the residual is 0.0088). The pressure changes are 1.7793; 0.0714; 0.2038; −1.8250 MPa against the exact 1.7110; 0.0189; 0.1437; −1.8736 MPa, so the error is up to 0.068 MPa.
The storage capacity multiplied by the sum of the changes is $S\sum\delta_i = 324\cdot0.2295 = 74$ m³. In the exact solution the sum is zero, because injection equals withdrawal. So the pressures that were not fully solved “created” 74 m³ of fluid, which is 1.9% of the 4000 m³ injected over the step.
To reach a residual of $10^{-4}$ we need 24 iterations. The pressure error is 0.0007 MPa, and the “created” volume is 0.8 m³, or 0.02% of the injection.

Answer: With a tolerance of

10^{-2}

the error of the pressures is up to 0.068 MPa, and the balance of the step is broken by 74 m³ (1.9% of the injection). With a tolerance of

10^{-4}

the error is up to 0.0007 MPa and the imbalance is 0.8 m³.

What it means. A residual of $10^{-2}$ sounds like “one percent”, but this is a percentage of the right-hand side, not of the pressures. The pressure change in cell 2, where the right-hand side is zero, came out 3.8 times larger than the exact one. The balance error builds up from step to step, and the simulator uses it to judge whether the equations are solved well enough. In the fully implicit scheme the tolerances of the linear solver are matched to the Newton iterations: at the first iterations it is enough to reduce the residual by one or two orders of magnitude, and at the last ones by more. The data in this example are assumed for illustration.

What it costs: why the solver decides the run time

The work of a direct and of an iterative solution

W_{\text{dense}} \approx \frac{n^3}{3},\qquad W_{\text{band}} \approx n\,b^2,\qquad W_{\text{iter}} \approx \nu\,c\,n

where:

$n$ number of unknowns, —
$b$ half-width of the band, —
$\nu$ number of iterations, —
$c$ operations per unknown per iteration: the multiplication of the matrix by a vector, the preconditioner, the dot products, —

Units: multiply–add operations

Holds when

leading terms of the estimates; for a seven-point matrix with ILU(0), $c$ is of the order of 20 to 30.

A direct method pays with a power of the size of the problem. An iterative method pays with the number of iterations, which also grows with the size, but more slowly.

Worked exampleA direct or an iterative solution for a three-dimensional model

Given (assumed data)

Grid: 100 × 100 × 20 cells, one pressure equation per cell: $n$ = 200,000
Iterative solver (assumption): 50 iterations of preconditioned conjugate gradients, about 25 operations per unknown per iteration
Computer (assumption): $10^9$ multiply–add operations per second on one core
Number of solutions per run (assumption): 4,000: hundreds of steps with several iterations each

Solution

A dense matrix: $n^3/3 = 2.67\cdot10^{15}$ operations, which is 31 d for one solution.
A band with natural numbering: $b = 100\cdot100 = 10{,}000$ . The work is $nb^2 = 2\cdot10^{13}$ , which is 5.6 h. The memory is $nb = 2\cdot10^{9}$ numbers of 8 bytes each, which is 16 GB.
Numbering along the vertical first: $b = 20\cdot100 = 2000$ . The work is $8\cdot10^{11}$ , which is 13 min per solution, and the memory is 3.2 GB.
An iterative solver: $50\cdot25\cdot200{,}000 = 2.5\cdot10^{8}$ operations, which is 0.25 s.
For 4000 solutions the band (with the better numbering) takes 37 d, and the iterative solver takes 17 min.

Answer: Banded elimination would take 13 min for one solution and 37 d for the run. The iterative solver takes 0.25 s and 17 min.

What it means. For three-dimensional models the choice is obvious, and the bigger the grid, the more obvious it is. The work of the banded method grows as $n^{7/3}$ (the width of the band is a layer of cells), and that of the iterative method roughly as $n$ times a slowly growing number of iterations. The speed of the computer and the number of operations here are assumptions, made to estimate the order of magnitude. The data in this example are assumed for illustration.

Why does the solver decide the run time? Building the system, which means calculating the properties, the flows and the derivatives, costs in proportion to the number of cells. Solving it costs the number of cells multiplied by the number of iterations. That number grows with the size of the grid, with contrasts in transmissibility, with the time step and with the nonlinearity of the problem. So in large models the linear solver usually takes the main share of the time, and its convergence matters no less than the choice of the grid.

If the solver converges slowly or does not converge

1
Read the log
see on which steps the number of linear iterations grows and where the step is cut. After a well starts or stops this is normal; on quiet stretches it is a warning sign.
2
Look for bad cells
a tiny pore volume next to large cells, transmissibilities that differ by many orders of magnitude, and non-neighbour connections with a huge transmissibility. They spoil the conditioning. Such cells are made inactive or are merged with their neighbours.
3
Check the wells
a large well index in a cell of small volume gives sharp changes of pressure and bad conditioning.
4
Change the settings with care
change the tolerances and the type of preconditioner only after comparing the balance and the results with a run at the previous settings.

Terms and symbols

Linear solver: The part of a simulator that solves the systems $\mathbf{A}\mathbf{x} = \mathbf{b}$ of the implicit steps and of the Newton iterations.
Sparse matrix: A matrix in which a row has as many non-zero numbers as the cell has links, plus the diagonal. All the other numbers are zeros.
Band width $b$: The largest difference between the numbers of linked unknowns. With natural numbering it is the number of cells in a row (2D) or in a layer (3D).
Fill-in: New non-zero elements that appear in the elimination: an eliminated unknown links its neighbours to each other.
Spectral radius $\rho$: The largest absolute value of an eigenvalue of the iteration matrix of an iterative method. The method converges if $\rho < 1$ .
Condition number $\kappa$: The ratio of the largest eigenvalue of a symmetric matrix to the smallest. It shows by how many times the relative error can exceed the relative residual.
Preconditioner: An easily solved rough copy of the matrix that reduces the number of iterations: ILU, nested factorization, CPR, multigrid methods.
Residual: $\mathbf{b} - \mathbf{A}\mathbf{x}$ . For the pressure equation, the imbalance of a cell in m³/d.

Common mistakes and misconceptions

Wrong

“To solve the system we must invert the matrix.”

In fact

The inverse of a sparse system is almost full: through the chain of neighbours each cell influences all the others, and for a model of a million cells that is a trillion numbers. We solve the system by elimination or by iterations, without building the inverse matrix.

Wrong

“The residual is small, so the pressures are accurate.”

In fact

The relative error can exceed the relative residual by a factor of $\kappa$ . In the example of the lesson, with a residual of $10^{-2}$ the pressure change in cell 2 came out almost four times larger than the exact one, and the balance of the step was broken by 1.9% of the injection.

Wrong

“If a run is slow, we can loosen the tolerance of the solver.”

In fact

A system that is not fully solved “creates” or “loses” fluid, and the balance error builds up from step to step. In the fully implicit scheme a loose tolerance also adds Newton iterations. We change the tolerance only after comparing the balance and the results.

Wrong

“Simple iterations suit any grid, we will just need more of them.”

In fact

Their number grows as the square of the number of cells along a side. In the 40 × 40 model problem Jacobi needs 3729 iterations, and conjugate gradients with ILU(0) need 30. Large models need Krylov methods with preconditioning.

In depth: advanced detail

Why simple iterations are slow on large grids

We write the matrix as $\mathbf{A} = \mathbf{D} - \mathbf{L} - \mathbf{U}$ : the diagonal, and the lower and upper parts with the opposite sign. The Jacobi method is $\mathbf{D}\mathbf{x}^{(\nu+1)} = \mathbf{b} + (\mathbf{L} + \mathbf{U})\,\mathbf{x}^{(\nu)}$ . We subtract it from the exact equality $\mathbf{D}\mathbf{x} = \mathbf{b} + (\mathbf{L} + \mathbf{U})\,\mathbf{x}$ :

The error of the iterations and the iteration matrix

\mathbf{e}^{(\nu+1)} = \mathbf{G}\,\mathbf{e}^{(\nu)},\qquad \mathbf{G}_J = \mathbf{D}^{-1}(\mathbf{L} + \mathbf{U}),\qquad \mathbf{G}_{GS} = (\mathbf{D} - \mathbf{L})^{-1}\,\mathbf{U}

where:

$\mathbf{e}^{(\nu)}$ error of the approximation $\mathbf{x} - \mathbf{x}^{(\nu)}$ , in the units of the unknowns
$\mathbf{D},\ \mathbf{L},\ \mathbf{U}$ the diagonal, and the lower and upper triangular parts of the matrix (with the opposite sign), in the units of the matrix
$\mathbf{G}_J,\ \mathbf{G}_{GS}$ iteration matrices of the Jacobi and Gauss–Seidel methods, —

Units: dimensionless

Holds when

a linear system with a constant matrix; a diagonal with no zeros.

After $\nu$ iterations the error is $\mathbf{G}^{\nu}\mathbf{e}^{(0)}$ : its components along the eigenvectors of $\mathbf{G}$ are multiplied by the eigenvalues. The component with the largest absolute eigenvalue, the spectral radius, dies out the slowest.

Take a one-dimensional strip of $N$ identical cells with a given pressure beyond the ends, in the simplest form (2 on the diagonal, −1 next to it). The eigenvectors of the Jacobi matrix are sine waves with $m = 1, \ldots, N$ half-waves, and the eigenvalues are $\cos[m\pi/(N + 1)]$ . The smoothest sine wave dies out the slowest: $\rho_J = \cos[\pi/(N + 1)] \approx 1 - \pi^2/[2(N + 1)^2]$ , and the number of iterations is $\nu \approx 2(N + 1)^2\ln(1/\varepsilon)/\pi^2$ , in proportion to the square of the number of cells. The jagged components, which change sign from cell to cell, are removed by the Gauss–Seidel method in a few iterations. Multigrid methods are built on this: simple iterations smooth the error, and its smooth part is corrected on a coarse grid, where it is no longer smooth.

How the conjugate gradient method works

The conjugate gradient algorithm

\begin{aligned} &\mathbf{r}^{(0)} = \mathbf{b} - \mathbf{A}\mathbf{x}^{(0)},\quad \mathbf{p}^{(0)} = \mathbf{r}^{(0)};\\ &\alpha_\nu = \frac{\mathbf{r}^{(\nu)\mathsf T}\mathbf{r}^{(\nu)}}{\mathbf{p}^{(\nu)\mathsf T}\mathbf{A}\,\mathbf{p}^{(\nu)}},\quad \mathbf{x}^{(\nu+1)} = \mathbf{x}^{(\nu)} + \alpha_\nu\,\mathbf{p}^{(\nu)},\quad \mathbf{r}^{(\nu+1)} = \mathbf{r}^{(\nu)} - \alpha_\nu\,\mathbf{A}\,\mathbf{p}^{(\nu)};\\ &\beta_\nu = \frac{\mathbf{r}^{(\nu+1)\mathsf T}\mathbf{r}^{(\nu+1)}}{\mathbf{r}^{(\nu)\mathsf T}\mathbf{r}^{(\nu)}},\quad \mathbf{p}^{(\nu+1)} = \mathbf{r}^{(\nu+1)} + \beta_\nu\,\mathbf{p}^{(\nu)}\end{aligned}

where:

$\mathbf{r}^{(\nu)}$ residual at iteration $\nu$ , in the units of the right-hand side
$\mathbf{p}^{(\nu)}$ search direction, in the units of the unknowns
$\alpha_\nu,\ \beta_\nu$ the length of the step along the direction and the coefficient of the new direction, —

Units: in the units of the system

Holds when

a symmetric positive definite matrix; with a preconditioner, in the directions and the dot products we use $\mathbf{z} = \mathbf{M}^{-1}\mathbf{r}$ in place of $\mathbf{r}$ .

One iteration takes one multiplication of the matrix by a vector and two dot products (each is a sum of products of the matching numbers of two vectors). The directions are conjugate: $\mathbf{p}^{(\mu)\mathsf T}\mathbf{A}\,\mathbf{p}^{(\nu)} = 0$ for $\mu \ne \nu$ , so a step along a new direction does not spoil what was reached along the earlier ones. In exact arithmetic, after $k$ iterations the error is the smallest possible in the Krylov subspace of dimension $k$ (the set of all combinations of the first $k$ vectors $\mathbf{b}$ , $\mathbf{A}\mathbf{b}$ , …). Rounding in the computer gradually destroys the conjugacy, and in practice the method works as an iterative one.

The systems of the fully implicit scheme and two-stage preconditioning

In the fully implicit scheme a cell has two or three unknowns, the pressure and the saturations, and the Newton matrix is made of blocks of 2 × 2 or 3 × 3 (Jacobian matrix). The parts of the problem behave differently. The pressure behaves almost “elliptically”: at an implicit step a change in the rate of a well is felt across the whole model, and the pressure error is smooth and shared. The saturations behave “hyperbolically”: the flow carries them, and they change mainly near fronts and wells, so their error is local.

Two-stage preconditioning CPR uses this difference. At the first stage the equations of the phases of a cell are combined into one equation for the pressure, much as the phase balances are added together in the pressure equation (pressure equation). This equation is solved by a strong method, usually algebraic multigrid. At the second stage a cheap incomplete LU decomposition of the whole system removes the local error that is left. Each stage does what it does best, and the number of iterations of the outer method (usually GMRES) stays small even for models of millions of cells.

Connections to other topics

Pressure equation for single-phase flow — the cell equations, the implicit scheme and the Thomas algorithm for a one-dimensional grid.
IMPES and fully implicit schemes — where the linear systems of the Newton iterations come from, and the convergence criteria.
Interblock flow and transmissibility — the transmissibilities, which are the coefficients of the matrix, and the non-neighbour connections.
Finite differences: from equations to computation — how the convergence of the iterations differs from the convergence of the difference scheme.
Grids in reservoir simulation models — the size and shape of the grid, on which the size of the system and the conditioning depend.
Two-phase flow equations — the pressure equation and the saturation equations of the two-phase model.

Check yourself

1
How many non-zero numbers are in a row of the pressure matrix for an interior cell of a three-dimensional Cartesian grid, with no non-neighbour connections and no wells?
1. A.7
2. B.3
3. C.27
4. D.5
Show the answer
Correct answer: A. 7
The diagonal and six neighbours across the faces: right, left, front, back, above and below. There would be 27 numbers if the cell were linked to all its neighbours across the corners and edges as well, as in some multi-point schemes.
2
A model has 200 × 100 × 30 cells. What is the band width with natural numbering (fastest along the long axis, then along the second axis, then by layers), and with numbering along the vertical first and then along the short horizontal axis? How many times does the work of banded elimination differ?
Show the answer
With natural numbering, $b = 200\cdot100 = 20{,}000$ . With numbering along the vertical first and then along the short axis, $b = 30\cdot100 = 3000$ . The work is proportional to $b^2$ : $(20{,}000/3000)^2 \approx 44$ . But even 3000 is too wide a band for six hundred thousand unknowns: $nb^2 \approx 5.4\cdot10^{12}$ operations for one solution.
3
The time step has been increased 10 times. How does the condition number of the pressure system of a closed reservoir with wells at given rates change, if the step is already large?
1. A.It does not change
2. B.It decreases about 10 times
3. C.It increases 100 times
4. D.It increases about 10 times
Show the answer
Correct answer: D. It increases about 10 times
The smallest eigenvalue is $S/\Delta t$ , so it decreases 10 times, and the largest hardly changes, because the links set it. In the example of the lesson, going from 10 to 100 d the condition number grew from 14.6 to 137, and Jacobi needed seven times more iterations.
4
The Jacobi method needed 600 iterations on a 30 × 30 grid, and conjugate gradients needed 50. About how many will each method need on a 60 × 60 grid of the same reservoir?
Show the answer
The number of Jacobi iterations grows as the square of the number of cells along a side: about $600\cdot4 = 2400$ . For conjugate gradients it grows as the square root of the condition number, that is, in proportion to the number of cells along a side: about 100. And each iteration on a 60 × 60 grid costs four times as much, so the total work grows 16 times and 8 times.
5
The relative residual of a solution is $10^{-3}$ , and the condition number of the system is about $10^4$ . What can we say about the relative error of the pressures?
1. A.It is not larger than $10^{-7}$
2. B.It can reach $10$ , that is, be larger than the solution itself
3. C.It is not larger than $10^{-3}$
4. D.It is zero, because the residual is small
Show the answer
Correct answer: B. It can reach $10$ , that is, be larger than the solution itself
The upper bound is $\kappa\cdot10^{-3} = 10$ : with bad conditioning a small residual does not guarantee even the right order of magnitude. In practice the error is usually smaller than the bound, but this is why, in badly conditioned problems, the tolerance is made stricter or the preconditioner is improved.
6
After the solution, the sum of the residuals of the pressure equations over all cells is +15 m³/d, and the step is 20 d. What does this mean for the material balance?
Show the answer
The sum of the residuals is the imbalance of the model: the flows between cells cancel in the sum, and the accumulation and the wells remain. Over the step the pressure equations “lost” or “created” $15\cdot20 = 300$ m³ of fluid, depending on the sign convention of the residual. This is the error of the material balance that the simulator prints. We compare it with the withdrawals and the pore volume, and if it is noticeable, we make the tolerance stricter.

Updated: 19 October 2026