Linear Solvers in OpenFOAM: A Technical Review and SIMPLE Convergence Study

1. Introduction

Computational Fluid Dynamics (CFD) is a part of fluid mechanics that uses numerical analysis to solve advection–diffusion–reaction type equations related to a wide range of applications involving fluid flow, and heat and mass transfer. Discretisation is a technique in numerical analysis that transforms a continuous problem into a discrete counterpart, which can be solved numerically by means of numerical solution methods.

It first requires that the geometry is divided into small, non-overlapping finite elements (or finite volumes), usually called cells. The discretisation method then transforms the partial differential equations into discrete algebraic equations, which are integrated over these discrete cells. After the discretisation of each transport equation, the resulting set of discrete equations needs to be linearised, which leads to a linear system of equations in the form of $Au=b$. The resulting systems of equations can be solved numerically by linear solvers.

Although every CFD simulation depends heavily on solving such linear systems, the underlying solution methods typically receive little attention in everyday CFD work and are generally treated as a ’black box’. Since the choice of linear solver does not affect the final discrete solution but only the speed at which it is reached, the common approach is to rely on default linear solver settings or to compensate slow convergence with additional computational hardware. However, different linear solvers exhibit very different performance characteristics. Understanding these strengths and weaknesses can lead to significant speed-ups, often far exceeding what can be achieved by simply increasing computational resources. This is especially true with OpenFOAM, where the linear solver options are very diverse.

This paper provides a detailed review of the linear solver theory in OpenFOAM, in order to improve the general understanding and appreciation of the mathematical models, and to aid in speeding up the convergence without the need for stronger hardware. Details are provided in Section 2. The theory in this work partially follows Refs. [1,2,3,4,5].

The performance of the solution methods in OpenFOAM is demonstrated for an industrial case in Section 3.

1.1. Decoupling the Velocity and Pressure in CFD

Before discussing the linear solvers used in this work, it is important to briefly revisit how the velocity and pressure fields are decoupled in CFD. Once the velocity field is determined, solving the transport of a particular scalar (such as energy and the variables for turbulence) is relatively easy, as the flow field of that scalar is ’frozen’, and the main challenge is usually to determine the source terms. However, to calculate the velocity field requires solving the discretised Navier-Stokes equations numerically, which is difficult. Not necessarily because it is non-linear, but due to the coupling of the velocity and pressure in the momentum equation. Consider the incompressible steady-state Navier-Stokes equations with neglected gravity

$$\nabla ·\mathit{U}=0,$$

(1)

$$\mathit{U}·\nabla \mathit{U}=-\nabla {\displaystyle \frac{p}{\rho }}+\nabla ·\nu \nabla \mathit{U},$$

(2)

where the vector $\mathit{U}$ (bold face) is the velocity, p is the pressure, $\rho$ is the density, and $\nu$ is the kinematic viscosity. These equations consist of 4 equations and 4 unknowns: $U_x$, $U_y$, $U_z$ and p (note that the density here is not a variable but a reference value). This would intuitively mean that it is directly solvable, even if there is no explicit equation for p. However, the continuity equation is not really the ’closing’ equation, but acts as a restriction to the momentum equations in order to satisfy mass conservation.

The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm decouples the pressure and velocity fields, allowing for calculating them separately, then using one field to correct the other field iteratively until convergence is achieved.

1.2. The SIMPLE Algorithm

The SIMPLE algorithm is developed by Patankar [6] to solve the steady-state incompressible Navier-Stokes equations. The main idea is to use the continuity equation to derive an equation for the pressure after splitting it from the momentum equations. The pressure equation will then act as a corrector for the modified momentum equation. It starts by by writing the momentum equations (Eq. 2) in general matrix form:

$$\mathit{A}\mathit{U}=-\nabla p,$$

(3)

where $\mathit{A}$ is the coefficient matrix resulting from discretising Eq. 2 using the Finite Volume Method. All these coefficients are known. Eq. 3 can be solved when the pressure is known from the previous iteration (starting from an initial guess at the first iteration), of which its gradient acts as the source term here. Also, the non-linearity appearing in the LHS of Eq. 2 is linearised by evaluating the mass fluxes using the velocity field from the previous iteration, resulting in a linear system for the velocity. The LHS of Eq. 3 is known in OpenFOAM as the velocity equation. The calculated velocity field does not yet satisfy continuity and is still a guess. This stage is called the momentum-predictor.

The derivation of the pressure equation starts from splitting the coefficient matrix $\mathit{A}$ by extracting its diagonal into matrix $\mathit{D}$:

$$\mathit{A}\mathit{U}=\mathit{D}\mathit{U}-\mathit{H}.$$

(4)

The reason is because a diagonal matrix can easily be inverted (${\mathit{D}}^{-1}$) which is very useful. $\mathit{H}$ is the residual vector containing the contributions of all the neighbour cells, and will be used as a source term for the pressure equation later. Combining Equations 3 and 4, and multiplying on both sides with ${\mathit{D}}^{-1}$ results in an equation for the velocity field which is an explicit function of the pressure field:

$$\mathit{U}={\mathit{D}}^{-1}\mathit{H}-{\mathit{D}}^{-1}\nabla p.$$

(5)

In the OpenFOAM source code, the diagonal matrix $\mathit{D}$ is defined as rAU, and the frequently recurring product ${\mathit{D}}^{-1}\mathit{H}$ is known as HbyA. Now we can use the continuity equation to derive an explicit equation for the pressure. Substituting Eq. 5 into Eq. 1 leads to following Poisson equation, which is known as the pressure corrector:

