kinKadath numerical method strategy
Newton iteration level
are the loop over the Newton corrections and the loop building up the Krylov subspace out
To solve this system of equations, a general Newton iteration scheme is used :
- Set \(u_0\) an initial guess
- For k = 0, 1, 2, …until convergence do:
- Solve J(\(u_k\)) \(\delta u_k\) = -F(\(u_k\))
- Set \(u_{k+1} = u_k + \lambda \, \delta u_k\), \(0 < \lambda \leq 1\)
- Tests for convergence
Where \(u_k\) is the k\(^{th}\) iterate to \(u\), and \(J(u) = F'(u)\) is the Jacobian system. A Newton correction \(\lambda \, \delta u_k\) is add to the previous value \(u_k\) to produce a new iterate \(u_{k+1}\). The overall loop stop when a convergence criteria is reached, for example when we obtain a required drop in the norm of the nonlinear residual or a sufficiently small Newton update.
At each stage of the iteration process, a linear system is solve in 2(a), based on the Jacobian system. The basic method implemented in Kadath use a direct method based on a Gauss algorithm with a \(LU\) factorization to solve the Jacobian at each step. Even if this method is robust, this is costly in memory, especially for very large systems.
Because Kadath use a spectral method to solve complex physical problems the system of equations \(F(x)=0\), the provided Jacobian matrix is generally a non-symmetric and dense matrix. Even for this worth case, iterative Krylov methods, combined with a newton loop can provide an optimize way to solve the given system of equations. In addition, this is strongly motivated by the idea that some improvement could be made to structure the coefficients of the spectral method, and thus better exploit the defined domains splitting to improve the numerical efficiency.
General Krylov methods
Krylov methods are iterative method based on Krylov subspace projection methods to solve a non singular large system of equations \(Ax=b\). Each Krylov subspace \(K_j\) is spanned by a iterative generation of vectors and defined by
where \(r_0=b-Ax_0\). here, \(A = J\), \(r_0 = -F(x_0) - J \delta u_0\). A sequence of approximate solutions \(x_m\) is seek from an affine subspace \(x_0 + K_m\), so that the corresponding residual \(r_m\) converge to zero in a finite number of iterations.
Different methods is used to construct this subspace basis. For example, for non-symmetric problems, Generalized Minimal RESidual method (GMRES) use an Arnoldi process and the Bi-Conjugate Gradient STABilized (BiCGSTAB) method use a Lanczos process [1], [4]. They are the main preferential method chosen by Kinkadath for our typical problem.
This gives a linear combination of the defined orthonormal basis vectors, as
Newton-Krylov Method
The Krylov method gives an approximation of the \(\delta u_0\) solution within a subspace of small dimension, but because the Newton correction \(\delta u_0\) drop to zero in the late Newton iteration, this gives a asymptotically reasonable guess in the Newton loop (starting with zero).
What we can observe to this, is that there is few need of Jacobian vector product and no explicit Jacobian computing is needed. We can also simplicity approximate this Jacobian vector product using a Fréchet derivatives. There is also others useful techniques as the automatic derivative. [54], [27].
Because we don’t need an exact solution, we can speed up the resolution by exploiting the iterative form of the Krylov method and approximately solve the Newton correction. We generally use a criteria as
holds. \(\eta_c\) is called the forcing term. We generally call this kind of method a Inexact Newton method.
The global convergence properties of theses algorithms are improved by additional technique, as the linesearch backtracking strategy [2], [3].
Two of the more popular methods for globalization of Newtons method are the line search method and the trust region method [54,84].
Preconditioning technique
preconditioning technique
The Newton Krylov methods gives there best efficiency when using preconditioning and a proper Jacobian information update strategy technique. We are using a right preconditioning to avoid changes in the norm of the residual of the convergence critera, but this is also possible to use a left preconditioning. When this is applied to the linear systems, we obtains the system
to solve, where \(P\) is the right preconditioning matrix. This can be solved in two step by adding the Preconditioning solve process. \(P\) can be solved approximately by a linear iterative method and only the inverse of the preconditioning by a vector operator need to be keep.
In practical case, the preconditioning \(P\) should be a linearized approximation of the system Jacobian matrix \(J\). We generally use a simplified version of the physical problem (physic based preconditionning). Combined this with a approximated linear iterative method gives the best performances. Additive Schwarz method gives one of the best performance when we can split the discretization in domains, that translate in subblocks in the Precondionning matrix.
Another performant approache, specially for parallel resolution is the Multigrid approaches. linear multigrid algorithm gives good performance to estimate a precondionning matrix [29,141].
the kinKadath preconditioning
An effective preconditioner is required for each iteration and determining such may be hard or expensive.
In pratique, because of the dense structure of the spectral method, this is very complicated to highlight a structure that can be exploited as a preconditionner. The best approache would be to reorganize the vector structures used in the spectral Kadath software and better exploit the domains. This can be also very dependant to the problem to solve.
The approach in kinKdath was to give a general precondionning that could be used by default in any Kadath problem, giving the opportunity to add better Precondionning technique dependant of each problem to solve.
The preconditioning matrix \(P\) use a pre-computation phase, called setup step, which store the pivots elements. The second step is to solve the Precondionning matrix using theses precomputed pivots elements. Thus, we can optimize in the Kinsol library the call of each step (see kinsol). To compute the preconditioning matrix \(P\), we are using a direct linear solvers, applied on the problem, and we will optimize the call of the two step described above. We will keep the privots and/or the computation of the Preconditioning inverse as long as there is still few differences. A block-diagonal approach can be used, but this will need a better domain partition to represent almost independent blocks.
first step is using a standard Gaussian elimination with partial pivoting to perform a LU factorization of the Jacobian matrix. The pivot information and this \(LU\) factorization are then used to solve the Jacobian system.
Optimize call strategy
The preconditioner matrix P in the case of iterative linear solvers) as infrequently as possible to balance the high costs of matrix operations against other costs. different tuning parameters that estimate if the Preconditionneur change a lot or not (and fore fixed number of iter).
for example; specifies the maximum number of nonlinear iterations that can be performed between calls to the preconditioner or Jacobian setup function. This is set to T = 10 by default. If we define default value of mbset, the number of nonlinear iterations after which a Jacobian information update is enforced, we can obtent the same newton algorithm.