Linear Solvers#
After discretization in space and time of the Navier-Stokes equations, we obtain a problem that needs to be solved using a Netwon non-linear solver that solves a linear system in every iteration. In the case of a stabilized approach, we need to solve a linear system that looks as follows:
where the matrix to the left-hand side is non-symmetric. In the literature this is frequently called a generalized saddle point problem.
In Lethe we support the following linear solvers:
Trilinos direct solvers.
GMRES preconditioned with ILU or AMG.
BiCGStab preconditioned with ILU.
Only coupled methods are used at the moment: a direct method and two iterative algorithms (Krylov subspace methods). We are mainly interested in solving problems that are both large and sparse, therefore, the two iterative methods will be explained in this section, along with their preconditioners. A direct method should only be used for tests and development of new features as it is not efficient for large problems. For more information on the direct solver supported by Trilinos see the deal.II documentation: TrilinosWrappers::SolverDirect. This is not by any means a detailed explanation of all the linear solvers; however, it should give you a general idea and it should point you to useful references in case you are more interested in this topic.
General Theory on Krylov Methods#
Suppose that \(x_0\) is an initial guess for the solution of \(x\), and define the residual \(r_0=b - \mathcal{A} x_0\). Then, Krylov subspace methods are methods whose \(k\)-th iterate \(x_k\) satisfies:
where \(\mathcal{K}_k(\mathcal{A},x_0)\) is the \(k\)-th Krylov subspace generated by \(\mathcal{A}\) and \(x_0\). Because of the \(k\) degrees of freedom in the choice of the iterate \(x_k\), \(k\) constraints are required to make \(x_k\) unique. This is achieved by imposing that the \(k\)-th residual \(r_k = b - \mathcal{A} x_k\) is orthogonal to a \(k\)-dimensional space \(\mathcal{C}_k\), called the constraints space:
By knowing the properties of \(\mathcal{A}\), it is possible to determine constraints spaces \(\mathcal{C}_k\) that lead to uniquely defined iterates \(x_k\) for \(k=1,2,\dotsc\). The choice of \(\mathcal{C}_k\) and subsequent implementation leads to a different method.
Generalized Minimal Residual Method (GMRES)#
If \(\mathcal{A}\) is non-symmetric positive definite, an approximate solution can be found by choosing \(\mathcal{C}_k = \mathcal{K}_k(\mathcal{A}, r_0)\), one implementation of this corresponds to the GMRES method.
The GMRES method is based on what is called an Arnoldi iteration. This iteration has a \(k\)-term recurrence, which means that its computational cost increases as \(k\) increases. The method is not very practical for very large \(k\) because of the increase on memory and computational requirements. Therefore, for these kinds of systems, usually a restart is needed after certain iterations to build a new Krylov subspace. The number of iterations depends on the application. In Lethe the default number of iterations is 100, however, for the problems tackled in Lethe 200 or 250 are required. This restart may affect the convergence of the method when the matrix is not positive definite; a remedy for this is the use of a preconditioner.
In Lethe, we use the Trilinos implementation of the GMRES method through deal.II: TrilinosWrappers::SolverGMRES. It is possible to specify the minimum residual, the solver residual, the maximum number of iterations and the maximum Krylov vectors. The latter is in fact the one that controls the number of iterations after which a restart should be done. The first two are also known in literature as absolute residual (\(||b - \mathcal{A}x||\)) and relative residual (\(||b - \mathcal{A}x||/||b||\)), and are used to define the linear solver tolerance. The first one monitors the maximum absolute error allowed in a solution, while the second one is relative to the solution value; the linear solver tolerance is defined as the maximum value of both to avoid numerical issues.
Biconjugate Gradient Stabilized Method (BiCGStab)#
One may choose instead \(\mathcal{C}_k = \mathcal{K}_k(\mathcal{A}^T, r_0)\), which leads to the implementation of the BiCG method. However, this is not well defined for a general non-symmetric matrix \(\mathcal{A}\), therefore, a similar implementation with a stabilization step leads to the BiCGStab method. It avoids irregular convergence patterns by maintaining the convergence speed.
In the case of the BiCGStab method, the iteration is known as a transpose-free bi-Lanczos iteration. This kind of iteration has a fixed cost per iteration and no restart is ever needed. This method has a three-term recurrence, however, it needs two matrix-vector multiplications per iteration as it uses the system matrix twice to avoid the use of its transpose.
In Lethe we use the Trilinos implentation of the BiCGStab method through deal.II: TrilinosWrappers::SolverBiCGStab. One can specify the same parameters as in the GMRES solver with the exception of the maximum Krylov vectors parameters as a restart is not needed.
Some Remarks#
In the literature, it is stated that the reduction of error in one iteration of a bi-Lanczos-based method is approximately equal to that of two iterations in an Arnoldi-based method. Since the latter requires only one matrix-vector multiplication per iteration in comparison to the two required by the former, the costs of both methods are comparable in terms of the number of matrix-vector multiplications.
There is no general agreement on which method is better than the other as it heavily depends on the structure of the system matrix given by the specific problem. However, the disadvantages of each method can be overcome by using an effective preconditioner. There are several preconditioners in literature that have been proposed for these methods. In this document we will focus on the two classes of preconditioners used in Lethe, which are both in the category of “black box” preconditioners.
For more information about the parameters for these solvers and some practical tips visit the Linear Solver Control page in the CFD parameters section.
Preconditioners#
The preconditioners are essential to accelerate convergence of the iterative methods. Both the efficiency and the robustness of the Krylov subspace methods can be improved by using preconditioning. In general, a preconditioner \(\mathcal{M}\) (a matrix) transforms the original linear system \(\mathcal{A} x = b\) into another system which has the same solution, but better properties for the iterative solution (i.e., it is easier to solve). The main requirement to choose this new matrix is, from the practical point of view, that it is cheap to solve linear systems of these type \(\mathcal{M}x = b\). The preconditioner can be applied either from the left or to the right of the system matrix:
Left preconditioning solves the following linear system:
using the Krylov subspace \(\mathcal{K}(\mathcal{M}^{-1} \mathcal{A}, \mathcal{M}^{-1} b)\) instead of \(\mathcal{K}(A,b)\).
Right preconditioning solves the following linear system:
using the Krylov subspace \(\mathcal{K}(\mathcal{A} M^{-1}, b)\) and \(x = M^{-1} y\).
If the preconditioner is factorized, i.e., \(\mathcal{M} = \mathcal{M}_L \mathcal{M}_R\), where the matrices to the right-hand side are typically triangular matrices, the preconditioner is split as follows:
The advantages between the left-, right- and factorized preconditioning are highly dependent on the problem and whether the system matrix is symmetric or not. The convergence of the preconditioned Krylov subspace method is then determined by the eigenvalues of \(\mathcal{M}^{-1}\mathcal{A}\) (\(\mathcal{A}\mathcal{M}^{-1}\)) and their eigenvectors. When implementing Krylov subspace methods, it is necessary to generate a basis of the Krylov subspace \(\mathcal{K}_k (\mathcal{A}, r_0)\) and for reasons of numerical stability this basis should be orthogonal. When creating preconditioners, this condition needs to be considered, which leads to the following requirements for the two iterative methods used in Lethe:
GMRES: requires a symmetric preconditioner \(\mathcal{M}\)
BiCGStab: can use any general preconditioner \(\mathcal{M}\)
Incomplete LU Factorization (ILU)#
For detailed information on the theory, algorithms and implementation of the different ILU versions, we recommend the book Iterative Methods for Sparse Linear Systems by Y. Saad. This method performs an approximate factorization of the system matrix that consists of a sparse lower-diagonal matrix \(L\) and a sparse upper-diagonal matrix \(U\):
Then, in the preconditioned Krylov solver, \(\mathcal{M}^{-1} y\) is computed by calculating \(z = L^{-1} y\) and then \(U^{-1} z\). Let us also define the following matrix:
The different versions of ILU arise with the different treatment of the non-zero entries of the matrix \(\mathcal{R}\):
The first ILU version is known as the no-fill or ILU(0) method. It is defined as any pair of matrices \(L\) and \(U\) so that the elements of \(\mathcal{R}\) are zero in the same locations as the system matrix \(\mathcal{A}\). This version may be insufficient to yield an adequate rate of convergence for certain problems. In practice this leads to an algorithm where the pattern of \(LU\) is equal to the zero pattern of \(\mathcal{A}\).
Another version is the ILU(k) which allows for fill-in elements and improves accuracy. It consists of zeroing out all the fills of level \(k+1\) or higher. At this point it is important to define the level of fill which is used when calculating the entries of the preconditioner; it is basically a way of identifying which elements should be preserved or dropped in the final preconditioner. If the level of fill is set to \(k`\) then all the fill-in elements whose level does not exceed \(k\) are kept. In other words, this method constructs \(LU\) by setting the pattern of this matrix to be the pattern obtained in the decomposition of ILU(k-1). The drawbacks of this approach are that the amount of fill-in and computational work are not predictable for \(k > 0\). Also, when the fill is large, this is closer to the traditional LU decomposition, which is expensive. In addition, one must keep in mind that if there are elements that are being dropped and have a significant value, this will lead to an inaccurate factorization and to a larger number of iterations.
The previous versions can be further improved by using a threshold that controls the fill-in according to the size of the elements and not only based on the sparsity pattern. One version is known as ILU with dual thresholding (ILUT).
Other versions try to reduce the effect of dropping some entries by taking into account the values in other ways. For example, one strategy is to add all the elements that have been dropped and subtract them from the diagonal entry of \(U\). This is another version called Modified ILU (MILU) factorization.
In Lethe we use the Trilinos implementation of the ILU method through deal.II: TrilinosWrappers::PreconditionILU. This is of type ILU(k) (where ILU(0) is also included). To use this preconditioner, one must specify the absolute and relative tolerance, and the fill level. If the iterative method fails, the fill level is increased by 1 and the solution is attempted again. The absolute tolerance (\(\alpha \geq 0\)) and relative tolerance (\(\beta \geq 1\)) change the diagonal of the matrix before factorization as follows: a diagonal entry \(a_{ii}\) is replaced by \(\alpha \: \text{sign}(a_{ii}) + \beta a_{ii}\). This strategy can improve the conditioning of the matrix in some cases. Default values in Lethe are 1e-12 and 1, respectively.
Algebraic Multigrid (AMG)#
The main idea of multigrid methods is to efficiently correct all components of the error of the solution on the fine level by using a hierarchy of levels composed by coarser grids. A good starting point to understand all the details of AMG is given in A Multigrid Tutorial. In the algebraic multigrid method, the levels are generated based on the algebraic structure of the system matrix \(\mathcal{A}\). It requires the definition of three key components:
Smoother: it is in charge of smoothing the solution of the residual equations at each level. In general, one uses stationary iterative methods, such as Jacobi or Gauss-Seidel or incomplete LU. The former two are known for having difficulties as preconditioners for saddle-point problems, therefore, in Lethe we use ILU(k).
Intergrid operators: operators that allow to move between the different levels. They are also known as prolongation (\(P\)) and restriction (\(R\)) and are used to define the matrices of each level \(l\): \(\mathcal{A}^{l-1} = R \mathcal{A}^l P\).
Coarse grid solver: performs a correction of the solution in the coarser level of the hierarchy. In Lethe we use also an ILU decomposition.
Some remarks:
AMG is more expensive than ILU in terms of setup time and cost per iteration, but it accelerates convergence more significantly.
There are two types of AMG methods: the classical AMG and smoothed aggregation AMG, which differ in their coarsening strategy that in turn leads to differences also in the prolongation and restriction operators.
In Lethe we use the Trilinos implementation of the AMG method through deal.II: TrilinosWrappers::PreconditionAMG. One must specify several parameters related to the number of cycles, the type of cycle and smoother parameters.
References#
[1] M. Benzi, G. H. Golub, and J. Liesen, “Numerical solution of saddle point problems,” Acta Numer., vol. 14, pp. 1–137, May 2005, doi: 10.1017/S0962492904000212.
[2] A. Ghai, C. Lu, and X. Jiao, “A comparison of preconditioned Krylov subspace methods for large-scale nonsymmetric linear systems,” Numer. Linear Algebra Appl., vol. 26, no. 1, p. e2215, 2019, doi: 10.1002/nla.2215.
[3] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed. SIAM, 2003.