6.2.5 Newton-Krylov Method Newton-Krylov方法[cfd-6-2-5]

First of all, let us rewrite the implicit scheme given by Eq. (6.28) as

where \(\bar{J}\) represents the implicit operator (system matrix). As we already saw, \(\bar{J}\) constitutes a large, sparse, and generally non-symmetric matrix. In the previous subsections, we discussed two methods that decompose \(\bar{J}\) into several factors which can be each more easily inverted than \(\bar{J}\) itself. However, due to the factorisation error (and approximate linearisation of \(\vec{R}^{n+1}\)), only a linear convergence to steady state can be achieved. In order to obtain the quadratic convergence of Newton's method for the solution of non-linear equations, four conditions must be fulfilled:

• the linearisation of the residual must be exact,
• \(\bar{J}\) must be accurately inverted,
• the time step has to be \(\Delta t \to \infty\),
• initial solution must be, in some sense, close to the final solution.

Obviously, the main obstacles that have to be overcome are the linearisation and the inversion of the full system matrix.

A particularly suitable class of iterative techniques for the solution of large linear equation systems are the so-called Krylov-subspace methods. Several were proposed for the inversion of matrices which arise in CFD. Examples are the Conjugate Gradient Squared (CGS) method [76], the Bi-Conjugate Gradient Stabilised (Bi-CGSTAB) scheme [77], or the Transpose-Free Quasi-Minimum Residual (TFQMR) approach [78]. However, the most successful Krylov subspace method became the Generalised Minimal Residual (GMRES) technique, which was originally suggested by Saad and Schulz [79], [80]. Since then, the GMRES method was improved and augmented by several researchers [81]-[84]. Because of its popularity, we shall focus on the GMRES approach in the following. Nevertheless, most of what we shall discuss also applies to the other Krylov subspace methods.

As we mentioned in Subsection 3.2.2 (see also Appendix A.12), the GMRES method minimises the norm of the global residual, i.e., \(\|\bar{J}\Delta\vec{W}^{n} + \vec{R}^{n}\|\) over a set of \(m\) orthonormal vectors (search directions), which span the Krylov subspace \(\mathcal{K}_{m}\) given by Eq. (3.10). The GMRES algorithm can be summarised as follows:

• guess a starting solution \(\Delta\vec{W}^{n}_{0}\) and evaluate the initial residual vector \(\vec{r}_{0} = \bar{J}\Delta\vec{W}^{n}_{0} + \vec{R}^{n}\),
• generate the \(m\) search directions (by Gram-Schmidt orthogonalisation),
• solve the minimisation problem,
• form an approximate solution of Eq. (6.61) as \(\Delta\vec{W}^{n} = \Delta\vec{W}^{n}_{0} + \vec{y}_{m}\).

Since the memory requirements increase linearly with the number of search directions, \(m\) is restricted to values between 10 and 40 in practice. This might not be sufficient for a converged solution \(\Delta\vec{W}^{n}\). Thus, the GMRES method has to be restarted, i.e., we set \(\Delta\vec{W}^{n}_{0} = \Delta\vec{W}^{n}\), compute \(\vec{r}_{0}\) and proceed with step 2. As pointed out in Ref. [85], instead of working with a constant number of search directions, \(m\) should be reduced if the norm of the global residual drops below a specified tolerance. In this way, a large number of operations can be saved in later stages of the Newton iteration.

GMRES and other Krylov subspace methods allow us to circumvent an explicit computation and storage of the flux Jacobian \(\partial\vec{R}/\partial\vec{W}\). The idea is based on the observation that the methods rely only on matrix-vector products of the form \(\bar{J}\Delta\vec{W}^{n}\) and do not need the matrix \(\bar{J}\) explicitly. The product of the flux Jacobian with the solution update can be approximated by a simple finite difference as

which requires only two evaluations of the residual. The stepsize \(h\) has to be chosen with some care, in order to minimise the numerical error [36]. One particularly suitable formulation reads [86]

where \(\epsilon\) denotes the machine accuracy, \(d\) the scalar product \(\vec{W}^{n}\cdot\Delta\vec{W}^{n}\), \(|\Delta\vec{W}^{n}|\) is the vector \(\Delta\vec{W}^{n}\) with all elements set to their absolute values, and finally \(\mathrm{typ}\,U\) represents a typical size of \(U\). Apart from saving memory and operations, there is an even more important advantage of the finite-difference approximation. Namely, numerically accurate linearisation of a high-order residual \(\vec{R}^{n}\) (including boundary conditions, limiters, source terms, etc.) can be easily achieved.