$$\nabla ·({\mathit{D}}^{-1}\nabla p)=\nabla ·\left({\mathit{D}}^{-1}\mathit{H}\right).$$

(6)

After solving for p in the above pressure equation, the velocity field (Eq. 5) can be corrected to satisfy continuity. Notice that the velocity is calculated explicitly and this is why this method is called Semi-Implicit. The problem now is that, once the velocity is corrected, the pressure equation is no longer satisfied because the matrix $\mathit{H}$ depends on the velocity which has been updated. Therefore this process is iterated until both the velocity and pressure fields no longer change after correction.

The SIMPLE algorithm can be extended to problems beyond steady incompressible flows. In such cases, most transport variables such as chemical species and enthalpy are solved after the momentum predictor step in OpenFOAM, and the turbulence transport equations are solved after the pressure corrector step.

2. Solution Methods

2.1. Introduction

For each transport equation, the discretisation and linearisation process leads to a system of linear equations, which is written as the following general matrix equation:

$$Au=b,$$

(7)

where A is the coefficient matrix resulting from the linearisation and mesh geometry, u is the vector that contains the unknown values of the dependent variable at the cell centroids which needs to be solved for, while b is the vector containing all the sources, constants and boundary conditions.

With the momentum-predictor in the previous section (Eq. 3), it is easy to see how it is translated into a general matrix equation:

• $\mathit{A}$$\equiv \mathit{A}$,
• $\mathit{U}$$\equiv u$,
• $-\nabla p$$\equiv b$,

while in the pressure-corrector (Eq. 6) the translation is a bit less apparent:

• ${\mathit{D}}^{*}$$\equiv \mathit{A}$,
• p$\equiv u$,
• $\nabla ·\left({\mathit{D}}^{-1}\mathit{H}\right)$$\equiv b$,

where ${\mathit{D}}^{*}$ is the resulting coefficient matrix after combining ${\mathit{D}}^{-1}$ with the discretised Laplace operator coefficients of p.

The linear systems are solved by linear solvers, for which different solution methods exist, which are generally grouped into direct and iterative methods, each consisting of many sub-groups. Due to the linearisation process with the finite volume method, the coefficient matrix is very large (many cells required for accuracy), very sparse, and diagonal dominant. Iterative methods therefore have been more popular because they are more suited for this type of applications as they typically require lower computational cost per iteration and lower memory. The (modified) direct methods can, however, be used as a preconditioner, in which they are used to replace the coefficient matrix by a matrix for which the corresponding linear system of equations is easier to solve, while keeping the same solution. The solution methods that are used in OpenFOAM are discussed here.

2.2. Direct Methods

Even though direct methods, as stand-alone solvers, are not efficient for solving sparse systems of linear equations due to their high computational cost, their discussion will lead the way for introducing efficient iterative methods in the next sections.

Direct methods apply some form of Gaussian elimination techniques to solve the above mentioned system of linear equations. The modern version is the LU factorisation, which computes the lower and upper triangular matrix L and U, such that

$$A=LU.$$

(8)

For diagonal dominant matrices, the upper triangular matrix U is constructed by performing Gaussian elimination on matrix A, while the elements of L consist of the Gauss factors by which the rows of A are multiplied to get U. Furthermore, the diagonal elements of L are set equal to one, so for a 3x3 matrix, Eq. 8 becomes

$$\begin{array}{c}\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ l_{21}& 1& 0\\ l_{31}& l_{32}& 1\end{array}\right]\left[\begin{array}{ccc}{u}_{11}& u_{12}& u_{13}\\ 0& u_{22}& u_{23}\\ 0& 0& u_{33}\end{array}\right].\end{array}$$

(9)

When L and U are found, the problem $Au=b$ is split into two problems that are much easier to solve

$$Ly=b,Uu=y.$$

(10)

To avoid rounding errors it might be necessary to apply some form of pivoting, i.e. reordering the rows of matrix A. However for diagonal dominant matrices this is not necessary [?]. If $\mathit{A}$ is Symmetric Positive Definite (SPD) the LU-factorisation reduces to the Cholesky factorisation of A.

$$A={CC}^{\mathsf{T}},$$

(11)

Where C is a lower triangular matrix. This results to even more memory and compuational savings as only one triangular matrix needs to be computed.

The efficiency of the direct methods however is lost due to the occurrence of fill-in after factorisation, where many zero entries in the sparse matrix become nonzero during the factorisation process. This significantly increases the memory requirements and computational cost, thereby disrupting the advantages of sparse matrix storage, especially for large practical CFD applications. This leaves an opening for iterative methods.

It will be shown that direct methods are very suitable as preconditioners for the Krylov subspace methods where the sparsity pattern is deliberately kept after factorisation, known as incomplete factorisation, so that $\mathit{A}\approx \mathit{LU}$ or $\mathit{A}\approx {\mathit{CC}}^{\mathsf{T}}$.

2.3. Basic Iterative Methods

The linear solvers in most CFD codes make use of iterative methods. Iterative methods are techniques created to obtain an approximate solution of linear systems. For the implementation of these methods, successive approximations are used. Starting from an initial guess $u_0$, a new approximation $u_k$ is obtained at each iteration k, until an approximated solution is found that is close enough to the exact solution u.

The error vector is defined as

$${\mathit{e}}_k=\mathit{u}-{\mathit{u}}_k.$$

