'Alternating direction implicit' preconditioner

[johnomotani]

The preconditioner within the Jacobian-free Newton Krylov (JFNK) solver used for implicit timestepping needs to 'invert a matrix'. The initial implementation does this by an LU-factorisation method, using Julia's LinearAlgebra package. This LU solver is serial. Some parallelisation can be achieved by constructing separate preconditioners in each distributed-MPI block, at the cost of dropping the elements of the preconditioner that would couple different blocks, but our shared-memory parallelism cannot be exploited by the LU solver, creating a bottleneck where for a 1D1V test with 32 elements in the z-direction, no speed up can be gained by using 128 cores on ARCHER2 compared to 8 cores.

An alternative which is less accurate but more parallelisable is based on the 'alternating direction implicit' (ADI) approach. The aim is to solve

$$J.\delta x = R$$

for $\delta x$.
In the ADI approach we:

* First identify all the terms that only couple in $w_\parallel$. Split the Jacobian matrix into an 'implicit' matrix $J_{I,w}$ containing the terms that couple in $w_\parallel$ - which is block diagonal because grid points at different $z$ values are uncoupled - and an 'explicit' matrix $J_{E,w}$ with the remainder

$$J = J_{I,w} + J_{E,w}$$

Take an inital guess $\delta x_0 = \vec{0}$ and find an updated guess $\delta\tilde{x}_1$ from

$$J_{I,w}.\delta\tilde{x}_1 = R - J_{E,w}.\delta x_0 = R$$

* Second identify all the terms that only couple in $z$, and do a similar thing

$$J = J_{I,z} + J_{E,w}$$

where $J_{I,z}$ is block diagonal, this time points at different $w_\parallel$ are uncoupled (also $g_e$ and $p_{e\parallel}$ are uncoupled). Solve for a new guess $\delta x_1$ from

$$J_{I,z}.\delta x_1 = R - J_{E,z}.\delta \tilde{x}_1$$

* Optionally iterate these two steps to get $\delta x_2$, $\delta x_3$, etc. which should be better approximations to the solution of $J.\delta x = R$. Initial tests suggest that a single iteration might be enough - using more iterations in the ADI solve would slightly reduce the number of JFNK iterations, but not by the factor of 2,3,4... that would be needed to cover the cost of the extra ADI iterations.

Todo:

• Add a test that checks that an ADI solver converges to the LU solution after enough iterations?

[johnomotani]

This PR is a draft for the moment, because I want to test whether it is possible and useful to add communication between different distributed-MPI blocks after each ADI iteration - could doing that and taking 2 or 3 ADI iterations speed up the JFNK convergence?

[johnomotani]

Initial performance testing for a 1D1V kinetic electron simulation on ARCHER2 with 32 elements in the $z$-direction split into 8 blocks. The grid in $z$ is uniform. Times and iteration counts are taken over 5 output steps, running from 0.004 to 0.009 $L_{ref}/c_{s,ref}$. This is still in an initial transient phase, which is maybe not ideal, but is as far as the 8-core runs (or 128 core LU) would get in 4 hours on 1 node of ARCHER2.

	Run time (s)	Total Krylov iterations	Time per preconditioner evaluation (ms)
LU, 8 cores	57.4	600093	4.60
LU, 128 cores	58.2	586143	4.24
ADI, 8 cores	115.8	885619	6.19
ADI, 128 cores	23.4	913841	0.397

This seems to agree with what we would expect. The ADI solver takes a few more iterations because it is only an approximate matrix inverse, whereas the LU is exact (within each block), but the LU solve does not scale with the number of shared-memory cores, whereas the ADI solve does - with almost perfect scaling as 6.19/0.397=15.6 times speed-up with 16 times more cores.

Suspect the run time for the 8 core ADI may be slightly high due to random fluctuation as based on the iteration count we would expect more like 90 s. This might make the perfect scaling mentioned above over-optimistic, depending whether the random slowdown happened (if it happened at all) during a call to the preconditioner or while the code was somewhere else.