Thus, the quadratic convergence of Newton's scheme can be realised at moderate costs. For this reason, we speak of such a scheme as of Newton-Krylov approach [86]-[88].

The efficiency of Krylov-subspace methods depends strongly on a good preconditioner. Its purpose is to cluster the eigenvalues of the system matrix \(\bar{J}\) around unity. Thus, instead of Equation (6.61), the left- or right-preconditioned system according to Eq. (3.11) is solved. Using Eqs. (6.61) and (6.62) together with the condition \(\Delta t \to \infty\), the Newton-Krylov method becomes

with left preconditioning, and

in the case of right preconditioning, respectively. The main difference between the two preconditioning methodologies is that left preconditioning scales the residual \(\vec{R}^{n}\) whereas right preconditioning does not. This has to be kept in mind when the convergence of the Krylov method is monitored.

Obviously, the preconditioner should be as close as possible to the inverse of the system matrix (\(\bar{P}_{L,R} \approx \bar{J}^{-1}\)). But on the other hand, it should be invertible with low numerical effort. Therefore, we have to find an optimal tradeoff between the convergence speed of the Krylov method and the time spend for inverting the preconditioning matrix. One of the most successful preconditioners is the Incomplete Lower Upper factorisation method [89], [90] with varying level of fill-in (mostly with zero, designated as ILU(0)). The ILU preconditioner is especially efficient in the case of viscous, turbulent flows, i.e., for stiff equations [91]. In order to obtain a good performance on unstructured grids, it is necessary to reorder the elements of the system matrix such that the bandwidth is reduced. We already discussed the RCM renumbering strategy [27], [28] in Section 6.2.1.

A serious disadvantage of the ILU preconditioning scheme is that elements of the matrix \(\bar{J}\) have to be computed (see Subsection 6.2.2) and stored. Therefore, some authors suggested to employ the LU-SGS scheme as a preconditioner [92], [93]. Hence, when Eq. (6.62) is applied, the formation and storage of \(\bar{J}\) is completely avoided. However, it was demonstrated in Ref. [91] on behalf of several 2-D cases that the GMRES method with LU-SGS is inferior to GMRES combined with ILU(0) in terms of the CPU-time. Nevertheless, the LU-SGS scheme (best when coupled with multigrid) still represents an attractive alternative, particularly in 3D.

The time step of the implicit Newton-Krylov method is infinitely large. However, it is advisable to use small time steps at the beginning of the Newton iteration process. The reason is that the flow solution is in general far from the steady state at the beginning of the solution process, i.e., the root of the nonlinear equation

and this may cause a breakdown of the Newton iteration. One possible remedy is the so-called Switched Evolution Relaxation (SER) technique [94]. Here, the term \(\Omega/\Delta t\) is retained in Eq. (6.61). The time step is evaluated in the same way as presented for the explicit scheme (Eq. (6.14) or Eq. (6.20) without the viscous eigenvalue). The CFL number \(\sigma\) is increased starting from a small initial value correspondingly to the reduction of the 2-norm of the residual, i.e.,

Hence, the convergence of the iteration procedure (6.61) will be at first linear, but it approaches the quadratic convergence of Newton's method for large CFL numbers. A further effect of the time term is the increased diagonal dominance of \(\bar{J}\) (inverse proportional to \(\Omega/\Delta t\)), which will help to stabilise the iteration. In Ref. [88], it was suggested to clip \(\sigma^{n+1}\) in Eq. (6.66) such that it increases by maximum factor of two and decreases less than factor of ten.

A further approach which can be used to overcome the start-up problems of Newton's method consist of grid sequencing, where the initial solution is obtained on a sequence of coarser grids and interpolated onto finer grids. Another possibility is to use a numerically cheap but robust iteration scheme for the initial guess. For example, we could start with a multigrid scheme driven by the LU-SGS method and then switch to GMRES with LU-SGS as preconditioner. This might be particularly interesting for viscous turbulent flows, where the convergence of a multigrid scheme usually slows down after the initial phase. However, the global flow solution is then already close to the steady state.