(12)

Solving the error vector is just as difficult as solving the exact solution of the discrete linear system $Au$ = b. A computable measure of the quality of the approximation is therefore obtained from the residual,

$${\mathit{r}}_k=\mathit{b}-{\mathit{Au}}_k,$$

(13)

which represents the imbalance of the approximation of the conservation laws, with $r_k=0$ being the exact solution. A common stopping criterion or tolerance ($ϵ$) for iterative methods is the relative residual, defined as the 2-norm of the residual of the k-th iteration divided by the 2-norm of the right-hand side,

$$r_k={\displaystyle \frac{\left|\right|{\mathit{r}}_k{\left|\right|}_2}{{\left|\right|\mathit{b}\left|\right|}_2}}\le ϵ.$$

(14)

The idea of an iterative method is that the matrix A is decomposed into two matrices, M and N. Such that $A=M-N$; and the original linear system $Au=b$ transforms into:

$$\mathit{Au}=(\mathit{M}-\mathit{N})\mathit{u}=\mathit{b}.$$

(15)

rearranging terms we obtain:

$$\mathit{Mu}=\mathit{Nu}+\mathit{b}=(\mathit{M}-\mathit{A})\mathit{u}+\mathit{b}.$$

(16)

The latter system is used to perform an iterative process, finding at each iteration (k) a more accurate solution. Most of the iterative methods are derived from the following recurrence relation:

$${\mathit{u}}_{k+1}={\mathit{M}}^{-1}(\mathit{M}-\mathit{A}){\mathit{u}}_k+{\mathit{M}}^{-1}\mathit{b},$$

(17)

Or after rearranging and using Eq. 13

$${\mathit{u}}_{k+1}={\mathit{u}}_k+{\mathit{M}}^{-1}{\mathit{r}}_k,$$

(18)

where the matrix M is chosen such that $M^{-1}$ can be determined easily, so for example a diagonal or triangular matrix.

Basic iterative methods (BIM) are obtained by decomposing the system matrix as $A=D-E-F$ (see Figure 1) . D being the diagonal of A, and E and F are the strictly lower and upper parts. Based on the choice of M and N, different iterative methods can be obtained. Some of the BIMs are presented in Table 1, of which Jacobi and Gauss-Seidel are the options for OpenFOAM’s BIM solver smoothSolver.

Out of presented methods, the Successive Over-Relaxation method is significantly faster than the others for optimal damping parameters $\omega$, which value is not straightforward to find. Nevertheless, BIMs in general are too slow for engineering problems as stand alone solvers, as they are only effective in eliminating the high frequency components of the error, but very slow with reducing the low frequency components (Figure 2). This gets worse for finer meshes.

Nevertheless, it is the ability of eliminating the high frequency error components (smoothing) why the BIMs serve as the second set of building blocks for more efficient iterative methods, in addition to the direct methods. As smoothers for the multigrid method or as preconditioners for the Krylov subspace methods.

2.4. Multigrid Methods

The convergence of a BIM can be accelerated with a multigrid method, which is known to be among the most efficient solvers for elliptic problems [7] , such as the steady-state transport equation of a particular scalar $\varphi$

$$\underset{\mathrm{advection}}{\underbrace{\nabla ·\left(\rho \mathit{U}\varphi \right)}}=\underset{\mathrm{diffusion}}{\underbrace{\nabla ·\Gamma \nabla \varphi }}+\underset{\mathrm{source}}{\underbrace{Q_{\varphi }}}.$$

(19)

The basic idea of a multigrid method is to first smooth the error with a BIM, then transferring the remaining low frequency error components (recursively) to a coarser grid without losing information, until the problem is small enough to solve directly, after which a correction can be transferred back to the fine grid. As a result, the number of iterations required for convergence becomes nearly independent of the mesh element size, and the total computational cost scales approximately linearly with the number of unknowns. This high efficiency explains why multigrid methods are widely used in many CFD codes [3,8,9] . Depending on the size of the 2D or 3D mesh, the speed-up can be one or multiple orders of magnitude as compared to the BIM [2].

The detailed procedure is a follows:

1.

Presmoothing - Starting from an initial guess ${\mathit{u}}_0$, apply a few number ${\nu }_1$ of iterations with a BIM to eliminate the high frequency error components (recall from Eq. 18).

$${\mathit{u}}_n^h={\mathit{u}}_{n-1}^h+{\mathit{M}}^{-1}{\mathit{r}}_{n-1}^h,$$

(20)

where $n=1,..,{\nu }_1$ and the superscript h denotes the fine grid size.

2.

Restriction - To transfer the residual to a coarser grid of grid size H, we use the restriction operator $I_h^H$

$${\mathit{r}}^H=I_h^H{\mathit{r}}^h,$$

(21)

where superscript H denotes the coarse grid size.

3.

Solve coarse grid error - The residual can now be used to determine the coarse grid error by combining Eqs. 12 and 13

$${\mathit{A}}^H{\mathit{e}}^H={\mathit{r}}^H.$$

(22)

This is known as the residual equation. As mentioned before, solving this equation is just as hard as solving the original system of linear equations. However, at the current grid level the problem is small enough to solve it directly and the obtained error vector can be used as an approximation at the fine grid in the following step.

4.

Prolongation - The error $e^H$ is projected into the fine grid using the prolongation operator $I_H^h$ (inverse of $I_h^H$), and the vector $u_{{\nu }_1}^h$ is corrected with the approximated error.

