Implement Poisson solver based on integrated Green functions

[aeriforme]

This uses the method from:
https://journals.aps.org/prab/abstract/10.1103/PhysRevSTAB.9.044204
https://journals.aps.org/prab/pdf/10.1103/PhysRevSTAB.10.129901

Because this uses the same FFT infrastructure as PSATD, right now, the code needs to be compiled with:
-DWarpX_PSATD=ON.

Also: since the solver uses a global FFT, right now only one MPI rank participates in the FFT solve (but this MPI rank can still use a full GPU) ; this could be improved later by using heFFTe. The rest of the PIC operations (charge deposition, gather, etc.) are still performed by the different MPI ranks in parallel.

Comparison with Basseti-Erskine formula (see the test script in this PR):

TODO

• Test on GPU
• Add back the relativistic factor in dx
• Overcome limitation with single MPI rank (use ParallelCopy)
• Fix bug in the multi-GPU implementation, in the case grid_type = collocated (this could be done with the regular Poisson solver instead of the relativistic one, so that we can output phi)
• Add timers to the FFT and the calculation of the IGF
• Finalize the automated test
• Decide on the user interface: introduce a new parameter warpx.poisson_solver that can be set to multigrid, fft_based, etc. + how to handle compatibility with boundary conditions?
  if electromagnetic simulation with PML + initialize_self_fields, then use IGF which implies open boundaries
  IGF only compatible with open and viceversa
  multigrid not compatible with open
• This should follow the same logic as the Maxwell solver. In the PICMI interface, there is already a “method” parameter for the Poisson solver.

more

• add new boundary condition “open” that is only valid for Poisson-like field solvers + organize documentation to group boundaries in similar types

more

• Abort if someone tries to use mesh refinement If mesh refinement is on and the user has selected the IGF solver, then use IGF only on the coarsest level (lev = 0) and MLMG in the refined patches
• Properly set guard cells of phi outside of the global domain (useful for grid_type = collocated)
• Save the FFT plans in a permanent object, so that we do not need to recompute them. Probably not needed. From the timers, it seems that the creation and destruction of FFT plan does not take any significant time.
• Clean-up ; move code to a different file, maybe IntegratedGreenFunctionSolver.H, in ablastr
• Handle dimensionality that are not 3D, abort if we are not in 3D
• Address warning

/Users/rlehe/Documents/codes/warpx_directory/warpx/Source/ablastr/fields/PoissonSolver.H:69:1: warning: function 'computePhiIGF' is not needed and will not be emitted [-Wunneeded-internal-declaration]
computePhiIGF ( amrex::MultiFab const & rho,

By calling the function in VectorPoissonSolver.H (for the magnetostatic solver)

• Fix ParallelCopy error that appears in DEBUG mode:

Assertion `boxArray().ixType() == src.boxArray().ixType()' failed, file "/global/u2/a/aforment/src/war

• Documentation

[ax3l]

Awesome 🚀 ✨

[RemiLehe]

@aeriforme Are you able to reproduce the build the Sphinx documentation on your local computer? I wonder why it fails here and not on other PRs.

[ax3l]

The openbc_poisson_solver test fails with a checksum issue right now.

Are you able to reproduce the build the Sphinx documentation on your local computer? I wonder why it fails here and not on other PRs.

Found it 🎉 Just a small typo, the bib file has invalid syntax. Will fix.

[RemiLehe]

I updated the checksum. It was a negligible change in Bz. (Bz is supposed to be almost 0 for a relativistic beam propagating along z.)

[ax3l]

Small update I would suggest:

• If FFTs were not compiled in, but is_solver_multigrid=false is passed to this function, then this silently falls back to a MLGM solve for the lowest level. This should raise a runtime error, because we cannot fulfill the API/user request.
  • An alternative, throwing a warning and "loud" fallback to MLMG, is not ideal in practice, because boundaries/padding & resolution would be quite different.
• I realized the is_solver_multigrid name is a but confusing, because we do the IGF only on level 0 and an MLMG on all others in MR. Maybe we find a good alternative name, e.g., is_solver_igf_on_lev0 (most explicit/good?) or is_mlmg_all_levels (kinda like now).

[aeriforme]

The first point does not really apply to WarpX, I believe.
There's an assert that checks if FFTs are compiled when the IGF solver is selected.
It's here.
Anyway, I added a similar assert in ablastr in the follow-up PR #4945.
Let me know if this solution is ok for you.

[ax3l]

Thanks! Yes, we need the check in ABLASTR directly because that is a library that we also call from ImpactX 🙏