nonlinear_solvers

Nonlinear solvers, using Jacobian-free Newton-Krylov methods.

These solvers use an outer Newton iteration. Each step of the Newton iteration requires a linear solve of the Jacobian. An 'inexact Jacobian' method is used, and the GMRES method (GMRES is a type of Krylov solver) is used to (approximately) solve the (approximate) linear system.

parallel_map() is used to apply elementwise functions to arbitrary numbers of arguments using shared-memory parallelism. We do this rather than writing the loops out explicitly so that newton_solve!() and linear_solve!() can work for arrays with any combination of dimensions.

Useful references: [1] V.A. Mousseau and D.A. Knoll, "Fully Implicit Kinetic Solution of Collisional Plasmas", Journal of Computational Physics 136, 308–323 (1997), https://doi.org/10.1006/jcph.1997.5736. [2] V.A. Mousseau, "Fully Implicit Kinetic Modelling of Collisional Plasmas", PhD thesis, Idaho National Engineering Laboratory (1996), https://inis.iaea.org/collection/NCLCollectionStore/Public/27/067/27067141.pdf. [3] https://en.wikipedia.org/wiki/Generalizedminimalresidualmethod [4] https://www.rikvoorhaar.com/blog/gmres [5] E. Carson , J. Liesen, Z. Strakoš, "Towards understanding CG and GMRES through examples", Linear Algebra and its Applications 692, 241–291 (2024), https://doi.org/10.1016/j.laa.2024.04.003. [6] Q. Zou, "GMRES algorithms over 35 years", Applied Mathematics and Computation 445, 127869 (2023), https://doi.org/10.1016/j.amc.2023.127869

gather_nonlinear_solver_counters!(nl_solver_params)

Where necessary, gather the iteration counters for the nonlinear solvers.

Where each solve runs in parallel using all processes, this is unnecessary as the count on each process already represents the global count. Where each solve uses only a subset of processes, the counters from different solves need to be added together to get the global total.

Apply the GMRES algorithm to solve the 'linear problem' J.δx^n = R(x^n), which is needed at each step of the outer Newton iteration (in newton_solve!()).

Uses Givens rotations to reduce the upper Hessenberg matrix to an upper triangular form, which allows conveniently finding the residual at each step, and computing the final solution, without calculating a least-squares minimisation at each step. See 'algorithm 2 MGS-GMRES' in Zou (2023) [https://doi.org/10.1016/j.amc.2023.127869].

newton_solve!(x, rhs_func!, residual, delta_x, rhs_delta, w, nl_solver_params;
              left_preconditioner=nothing, right_preconditioner=nothing, coords)

x is the initial guess at the solution, and is overwritten by the result of the Newton solve.

rhs_func!(residual, x) is the function we are trying to find a solution of. It calculates

\[\mathtt{residual} = F(\mathtt{x})\]

where we are trying to solve $F(x)=0$.

residual, delta_x, rhs_delta and w are buffer arrays, with the same size as x, used internally.

left_preconditioner or right_preconditioner apply preconditioning. They should be passed a function that solves $P.x = b$ where $P$ is the preconditioner matrix, $b$ is given by the values passed to the function as the argument, and the result $x$ is returned by overwriting the argument.

coords is a NamedTuple containing the coordinate structs corresponding to each dimension in x.

Tolerances

Note that the meaning of the relative tolerance rtol and absolute tolerance atol is very different for the outer Newton iteration and the inner GMRES iteration.

For the outer Newton iteration the residual $R(x^n)$ measures the departure of the system from the solution (at each grid point). Its size can be compared to the size of the solution x, so it makes sense to define an `error norm' for $R(x^n)$ as

\[E(x^n) = \left\lVert \frac{R(x^n)}{\mathtt{rtol} x^n \mathtt{atol}} \right\rVert_2\]

where $\left\lVert \cdot \right\rVert$ is the 'L2 norm' (square-root of sum of squares). We can further try to define a grid-size independent error norm by dividing out the number of grid points to get a root-mean-square (RMS) error rather than an L2 norm.

\[E_{\mathrm{RMS}}(x^n) = \sqrt{ \frac{1}{N} \sum_i \frac{R(x^n)_i}{\mathtt{rtol} x^n_i \mathtt{atol}} }\]

where $N$ is the total number of grid points.

In contrast, GMRES is constructed to minimise the L2 norm of $r_k = b - A\cdot x_k$ where GMRES is solving the linear system $A\cdot x = b$, $x_k$ is the approximation to the solution $x$ at the $k$'th iteration and $r_k$ is the residual at the $k$'th iteration. There is no flexibility to measure error relative to $x$ in any sense. For GMRES, a `relative tolerance' is relative to the residual of the right-hand-side $b$, which is the first iterate $x_0$ (when no initial guess is given). [Where a non-zero initial guess is given it might be better to use a different stopping criterion, see Carson et al. section 3.8.]. The stopping criterion for the GMRES iteration is therefore

\left\lVert r_k \right\rVert < \max(\mathtt{linear\_rtol} \left\lVert r_0 \right\rVert, \mathtt{linear\_atol}) = \max(\mathtt{linear\_rtol} \left\lVert b \right\rVert, \mathtt{linear\_atol})

As the GMRES solve is only used to get the right direction' for the next Newton step, it is not necessary to have a very tightlinear_rtol` for the GMRES solve.

reset_nonlinear_per_stage_counters!(nl_solver_params::Union{nl_solver_info,Nothing})

Reset the counters that hold per-step totals or maximums in nl_solver_params.

Also increment nl_solver_params.stage_counter[].

coords is a NamedTuple of coordinates corresponding to the dimensions of the variable that will be solved. The entries in coords should be ordered the same as the memory layout of the variable to be solved (i.e. fastest-varying first).

The nonlinear solver will be called inside a loop over outer_coords, so we might need for example a preconditioner object for each point in that outer loop.