$${\mathit{u}}_{{\nu }_1}^h={\mathit{u}}_{{\nu }_1}^h+I_H^h{e}^H$$

(23)

5.

Postsmoothing - Apply a few number of postsmoothing ${\nu }_2$ iterations as described in step 1 to obtain ${\mathit{u}}_1=:{\mathit{u}}_{{\nu }_1+{\nu }_2}^h$

In the above procedure, the sequence described by steps 2 to 4 is known as the defect correction, or coarse grid correction (CGC) and can be summarised by the following equation

$${\mathit{u}}_n^h=\left(\mathit{I}-I_H^h{\left({\mathit{A}}^H\right)}^{-1}{I}_h^H\mathit{A}\right){\mathit{u}}_{n-1}^h.$$

(24)

Although the BIM and the CGC are individually not efficient, the combination of the two as applied above is very efficient as they effectively reduce all the frequency components of the error. The described procedure above is essentially a two-grid cycle when going through the steps once. When steps 1-2 (hence also steps 4-5) are done recursively within one loop we arrive at the multigrid method. Most of the times this is necessary as the long wavelengths of the error modes become shorter in the coarser mesh and therefore presmoothing is needed to eliminate them quickly. Applying a certain number of the restriction-step will lead to a negligible cost of solving the residual equation compared to a smoothing sweep at the finest grid.

Applying step 3 only once is called a V-cycle, and applying it multiple times can result in a W-cycle or an F-cycle, which are more accurate but require also more work. A schematic is shown in Figure 3. There are many strategies and choices to tune the multigrid method, such as the amount of smoothing sweeps, grid-levels, and cycles, with varying trade-offs between speed of solving a single iteration and the overall rate of convergence.

One of the approaches to construct the coarser grid is by applying the restriction operator directly to the system matrix, which is known as the algebraic multigrid (AMG). This method is relatively easy to implement, but works as a black-box solver as it does not have any geometric interpretation of the mesh. The geometric multigrid (GMG) on the other hand obtains the coarser grid by clustering (or agglomerating) the cells of the mesh.

OpenFOAM works with the generalised method of geometric-algebraic multigrid (GAMG) which by default operates as a GMG solver, but also gives the option to enable the older AMG implementation (of its predecessor) instead [?]. The GMG solver agglomerates the cells by pairing each cell with an unpaired neighbour cell with the largest shared face area. This is repeated for every coarse grid level. An examples is shown in Figure 4. If a cells does not have a match, it will be added to the neighbouring group. The system matrix is adjusted accordingly.

The GAMG is a V-cycle method where the user can select a smoother (Gauss Seidel being most favoured) and adjust certain parameters, such as the mesh size at the coarsest level, the amount of presmoothing and postsmoothing sweeps. The solver calculates automatically the amount of intermediate coarse grid levels. A practical example of how the GAMG solver is set up is shown in Sec. Section 3.

Next to multigrid methods, another class of solvers called Krylov subspace methods can be adopted.

2.5. Krylov Subspace Methods

In this section a few other solvers are presented that are currently used in OpenFOAM, which are also one of the most efficient methods for solving large linear systems [1]. These methods are based on projection processes onto Krylov subspaces ${\mathcal{K}}_i$, that is

$${\mathcal{K}}_i(\mathit{A},{\mathit{r}}_0)=span\{{\mathit{r}}_0,{\mathit{Ar}}_0,{\mathit{A}}^2{\mathit{r}}_0,...,{\mathit{A}}^{i-1}{\mathit{r}}_0\},$$

(25)

where $\mathcal{K}$ is called the Krylov space of dimension i corresponding to matrix $\mathit{A}$ and initial residual ${\mathit{r}}_0$ (Eq. 13). The basic idea of Krylov subspace methods is that, for large and sparse matrices, it is more efficient to use A only in matrix–vector products. In this way, the approximate solution is constructed as a polynomial in A applied to the initial residual, without explicitly forming or inverting the matrix [10].

The Conjugate Gradient (CG) method is a Krylov subspace method that is developed to minimize the following quadratic equation

$$F\left(u\right)={\displaystyle \frac{1}{2}}{u}^{\mathsf{T}}Au-b^{\mathsf{T}}u,$$

(26)

where $\mathit{A}$ is an SPD matrix. Solving this minimization problem (which is finding ${\mathit{u}}^{*}$ such that $\nabla F\left(\mathit{u}\right)=0$ ) is equivalent to solving $Au=b$, and falls under the type of gradient methods (see Figure 5).

The solution ${\mathit{u}}^{*}$ is approached recursively via the relation

$$u_{k+1}=u_k+{\mathit{\alpha }}_k{d}_k,$$

(27)

where ${\mathit{d}}_k$ is the search direction for the next iteration. One way is to choose the residual vector as the search direction, which is known as the $\mathit{steepestdescent}$ method. However, this can lead to the same search directions being computed more than once and oscillations occurring around the local minima, leading to very slow convergence. Therefore in the CG method, every search direction should be in a unique direction. This is accomplished by selecting a set of search directions that are $\mathit{A}$-orthogonal (or $\mathit{A}$-$\mathit{conjugate}$) to the previous directions. This means that they satisfy the condition $d_k^{\mathsf{T}}A{d}_j=0$, for $k\ne j$.

The CG method starts with the residual vector being chosen as the first search direction,

$$d_0=r_0=b-{Au}_0.$$

(28)

The factor ${\alpha }_k$ in Eq. 27 ensures that the minimum point is found along the current search direction, which can be derived to the following expression

