Improve spectral solver on staggered grids (#1468)

* Implement Galilean PSATD equations on staggered grids

* Implement high-order interpolation in 2D/3D

* Include missing header file and small clean-up

* Fix bug for FDTD build

* Small clean-up

* Modify current correction for staggered grids

* Implement comoving PSATD scheme (formulation with rho)

* Fix single-precision builds

* Do not implement rho-free formulation for comoving PSATD yet

* Invert sign of comoving velocity to match Galilean convention

* Fix two bugs in comoving PSATD algorithm

* Update benchmark of CI test momentum-conserving-gather

* Update benchmark of CI test PlasmaAccelerationMR

* Update documentation

* Clean up comoving PSATD class

* Clean up comoving PSATD class (more)

* Clean up comoving PSATD class (more)

* Implement changes requested in PR review

* Add 2D regression test for staggered Galilean PSATD

* Add 2D regression test for staggered comoving PSATD

* Unify input files for new CI tests to avoid duplication

* Fully rebase benchmarks changes on development

* Update benchmark of Galilean hybrid test after #1536

1 files changed, 14 insertions, 6 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index 0cfd899df..4c9cb6967 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -10,6 +10,7 @@
 #include "SpectralAlgorithms/GalileanAlgorithm.H"
 #include "SpectralAlgorithms/AvgGalileanAlgorithm.H"
 #include "SpectralAlgorithms/PMLPsatdAlgorithm.H"
+#include "SpectralAlgorithms/ComovingPsatdAlgorithm.H"
 #include "WarpX.H"
 #include "Utils/WarpXProfilerWrapper.H"
 #include "Utils/WarpXUtil.H"
@@ -39,6 +40,7 @@ SpectralSolver::SpectralSolver(
                 const int norder_x, const int norder_y,
                 const int norder_z, const bool nodal,
                 const amrex::Array<amrex::Real,3>& v_galilean,
+                const amrex::Array<amrex::Real,3>& v_comoving,
                 const amrex::RealVect dx, const amrex::Real dt,
                 const bool pml, const bool periodic_single_box,
                 const bool update_with_rho,
@@ -64,15 +66,21 @@ SpectralSolver::SpectralSolver(
                 k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt);
         }
         else {
-            if ((v_galilean[0]==0) && (v_galilean[1]==0) && (v_galilean[2]==0)){
-                // v_galilean is 0: use standard PSATD algorithm
-                algorithm = std::make_unique<PsatdAlgorithm>(
-                   k_space, dm, norder_x, norder_y, norder_z, nodal, dt, update_with_rho);
-            }
-            else {
+            // Galilean PSATD algorithm
+            if (v_galilean[0] != 0. || v_galilean[1] != 0. || v_galilean[2] != 0.) {
                 algorithm = std::make_unique<GalileanAlgorithm>(
                     k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt, update_with_rho);
             }
+            // Comoving PSATD algorithm
+            else if (v_comoving[0] != 0. || v_comoving[1] != 0. || v_comoving[2] != 0.) {
+                algorithm = std::make_unique<ComovingPsatdAlgorithm>(
+                    k_space, dm, norder_x, norder_y, norder_z, nodal, v_comoving, dt, update_with_rho);
+            }
+            // Standard PSATD algorithm
+            else {
+                algorithm = std::make_unique<PsatdAlgorithm>(
+                    k_space, dm, norder_x, norder_y, norder_z, nodal, dt, update_with_rho);
+            }
         }
     }