Add a factored-Laplacian MLMG preconditioner for the Darwin GMRES solve

[roelof-groenewald]

On its solenoidal subspace the semi-implicit Darwin field operator nabla^4(Z) + curl(chi curl(Z)) reduces, in the constant-susceptibility limit, to (-lap)(-lap + chi). Preconditioning each Z component with two successive cell-centered MLMG solves - a Poisson solve then a Helmholtz solve using the local chi(x) built from the diagonal mass-matrix row-sums therefore captures the bulk of the operator's spectrum.

Measured on the 2D ES-coupled EM-modes CI case (Nz=128, 50 steps): GMRES drops from a mean of 736 iterations per solve (min 382, max 1834, with ~25 restart cycles each) to 7.7 (min 7, max 8, a single cycle), the GMRES share of the run from 10.6 s to 1.7 s, and total wall time from 16.9 s to 7.7 s. Both variants converge the same relative residual (~5e-5), and the resulting fields agree to 6e-5 in B and 5e-6 in the particle quantities.

The B-staggered components are collocated onto the cell-centered grid by index identification rather than averaging: pair-averaging annihilates each component's nodal-dimension Nyquist planes, leaving GMRES unable to reduce the deposited-current noise there (measured as a hard stall at a few 1e-2 relative residual). The index shift is spectrally exact in periodic dimensions, and any residual approximation only affects preconditioner quality, not correctness, since amrex::GMRES applies the preconditioner on the right and always converges the true residual.

Curl-free content is over-damped by (1 + chi/k^2), because the true operator reduces to nabla^4 alone there; a grad-div (Coulomb-gauge penalty) term in the operator would make the factorization exact on that subspace as well.

Enabled via amrex_gmres.pc_type = pc_darwin_mlmg (PICMI: GMRESLinearSolver(pc_type=DarwinMLMGPreconditioner())); off by default. The 2D Darwin EM-modes CI test now exercises it, so its checksum needs regenerating; the 1D test keeps covering the unpreconditioned path.

[roelof-groenewald]

We should decide whether to change the default Darwin behavior to use the preconditioner. I think we should since it is (so far) always a net win in performance without any observed loss in accuracy.

[RemiLehe]

Thanks for this PR!

Regarding setting the number of MLMG iterations: could you make this PR consistent with: #7199

[RemiLehe]

Is there a reason to choose MLABecLaplacian (cell-centered) rather than MLEBNodeFDLaplacian (node-centered). Note that MLEBNodeFDLaplacian is used in the rest of WarpX, including for the magnetostatic solver, which has a very similar issue (i.e. similar to the fact that Z is B-staggered, A in the magnetostatic solver is E-staggered ; however in the case of the magnetostatic solver, we use the averaging strategy which is discouraged here ; it would be nice to have a consistent strategy for these different solvers)

[roelof-groenewald]

It is mostly due to the fact that the cell-centered grid is contained within the staggered E-grid, so the cell label identification can be done - for a nodal grid extra grid points are needed and therefore require some interpolation or extension. When / if we switch to a nodal mass matrix we can switch this to also use a nodal MLMG solver.

[roelof-groenewald]

Regarding setting the number of MLMG iterations: could you make this PR consistent with: #7199

Thanks for flagging this @RemiLehe. Now updated.

[RemiLehe]

In this test, could you add an assert that checks the average number of GMRes iterations that were needed to converge, during the simulation.

This would guard against any future bug or regression of the preconditioner. (This type of bug might otherwise be undetectable, since the outer GMRes algorithm would make sure that we get the correct physics nonetheless. But an incorrect preconditioner might result in many more iterations to get there.)

[roelof-groenewald]

Good idea! We would need to read the stdout log to get that information, right? Do you know of a different test that currently does that so I can match the way it is done?

[RemiLehe]

Hm, now that I look into it, I don't think that there is an exact precedent.
It seems that https://github.com/BLAST-WarpX/warpx/blob/development/Examples/Tests/implicit/analysis_petsc_matrix.py#L18 does something similar, but relies on this diagnostic machinery:
https://github.com/BLAST-WarpX/warpx/blob/development/Source/NonlinearSolvers/NewtonSolver.H#L243
which only exists for non-linear solvers.

Maybe this should to be addressed in a follow-up PR.

[RemiLehe]

I guess this assert should be equivalent to "Is this driving the GMRes algorithm through a non-linear solver". However, the corresponding condition here (pc_type.name is not None) does not make this very clear. Could you explain this condition?

[roelof-groenewald]

Good point. I added a new class attribute for PreconditionerBase called supports_direct_gmres which is only True for the DarwinMLMGPreconditioner(). The assert then becomes the much clearer:

            assert pc_type.supports_direct_gmres, (
                f"{type(pc_type).__name__} cannot be selected directly on the "
                "GMRES solver; pass it to the nonlinear solver instead"
            )

This change also made me realize that we can collect preconditioner_type_initialize_inputs() into the parent class which reduced code duplication.

[RemiLehe]

I started adding a few comments. Would you have time to take a look.

[RemiLehe]

My understanding is that, starting from the full mass matrix $$\partial \boldsymbol{J}/\partial \boldsymbol{E}$$ this picks only the coupling $$\partial J_x/\partial E_x$$, $$\partial J_y/\partial E_y$$, $$\partial J_z/\partial E_z$$ and neglects the other couplings (e.g. $$\partial J_y/\partial E_x$$).
However, these couplings are probably not negligible for strongly-magnetized species.

If the above is correct, it might be worth adding a comment here, as well as in the documentation in parameters.rst.

Also: we could potentially lift this approximation in the future (in a follow-up PR).

[roelof-groenewald]

It is even a stronger reduction that that. We can think of the full mass matrix as a 3x3 block tensor, each block (9 in total) is a matrix holding the component couplings, i.e., M_{xx}, M_{xy}, M_{xz}, etc. Each M_{ij} is a matrix of (roughly) size N_cells x N_cells (the "roughly" is due to differences in grid staggering for the different components).

The reduction here is to take that set of matrices and reduce it to a single matrix. This is done by

• discarding the cross-component couplings, i.e., throw away M_{xy}, M_{xz}, M{yx}, etc. This leaves only M_{xx}, M_{yy}, and M_{zz}.
• Sum the rows in each M_{ii}. This leaves a three numbers per cell / grid point from each M_{ii} matrix.
• Average together those three numbers to get a single scalar value per cell center to use as the sigma for the Helmholtz solve.

Relaxing this could be done in a number of stages. For example, one level higher would be to not average the row summed M_{ii} values, but instead use different sigma MFs for the different components in the Helmholtz solve.

[roelof-groenewald]

The first stage suggestion above has been implemented in #7239.

[RemiLehe]

I don't think I understand why it is needed for the Multifab to be nodal in at most one dimension. Could you explain?

[roelof-groenewald]

It isn't strictly necessary since the node -> cell association would work even for fully nodal MFs. With the current implementation only supporting periodic boundaries, I don't think there would be any issue passing a MF with arbitrary staggering to the preconditioner. In the future, when we add support for non-periodic boundaries, I expect we'll need to be more careful about how the dropped plane in nodal directions are repopulated. At that point we will specialize the implementation for B-staggered MFs. There is already an assert later in this file that the a_x and b_x have B-field staggering (similar to the asserts in the other PC classes).
I'll remove this comment to avoid confusion in the future.