$${\alpha }_k={\displaystyle \frac{d_k^{\mathsf{T}}{r}_k}{d_k^{\mathsf{T}}A{d}_k}}.$$

(29)

After obtaining the new value for $\mathit{u}$ (Eq. 27), the new residual is calculated,

$$r_{k+1}=r_k-{\mathit{\alpha }}_k{Ad}_k.$$

(30)

To make sure that the next search direction is $\mathit{A}$-conjugate to the previous, the following coefficient is used which is simplified to,

$${\beta }_k={\displaystyle \frac{r_{k+1}^{\mathsf{T}}{r}_{k+1}}{r_k^{\mathsf{T}}{r}_k}},$$

(31)

so that finally the next search direction can be calculated, which completes the algorithm,

$$d_{k+1}=r_{k+1}+{\mathit{\beta }}_k{d}_k.$$

(32)

If A is SPD and the denominators of Eq. 29 and Eq. 31 are equal to zero, it means that the CG method breaks down when the problem is already solved (${\mathit{r}}_k$= 0), which makes the method robust. The rate of convergence of the CG methods depends on the spectral properties (i.e. the eigenvalue distribution) of matrix A and is related to the condition number of matrix $\mathit{A}$,

$${\kappa }_2(\mathit{A})={\displaystyle \frac{{\lambda }_{\max }\left(\mathit{A}\right)}{{\lambda }_{\min }\left(\mathit{A}\right)}},$$

(33)

where ${\lambda }_{max}$ and ${\lambda }_{min}$ are the largest and smallest eigenvalues of $\mathit{A}$, respectively. The more clustered the eigenvalues are, the smaller the value of ${\kappa }_2\left(\mathit{A}\right)$, hence faster convergence. For CFD applications, the condition number is about the square of the maximum number of grid points, hence the standard CG method is slow [2].

2.5.1. Preconditioners

Preconditioners are used to accelerate the convergence by transforming the problem into a similar one, but with a more clustered eigenvalue distribution. This is done by multiplying the original system of equations by the inverse of the preconditioned matrix $\mathit{P}$ as follows

$$P^{-1}Au=P^{-1}b.$$

(34)

In order for the preconditioned matrix $\mathit{P}$ to be effective it must fulfill the following requirements. $\mathit{P}$ should approximate $\mathit{A}$, the computation of ${\mathit{P}}^{-1}$ should be cheap, and the condition number of the transformed system should be smaller than that of the original system of equations,

$${\kappa }_2(P^{-1}\mathit{A})={\displaystyle \frac{{\mathit{\lambda }}_{\max }(P^{-1}\mathit{A})}{{\mathit{\lambda }}_{\min }(P^{-1}\mathit{A})}}<<{\mathit{\kappa }}_2(\mathit{A}).$$

(35)

For the CG method, $\mathit{P}$ must also be an SPD matrix. Therefore, a good and common choice is the Cholesky factorisation (Eq. 11), which yields $\mathit{P}={\mathit{CC}}^{\mathsf{T}}$. But due to fill-in, it takes large amounts of work and memory to construct $\mathit{C}$. Hence, the non-zero fill-in elements are disregarded, leading to the incomplete Cholesky factorisation of $\mathit{P}$.

The algorithm of the CG method (Eqs. 27 - 32) can be transformed into that of the preconditioned CG method by replacing ${\mathit{r}}_k$ with ${\mathit{z}}_k={\mathit{P}}^{-1}{\mathit{r}}_k$. The modification of the CG method is summarized in Algorithm 1.

Algorithm 1 Preconditioned Conjugate Gradient (PCG) method

$r_0=b-{Au}_0$ and $d_0=P^{-1}{r}_0$▹ Choose starting direction for $k=0,1,...,$ until convergence do     $z_k=P^{-1}{r}_k$     ${\alpha }_k={\displaystyle \frac{d_k^{\mathsf{T}}{z}_k}{d_k^{\mathsf{T}}A{d}_k}}$     $u_{k+1}=u_k+{\mathit{\alpha }}_k{d}_k$     $r_{k+1}=r_k-{\mathit{\alpha }}_k{Ad}_k$     ${\beta }_k={\displaystyle \frac{r_{k+1}^{\mathsf{T}}{z}_{k+1}}{r_k^{\mathsf{T}}{z}_k}}$     $d_{k+1}=z_{k+1}+{\mathit{\beta }}_k{d}_k$ end for

Due to its superlinear convergence and for being robust, the preconditioned CG (PCG) method is one of the best methods for solving the heat equation and the pressure corrector equation (Eq. 6), as the coefficient matrices $\mathit{A}$ that result from discretising these diffusion equations on the mesh are SPD.

For the general transport equations which contain the advection term, the coefficient matrix $\mathit{A}$ is not symmetric. Therefore the PCG method is not applicable. A significant amount of Krylov subspace methods are developed for general non-SPD matrices A, but none have all the following three nice properties which the CG method has for SPD matrices:

1.

It is a Krylov subspace method: ${\mathit{u}}_k\in {\mathcal{K}}_k(\mathit{A},{\mathit{r}}_0)$

2.

Optimality: the minimizer of the quadratic function (Eq. 26) exists, hence convergence is guaranteed.

3.

Short recurrences, hence requiring low computing work and storage.

2.5.2. General matrices

The Bi-Conjugate Gradient (BiCG) method is a Krylov subspace method for general matrices that transforms the unsymmetrical system into a symmetrical one so that it can by solved using a modified CG method. This is done by rewriting $\mathit{Au}=\mathit{b}$ as

$$\begin{array}{c}\left[\begin{array}{cc}0& \mathit{A}\\ {\mathit{A}}^{\mathsf{T}}& 0\end{array}\right]\left[\begin{array}{c}\widehat{\mathit{u}}\\ \mathit{u}\end{array}\right]=\left[\begin{array}{c}\mathit{b}\\ 0\end{array}\right],\end{array}$$

(36)

where $\widehat{\mathit{u}}$ is a dummy variable which is not solved for, but only used to convert the system. The update relations in the CG method for the residuals and search directions are augmented by their shadow equivalents which are based on ${\mathit{A}}^{\mathsf{T}}$,

$${\widehat{r}}_{k+1}={\widehat{r}}_k-{\alpha }_k{A}^{\mathsf{T}}{\widehat{d}}_k$$

(37)

and

$${\widehat{d}}_{k+1}={\widehat{r}}_{k+1}+{\beta }_k{\widehat{d}}_k.$$

(38)

The orthogonality of the residual and search direction with their shadow counterparts is ensured via the relation

$${\widehat{\mathit{r}}}_i^{\mathsf{T}}{\mathit{r}}_j={\widehat{\mathit{d}}}_i^{\mathsf{T}}{\mathit{Ad}}_j=0,i\ne j,$$

(39)

hence the name bi-orthogonal or bi-conjugate. The BiCG method is very close to the CG method and generates the same solution as the CG method if A is SPD. It does however require two matrix-vector multiplications (with $\mathit{A}$ and with ${\mathit{A}}^{\mathsf{T}}$) instead of one, such that every iteration is nearly twice as expensive as compared to the CG method.

Although the BiCG method shares two of the nice properties of the CG method, there is no optimality for general matrices $\mathit{A}$, hence convergence is not guaranteed. Moreover, the irregular convergence can lead to serious break-down due to large rounding errors. The Biconjugate Gradient STABilized (Bi-CGSTAB) method [11] is a more robust variant where the recurrence relations ${\tilde{\mathit{r}}}_k=P_k\left(\mathit{A}\right){\mathit{r}}_0$ are modified to a form

$${\tilde{\mathit{r}}}_k=Q_k\left(\mathit{A}\right)P_k\left(\mathit{A}\right){\mathit{r}}_0$$

(40)

such that the multiplication with ${\mathit{A}}^{\mathsf{T}}$ is avoided. For this recurrence, a $\mathit{k}$-th degree polynomial is taken of the form

$$Q_k\left(\mathit{A}\right)=(\mathit{I}-{\omega }_1\mathit{A})(\mathit{I}-{\omega }_2\mathit{A})...(\mathit{I}-{\omega }_k\mathit{A})$$

(41)

with suitable constants ${\omega }_k$ in the k-th iteration that minimizes residual $r_k$ with respect to ${\omega }_k$ , giving the method a semi-optimality property and smoother convergence behaviour as compared to the BiCG method. The Bi-CGSTAB method is currently the most recent and advanced Krylov subspace method that is implemented in the standard library of OpenFOAM, and can be preconditioned using e.g. the incomplete LU factorisation which is appropriate for general matrices. It is much faster than its predecessor, but sudden breakdown due to rounding errors still occur for the more difficult problems in this paper.

2.5.3. Other Krylov Subspace Methods

Although not used in this work, it is worth to mention the Generalized Minimal RESidual (GMRES) method [1], which falls under another class of Krylov subspace methods for solving general systems of linear equations. It is a stable method with respect to rounding errors, and has an optimality property which guarantees convergence. While the CG method aims at minimizing the A-norm of the error $\mathit{u}-{\mathit{u}}_k$ (only possible for SPD matrices), the GMRES method finds orthogonal vectors (search directions) that minimize the Euclidean norm of the residual ${\mathit{r}}_k$, leading to a comparable superlinear convergence behavior as the CG method. The main drawback however is that the whole sequence of search directions from the previous iterates have to be stored and multiplied with, such that the computational and memory requirements quickly become prohibitive if the method does not converge after a few iterations. Restarting the method after a certain number of iterations overcomes this limitation of long recurrences, but comes at the cost of losing the nice convergence properties due to throwing away all search information [12]. The generalized conjugate residual (GCR) method is similar to the GMRES method but with somewhat more floating point operations per iteration. However, the GCR method gives the ability to truncate to the last few search directions that are saved for the next iteration, which in return performs better than restarting the GMRES method.

Hybrid methods provide a balance between the short recurrences of the BiCG-type methods on one side and the optimality property of the GMRES-type methods on the other. First on the list is the GMRES Recursive (GMRESR) method [12] that consists of an inner and outer loop. The solution is first approximated in an inner-loop with the restarted GMRES method, after which the found search directions are condensed to the outer-loop, where they are used to approximate the solution with the truncated GCR method. The GMRESR method generally speeds up the convergence rate of the restarted GMRES method. Later improvements of the GMRESR-type methods are the GCRO-type methods, where the inner iteration loop takes place in a Krylov subspace orthogonal to the subspace of the outer loop (so-called subspace recycling), yielding further acceleration of the convergence rate[13] [14].

The final hybrid method worth mentioning, which leans more towards the BiCG-type methods, is the Induced Dimension Reduction (IDR) method. This method was already proposed before the Bi-CGSTAB method but has been revived to a more generalized variant IDR(s) [15]. The IDR(s) method is closely related to the Bi-CGSTAB method in the sense that the matrix-vector multiplication with ${\mathit{A}}^{\mathsf{T}}$ is avoided via a minimizing polynomial. The main difference is that the generated residuals are forced to be in subspaces of decreasing order until an s-dimensional space is reached. For the right value of s (typically $s\le$ 10.), this robust and efficient method is at least as fast as the Bi-CGSTAB method when s = 1 (original IDR method) and significantly faster for increased values of s, especially for more difficult problems. When s is large enough, the rate of convergence is nearly as good as full GMRES [16]. Hybrid methods such as the IDR(s) method and the GCRO-type methods are good candidates for the next generation of linear solvers for general matrices to implement in OpenFOAM.

3. Impact Linear Solvers on Convergence of the SIMPLE Algorithm

The speed of convergence towards a solution depends not only on the type of linear solvers used, but also on the type of problem that is being solved. This section evaluates the performance of the advanced linear solvers in OpenFOAM by conducting several numerical experiments on different cases.

3.1. Relative Tolerance

The solver algorithms of OpenFOAM are constructed in ’inner’ and ’outer’ loops. In the inner loops, the linear systems of equations are solved sequentially, and this is where the linear solvers do the work, which are referred to as ’inner’ iterations. The algorithm proceeds to the next linear system of equations in the inner loop once the current relative tolerance target is reached. In the outer loop, the flow develops to the next (pseudo) time step, which may or may not require a time integration, depending on whether the SIMPLE algorithm is used for steady-state problems, or the PISO algorithm for transient problems, or combinations thereof (PIMPLE). As the algorithms are iterative themselves, the linear solvers do not need to iterate to the smallest possible error in each loop before the algorithm progresses.

The 2D steady-state heat equation is a problem which does not require time integration or a flow to be developed. As such, this problem is solved within 1 outer iteration, and even up to 1 inner iteration with the most efficient solver, as is shown in Figure 6 and Figure 7.

In a fluid-related problem, however, the flow is physically developing on each outer iteration. Aiming for a low relative tolerance per inner loop can be more stable, but it also makes it much more time-consuming. The famous Pitz and Daily wind tunnel tutorial case [5] demonstrates that using a relative tolerance of ${10}^{-1}$ instead of ${10}^{-3}$ for all equations will lead to about twice as fast convergence, while a relative tolerance of ${10}^{-3}$ leads to smoother convergence with less oscillations, as is shown in Figure 8.

3.2. Numerical Experiment: Non-Reacting Flow

This section evaluates the performance of the linear solvers on a three-dimensional full scale, non-reacting and non-rotating kiln. The used solver is rhoSimpleFOAM, which is the standard solver for turbulent non-reacting steady-state compressible flows.

The geometry is shown in Figs. Figure 9 and Figure 10, where the air flow enters the domain at atmospheric pressure through the primary and secondary inlets with 0.15 kg/s and 1.6 kg/s, respectively, while the fuel inlet flow is excluded. Further boundary conditions are shown in Table 2. The mesh contains 2.3 million cells and all simulations are run for 5000 outer iterations on 20 cores in a single node.

Figure 11 shows the basic stream pattern in the longitudinal cross-section, where the annular shaped recirculation zone surrounding the burner is spotted.

3.2.1. Solver Combinations

The most advanced standard linear solvers of OpenFOAM have been tested on this applications, in different combinations which are shown in Table 3. When a Krylov subspace method is used, the solver name is preceded with a preconditioner. Instead of the IC and ILU factorisations, OpenFOAM implements the simplified diagonal-based variants as preconditioner, respectively the DIC and DILU, where off-diagonal terms are dropped. The reciprocal of the preconditioned diagonal is calculated and stored. The MG solver can also be used as a preconditioner. The settings for OpenFOAM’s MG solver (GAMG) is shown in Appendix A, and for these type of problems with complex shape and flows, the geometric variant of the GAMG method (agglomerator: faceAreaPair) is preferred, which is about twice as fast as the algebraic variant.

3.2.2. Results

Before discussing the linear solver performance, it is observed that by using the advanced solvers of OpenFOAM usually leads to 1 or 2 inner iterations per outer iteration for all variables, except with the pressure equation, which usually requires 1 or more orders of magnitude of inner iterations, before moving to the next outer iteration. Therefore the pressure-corrector equation is the major bottleneck that delays the solution to converge quickly. This occurs for most non-reacting flow cases. Figure 12 shows that the CG method with the DIC preconditioner is the least efficient method for the pressure equation, requiring 2 orders of magnitude of inner iterations on average to reach the relative target of ${10}^{-3}$. The MG method performs significantly better, and the most efficient method is the combination of the former two: the CG method with the MG as preconditioner.

Figure 13 shows the overall performance of the different linear solver combinations. Three critera are used to determine the speed of convergence:

• The time (in hours) to reach 5000 outer iterations,
• The amount of outer iterations to reach convergence,
• The time (in hours) to reach convergence.

The time needed to reach 5000 iterations with the different combinations is in line with the amount of required inner iterations for the pressure equation as is discussed above, since the pressure equation is the bottleneck.

When looking at the amount of outer iterations required to converge, the CG method with MG preconditioner also requires the least amount of iterations. This further accelerates the convergence time, which is a factor 2 to 3 times faster than the other solver combinations. However, if the solver is forced to a strictly compressible flow environment by setting the option transonic to on, the pressure equation becomes non-symmetric and the MG method will not work.

With regard to the other variables, the fastest solver is the Bi-CGSTAB method with the DILU preconditioner. The MG method, as solver or as preconditioner, is a bit slower because the other variables only require 1 or 2 inner iterations to reach the relative target of ${10}^{-3}$, and one MG operation is computationally more expensive than one Bi-CGSTAB operation, and much more expensive than one DILU operation.

The simulation can be accelerated further by relaxing the relative target from ${10}^{-3}$ to ${10}^{-1}$ , as is shown in Figure 14. This however may lead to less stability. In fact when using the DICCG-DILUBiCGSTAB method, the solution will even quickly diverge, as can be seen in Figure 15, while the MGCG-DILUBiCGSTAB method remains stable and the convergence is accelerated with 25 % as compared with the relative tolerance of ${10}^{-3}$ . Comparing these two extremes, with on one side the slowest of the selected methods that is restricted to a relative target of ${10}^{-3}$ , and on the other side the fastest method that is stable at a relative target of ${10}^{-1}$, a speed-up of nearly a factor 7 is reached. This speed up does not even consider the basic iterative solvers of OpenFOAM (smoothSolver).

Table 4. Simulation time of the different methods w.r.t time and outer iterations, using a relative tolerance of ${10}^{-3}$

Method	Time to 5000 iter	Convergence	Time to Convergence
MG-DILUBiCGSTAB	3.4 h	3000 iter	2.2 h
MG-MG	4.3 h	3500 iter	3.1 h
DICCG-MG	5.3 h	3000 iter	3.6 h
DICCG-DILUBiCGSTAB	6.4 h	5000 iter	6.4 h
MGCG-DILUBiCGSTAB	2.3 h	2500 iter	1.2 h
MGCG-MGBiCGSTAB	2.7 h	2500 iter	1.4 h

Table 4. Simulation time of the different methods w.r.t time and outer iterations, using a relative tolerance of ${10}^{-3}$

Method	Time to 5000 iter	Convergence	Time to Convergence
MG-DILUBiCGSTAB	3.4 h	3000 iter	2.2 h
MG-MG	4.3 h	3500 iter	3.1 h
DICCG-MG	5.3 h	3000 iter	3.6 h
DICCG-DILUBiCGSTAB	6.4 h	5000 iter	6.4 h
MGCG-DILUBiCGSTAB	2.3 h	2500 iter	1.2 h
MGCG-MGBiCGSTAB	2.7 h	2500 iter	1.4 h

Table 5. Simulation time of the different methods w.r.t time and outer iterations, using a relative tolerance of ${10}^{-1}$

Method	Time to 5000 iter	Convergence	Time to Convergence
MG-DILUBiCGSTAB	2.0 h	3000 iter	1.2 h
MG-MG	2.5 h	4000 iter	2.0 h
DICCG-MG	2.7 h	3000 iter	1.8 h
DICCG-DILUBiCGSTAB	–	Diverged	–
MGCG-DILUBiCGSTAB	1.5 h	3000 iter	0.9 h

Table 5. Simulation time of the different methods w.r.t time and outer iterations, using a relative tolerance of ${10}^{-1}$

Method	Time to 5000 iter	Convergence	Time to Convergence
MG-DILUBiCGSTAB	2.0 h	3000 iter	1.2 h
MG-MG	2.5 h	4000 iter	2.0 h
DICCG-MG	2.7 h	3000 iter	1.8 h
DICCG-DILUBiCGSTAB	–	Diverged	–
MGCG-DILUBiCGSTAB	1.5 h	3000 iter	0.9 h

3.2.3. Further Advancements

As discussed in the end of Sec. Section 2, faster linear solvers do exist than the ones provided by OpenFOAM. The state-of-the-art linear solvers and preconditioners can be found in the open source library Portable, Extensible Toolkit for Scientific computation (PETSc), and a plug-in into OpenFOAM is developed in [17]. Their work also discusses the very poor scalability of the MG preconditioner of OpenFOAM for massively parallel clusters (10³ cores), whereas the DIC preconditioner is shown to have outstanding superlinear scalability. It is reported that the DIC preconditioner overtakes the MG preconditioner in convergence time when using more than 2000 cores for a 3D laminar lid driven cavity flow problem consisting of 64 millions cells. The PETSc counterpart of OpenFOAM’s MG preconditioner solves the scalability issue.

4. Conclusion

This work has demonstrated the large differences in speed and robustness between the linear solvers available in OpenFOAM, with convergence times varying by more than an order of magnitude. The pressure-corrector equation was the major bottleneck that delayed the solution from converging quickly. Several of the more advanced standard linear solvers were tested with a 20-core high-performance computer on a 3D full-scale kiln, consisting of 2.3 million cells.

The most efficient and stable solution method turned out to be the Conjugate Gradient (CG) solver, combined with the Generalised Geometric-Algebraic MultiGrid (GAMG) as a preconditioner. A speed-up factor of 7 was reached as compared with the slowest of the advanced methods. Even higher speed-ups were reached compared with the basic linear solvers of OpenFOAM. Also, the stability is enhanced when the GAMG is used as preconditioner for the Bi-CG Stabilized (Bi-CGSTAB) method instead of the simple diagonal-based preconditioners.

Finally, due to the iterative nature of the SIMPLE algorithm, a faster linear solver than the MG-CG combination was redundant for this case as each time step required 1 or 2 iterations.