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

8 files changed, 801 insertions, 80 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
index a050bd5c2..ab8e439ed 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
@@ -5,6 +5,7 @@ target_sources(WarpX
     PsatdAlgorithm.cpp
     SpectralBaseAlgorithm.cpp
     AvgGalileanAlgorithm.cpp
+    ComovingPsatdAlgorithm.cpp
 )
 
 if(WarpX_DIMS STREQUAL RZ)
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H
new file mode 100644
index 000000000..2d9962e4d
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H
@@ -0,0 +1,92 @@
+#ifndef WARPX_COMOVING_PSATD_ALGORITHM_H_
+#define WARPX_COMOVING_PSATD_ALGORITHM_H_
+
+#include "SpectralBaseAlgorithm.H"
+
+#if WARPX_USE_PSATD
+
+/* \brief Class that updates the field in spectral space and stores the coefficients
+ * of the corresponding update equation, according to the comoving spectral scheme.
+ */
+class ComovingPsatdAlgorithm : public SpectralBaseAlgorithm
+{
+    public:
+
+        /**
+         * \brief Class constructor
+         */
+        ComovingPsatdAlgorithm (const SpectralKSpace& spectral_kspace,
+                                const amrex::DistributionMapping& dm,
+                                const int norder_x,
+                                const int norder_y,
+                                const int norder_z,
+                                const bool nodal,
+                                const amrex::Array<amrex::Real,3>& v_comoving,
+                                const amrex::Real dt,
+                                const bool update_with_rho);
+
+        /**
+         * \brief Override the update equations in Fourier space
+         */
+        virtual void pushSpectralFields (SpectralFieldData& f) const override final;
+
+        virtual int getRequiredNumberOfFields () const override final
+        {
+            return SpectralFieldIndex::n_fields;
+        };
+
+        /* \brief Initialize the coefficients needed in the update equations
+         */
+        void InitializeSpectralCoefficients (const SpectralKSpace& spectral_kspace,
+                                             const amrex::DistributionMapping& dm,
+                                             const amrex::Real dt);
+
+        /**
+         * \brief Virtual function for current correction in Fourier space.
+         * This function overrides the virtual function \c CurrentCorrection in the
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
+         *
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
+         */
+        virtual void CurrentCorrection (SpectralFieldData& field_data,
+                                        std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+                                        const std::unique_ptr<amrex::MultiFab>& rho) override final;
+
+        /**
+         * \brief Virtual function for Vay current deposition in Fourier space.
+         * This function overrides the virtual function \c VayDeposition in the
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
+         *
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         */
+        virtual void VayDeposition (SpectralFieldData& field_data,
+                                    std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
+
+    private:
+
+        // Real and complex spectral coefficients
+        SpectralRealCoefficients C_coef, S_ck_coef;
+        SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef;
+
+        // k vectors
+        KVectorComponent kx_vec;
+#if (AMREX_SPACEDIM==3)
+        KVectorComponent ky_vec;
+#endif
+        KVectorComponent kz_vec;
+
+        // Additional member variables
+        amrex::Array<amrex::Real,3> m_v_comoving;
+        amrex::Real m_dt;
+        bool m_update_with_rho;
+};
+
+#endif // WARPX_USE_PSATD
+#endif // WARPX_COMOVING_PSATD_ALGORITHM_H_
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp
new file mode 100644
index 000000000..ca82d5a8c
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp
@@ -0,0 +1,502 @@
+#include "ComovingPsatdAlgorithm.H"
+#include "Utils/WarpXConst.H"
+
+#if WARPX_USE_PSATD
+
+using namespace amrex;
+
+ComovingPsatdAlgorithm::ComovingPsatdAlgorithm (const SpectralKSpace& spectral_kspace,
+                                                const DistributionMapping& dm,
+                                                const int norder_x, const int norder_y,
+                                                const int norder_z, const bool nodal,
+                                                const amrex::Array<amrex::Real, 3>& v_comoving,
+                                                const amrex::Real dt,
+                                                const bool update_with_rho)
+     // Members initialization
+     : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal),
+       // Initialize the infinite-order k vectors (the argument n_order = -1 selects
+       // the infinite order option, the argument nodal = false is then irrelevant)
+       kx_vec(spectral_kspace.getModifiedKComponent(dm, 0, -1, false)),
+#if (AMREX_SPACEDIM==3)
+       ky_vec(spectral_kspace.getModifiedKComponent(dm, 1, -1, false)),
+       kz_vec(spectral_kspace.getModifiedKComponent(dm, 2, -1, false)),
+#else
+       kz_vec(spectral_kspace.getModifiedKComponent(dm, 1, -1, false)),
+#endif
+       m_v_comoving(v_comoving),
+       m_dt(dt),
+       m_update_with_rho(update_with_rho)
+{
+    const BoxArray& ba = spectral_kspace.spectralspace_ba;
+
+    // Allocate arrays of real spectral coefficients
+    C_coef    = SpectralRealCoefficients(ba, dm, 1, 0);
+    S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+
+    // Allocate arrays of complex spectral coefficients
+    X1_coef     = SpectralComplexCoefficients(ba, dm, 1, 0);
+    X2_coef     = SpectralComplexCoefficients(ba, dm, 1, 0);
+    X3_coef     = SpectralComplexCoefficients(ba, dm, 1, 0);
+    X4_coef     = SpectralComplexCoefficients(ba, dm, 1, 0);
+    Theta2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+
+    // Initialize real and complex spectral coefficients
+    InitializeSpectralCoefficients(spectral_kspace, dm, dt);
+}
+
+void
+ComovingPsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
+{
+    // Loop over boxes
+    for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi){
+
+        const amrex::Box& bx = f.fields[mfi].box();
+
+        // Extract arrays for the fields to be updated
+        amrex::Array4<Complex> fields = f.fields[mfi].array();
+
+        // Extract arrays for the coefficients
+        amrex::Array4<const amrex::Real> C_arr    = C_coef   [mfi].array();
+        amrex::Array4<const amrex::Real> S_ck_arr = S_ck_coef[mfi].array();
+        amrex::Array4<const Complex>     X1_arr   = X1_coef  [mfi].array();
+        amrex::Array4<const Complex>     X2_arr   = X2_coef  [mfi].array();
+        amrex::Array4<const Complex>     X3_arr   = X3_coef  [mfi].array();
+        amrex::Array4<const Complex>     X4_arr   = X4_coef  [mfi].array();
+
+        // Extract pointers for the k vectors
+        const amrex::Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr();
+#if (AMREX_SPACEDIM==3)
+        const amrex::Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+#endif
+        const amrex::Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr();
+
+        // Loop over indices within one box
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            // Record old values of the fields to be updated
+            using Idx = SpectralFieldIndex;
+            const Complex Ex_old = fields(i,j,k,Idx::Ex);
+            const Complex Ey_old = fields(i,j,k,Idx::Ey);
+            const Complex Ez_old = fields(i,j,k,Idx::Ez);
+            const Complex Bx_old = fields(i,j,k,Idx::Bx);
+            const Complex By_old = fields(i,j,k,Idx::By);
+            const Complex Bz_old = fields(i,j,k,Idx::Bz);
+
+            // Shortcuts for the values of J and rho
+            const Complex Jx = fields(i,j,k,Idx::Jx);
+            const Complex Jy = fields(i,j,k,Idx::Jy);
+            const Complex Jz = fields(i,j,k,Idx::Jz);
+            const Complex rho_old = fields(i,j,k,Idx::rho_old);
+            const Complex rho_new = fields(i,j,k,Idx::rho_new);
+
+            // k vector values
+            const amrex::Real kx_mod = modified_kx_arr[i];
+#if (AMREX_SPACEDIM==3)
+            const amrex::Real ky_mod = modified_ky_arr[j];
+            const amrex::Real kz_mod = modified_kz_arr[k];
+#else
+            constexpr amrex::Real ky_mod = 0._rt;
+            const     amrex::Real kz_mod = modified_kz_arr[j];
+#endif
+
+            // Physical constant c**2 and imaginary unit
+            constexpr amrex::Real c2 = PhysConst::c*PhysConst::c;
+            constexpr Complex I = Complex{0._rt,1._rt};
+
+            // The definition of these coefficients is explained in more detail
+            // in the function InitializeSpectralCoefficients below
+            const amrex::Real C    = C_arr(i,j,k);
+            const amrex::Real S_ck = S_ck_arr(i,j,k);
+            const Complex     X1   = X1_arr(i,j,k);
+            const Complex     X2   = X2_arr(i,j,k);
+            const Complex     X3   = X3_arr(i,j,k);
+            const Complex     X4   = X4_arr(i,j,k);
+
+            // Update E
+            fields(i,j,k,Idx::Ex) = C*Ex_old + S_ck*c2*I*(ky_mod*Bz_old - kz_mod*By_old)
+                + X4*Jx - I*(X2*rho_new - X3*rho_old)*kx_mod;
+
+            fields(i,j,k,Idx::Ey) = C*Ey_old + S_ck*c2*I*(kz_mod*Bx_old - kx_mod*Bz_old)
+                + X4*Jy - I*(X2*rho_new - X3*rho_old)*ky_mod;
+
+            fields(i,j,k,Idx::Ez) = C*Ez_old + S_ck*c2*I*(kx_mod*By_old - ky_mod*Bx_old)
+                + X4*Jz - I*(X2*rho_new - X3*rho_old)*kz_mod;
+
+            // Update B
+            fields(i,j,k,Idx::Bx) = C*Bx_old - S_ck*I*(ky_mod*Ez_old - kz_mod*Ey_old)
+                + X1*I*(ky_mod*Jz - kz_mod*Jy);
+
+            fields(i,j,k,Idx::By) = C*By_old - S_ck*I*(kz_mod*Ex_old - kx_mod*Ez_old)
+                + X1*I*(kz_mod*Jx - kx_mod*Jz);
+
+            fields(i,j,k,Idx::Bz) = C*Bz_old - S_ck*I*(kx_mod*Ey_old - ky_mod*Ex_old)
+                + X1*I*(kx_mod*Jy - ky_mod*Jx);
+        });
+    }
+}
+
+void ComovingPsatdAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& spectral_kspace,
+                                                             const amrex::DistributionMapping& dm,
+                                                             const amrex::Real dt)
+{
+    const amrex::BoxArray& ba = spectral_kspace.spectralspace_ba;
+
+    // Loop over boxes and allocate the corresponding coefficients for each box
+    for (amrex::MFIter mfi(ba, dm); mfi.isValid(); ++mfi) {
+
+        const amrex::Box& bx = ba[mfi];
+
+        // Extract pointers for the k vectors
+        const amrex::Real* kx_mod = modified_kx_vec[mfi].dataPtr();
+        const amrex::Real* kx     = kx_vec[mfi].dataPtr();
+#if (AMREX_SPACEDIM==3)
+        const amrex::Real* ky_mod = modified_ky_vec[mfi].dataPtr();
+        const amrex::Real* ky     = ky_vec[mfi].dataPtr();
+#endif
+        const amrex::Real* kz_mod = modified_kz_vec[mfi].dataPtr();
+        const amrex::Real* kz     = kz_vec[mfi].dataPtr();
+
+        // Extract arrays for the coefficients
+        amrex::Array4<amrex::Real> C    = C_coef     [mfi].array();
+        amrex::Array4<amrex::Real> S_ck = S_ck_coef  [mfi].array();
+        amrex::Array4<Complex>     X1   = X1_coef    [mfi].array();
+        amrex::Array4<Complex>     X2   = X2_coef    [mfi].array();
+        amrex::Array4<Complex>     X3   = X3_coef    [mfi].array();
+        amrex::Array4<Complex>     X4   = X4_coef    [mfi].array();
+        amrex::Array4<Complex>     T2   = Theta2_coef[mfi].array();
+
+        // Store comoving velocity
+        const amrex::Real vx = m_v_comoving[0];
+#if (AMREX_SPACEDIM==3)
+        const amrex::Real vy = m_v_comoving[1];
+#endif
+        const amrex::Real vz = m_v_comoving[2];
+
+        // Loop over indices within one box
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            // Calculate norm of finite-order k vector
+            const amrex::Real knorm_mod = std::sqrt(
+                std::pow(kx_mod[i], 2) +
+#if (AMREX_SPACEDIM==3)
+                std::pow(ky_mod[j], 2) +
+                std::pow(kz_mod[k], 2));
+#else
+                std::pow(kz_mod[j], 2));
+#endif
+            // Calculate norm of infinite-order k vector
+            const amrex::Real knorm = std::sqrt(
+                std::pow(kx[i], 2) +
+#if (AMREX_SPACEDIM==3)
+                std::pow(ky[j], 2) +
+                std::pow(kz[k], 2));
+#else
+                std::pow(kz[j], 2));
+#endif
+            // Physical constants c, c**2, and epsilon_0, and imaginary unit
+            constexpr amrex::Real c   = PhysConst::c;
+            constexpr amrex::Real c2  = c*c;
+            constexpr amrex::Real ep0 = PhysConst::ep0;
+            constexpr Complex     I   = Complex{0._rt, 1._rt};
+
+            // Auxiliary coefficients used when update_with_rho=false
+            const amrex::Real dt2 = dt * dt;
+
+            // Calculate dot product of k vector with comoving velocity
+            const amrex::Real kv = kx[i]*vx +
+#if (AMREX_SPACEDIM==3)
+                ky[j]*vy + kz[k]*vz;
+#else
+                kz[j]*vz;
+#endif
+
+            if (knorm_mod != 0. && knorm != 0.) {
+
+                // Auxiliary coefficients
+                const amrex::Real om_mod  = c * knorm_mod;
+                const amrex::Real om2_mod = om_mod * om_mod;
+                const amrex::Real om  = c * knorm;
+                const amrex::Real om2 = om * om;
+                const Complex tmp1 = amrex::exp(  I * om_mod * dt);
+                const Complex tmp2 = amrex::exp(- I * om_mod * dt);
+                const Complex tmp1_sqrt = amrex::exp(  I * om_mod * dt * 0.5_rt);
+                const Complex tmp2_sqrt = amrex::exp(- I * om_mod * dt * 0.5_rt);
+
+                C   (i,j,k) = std::cos(om_mod * dt);
+                S_ck(i,j,k) = std::sin(om_mod * dt) / om_mod;
+
+                const amrex::Real nu = - kv / om;
+                const Complex theta      = amrex::exp(  I * nu * om * dt * 0.5_rt);
+                const Complex theta_star = amrex::exp(- I * nu * om * dt * 0.5_rt);
+
+                T2(i,j,k) = theta * theta;
+
+                if ( (nu != om_mod/om) && (nu != -om_mod/om) && (nu != 0.) ) {
+
+                    Complex x1 = om2 / (om2_mod - nu * nu * om2)
+                        * (theta_star - theta * C(i,j,k) + I * nu * om * theta * S_ck(i,j,k));
+
+                    // X1 multiplies i*(k \times J) in the update equation for B
+                    X1(i,j,k) = x1 / (ep0 * om2);
+
+                    // X2 multiplies rho_new in the update equation for E
+                    // X3 multiplies rho_old in the update equation for E
+                    X2(i,j,k) = c2 * (x1 * om2_mod - theta * (1._rt - C(i,j,k)) * om2)
+                        / (theta_star - theta) / (ep0 * om2 * om2_mod);
+                    X3(i,j,k) = c2 * (x1 * om2_mod - theta_star * (1._rt - C(i,j,k)) * om2)
+                        / (theta_star - theta) / (ep0 * om2 * om2_mod);
+
+                    // X4 multiplies J in the update equation for E
+                    X4(i,j,k) = I * nu * om * X1(i,j,k) - theta * S_ck(i,j,k) / ep0;
+                }
+
+                // Limits for nu = 0
+                if (nu == 0.) {
+
+                    // X1 multiplies i*(k \times J) in the update equation for B
+                    X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2_mod);
+
+                    // X2 multiplies rho_new in the update equation for E
+                    // X3 multiplies rho_old in the update equation for E
+                    X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2_mod);
+                    X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2_mod);
+
+                    // Coefficient multiplying J in update equation for E
+                    X4(i,j,k) = - S_ck(i,j,k) / ep0;
+                }
+
+                // Limits for nu = omega_mod/omega
+                if (nu == om_mod/om) {
+
+                    // X1 multiplies i*(k \times J) in the update equation for B
+                    X1(i,j,k) = tmp1_sqrt * (1._rt - tmp2 * tmp2 - 2._rt * I * om_mod * dt) / (4._rt * ep0 * om2_mod);
+
+                    // X2 multiplies rho_new in the update equation for E
+                    // X3 multiplies rho_old in the update equation for E
+                    X2(i,j,k) = c2 * (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * om_mod * dt * tmp1)
+                        / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt));
+                    X3(i,j,k) = c2 * (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * om_mod * dt * tmp1)
+                        / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt));
+
+                    // Coefficient multiplying J in update equation for E
+                    X4(i,j,k) = tmp1_sqrt * (I - I * tmp2 * tmp2 - 2._rt * om_mod * dt) / (4._rt * ep0 * om_mod);
+                }
+
+                // Limits for nu = -omega_mod/omega
+                if (nu == -om_mod/om) {
+
+                    // X1 multiplies i*(k \times J) in the update equation for B
+                    X1(i,j,k) = tmp2_sqrt * (1._rt - tmp1 * tmp1 + 2._rt * I * om_mod * dt) / (4._rt * ep0 * om2_mod);
+
+                    // X2 multiplies rho_new in the update equation for E
+                    // X3 multiplies rho_old in the update equation for E
+                    X2(i,j,k) = c2 * (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * om_mod * dt)
+                        / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt));
+                    X3(i,j,k) = c2 * (3._rt - 2._rt * tmp2 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * om_mod * dt)
+                        / (4._rt * ep0 * om2_mod * (tmp1 - 1._rt));
+
+                    // Coefficient multiplying J in update equation for E
+                    X4(i,j,k) = tmp2_sqrt * (- I + I * tmp1 * tmp1 - 2._rt * om_mod * dt) / (4._rt * ep0 * om_mod);
+                }
+            }
+
+            // Limits for omega = 0 only
+            else if (knorm_mod != 0. && knorm == 0.) {
+
+                const amrex::Real om_mod  = c * knorm_mod;
+                const amrex::Real om2_mod = om_mod * om_mod;
+
+                C   (i,j,k) = std::cos(om_mod * dt);
+                S_ck(i,j,k) = std::sin(om_mod * dt) / om_mod;
+                T2(i,j,k) = 1._rt;
+
+                // X1 multiplies i*(k \times J) in the update equation for B
+                X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2_mod);
+
+                // X2 multiplies rho_new in the update equation for E
+                // X3 multiplies rho_old in the update equation for E
+                X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2_mod);
+                X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2_mod);
+
+                // Coefficient multiplying J in update equation for E
+                X4(i,j,k) = - S_ck(i,j,k) / ep0;
+
+            }
+
+            // Limits for omega_mod = 0 only
+            else if (knorm_mod == 0. && knorm != 0.) {
+
+                const amrex::Real om  = c * knorm;
+                const amrex::Real om2 = om * om;
+                const amrex::Real nu = - kv / om;
+                const Complex theta      = amrex::exp(I * nu * om * dt * 0.5_rt);
+                const Complex theta_star = amrex::exp(- I * nu * om * dt * 0.5_rt);
+
+                C(i,j,k) = 1._rt;
+                S_ck(i,j,k) = dt;
+                T2(i,j,k) = theta * theta;
+
+                if (nu != 0.) {
+                    // X1 multiplies i*(k \times J) in the update equation for B
+                    X1(i,j,k) = (-theta_star + theta - I * nu * om * dt * theta)
+                        / (ep0 * nu * nu * om2);
+
+                    // X2 multiplies rho_new in the update equation for E
+                    // X3 multiplies rho_old in the update equation for E
+                    X2(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om * dt * T2(i,j,k)
+                        + 0.5_rt * nu * nu * om2 * dt * dt * T2(i,j,k))
+                        / (ep0 * nu * nu * om2 * (T2(i,j,k) - 1._rt));
+                    X3(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om * dt * T2(i,j,k)
+                        + 0.5_rt * nu * nu * om2 * dt * dt)
+                        / (ep0 * nu * nu * om2 * (T2(i,j,k) - 1._rt));
+
+                    // Coefficient multiplying J in update equation for E
+                    X4(i,j,k) = I * (theta - theta_star) / (ep0 * nu * om);
+                }
+
+                else {
+                    // X1 multiplies i*(k \times J) in the update equation for B
+                    X1(i,j,k) = dt2 / (2._rt * ep0);
+
+                    // X2 multiplies rho_new in the update equation for E
+                    // X3 multiplies rho_old in the update equation for E
+                    X2(i,j,k) = c2 * dt2 / (6._rt * ep0);
+                    X3(i,j,k) = - c2 * dt2 / (3._rt * ep0);
+
+                    // Coefficient multiplying J in update equation for E
+                    X4(i,j,k) = -dt / ep0;
+                }
+            }
+
+            // Limits for omega_mod = 0 and omega = 0
+            else if (knorm_mod == 0. && knorm == 0.) {
+
+                C(i,j,k) = 1._rt;
+                S_ck(i,j,k) = dt;
+                T2(i,j,k) = 1._rt;
+
+                // X1 multiplies i*(k \times J) in the update equation for B
+                X1(i,j,k) = dt2 / (2._rt * ep0);
+
+                // X2 multiplies rho_new in the update equation for E
+                // X3 multiplies rho_old in the update equation for E
+                X2(i,j,k) = c2 * dt2 / (6._rt * ep0);
+                X3(i,j,k) = - c2 * dt2 / (3._rt * ep0);
+
+                // Coefficient multiplying J in update equation for E
+                X4(i,j,k) = -dt / ep0;
+            }
+        });
+    }
+}
+
+void
+ComovingPsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
+                                           std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+                                           const std::unique_ptr<amrex::MultiFab>& rho)
+{
+    // Profiling
+    BL_PROFILE("ComovingAlgorithm::CurrentCorrection");
+
+    using Idx = SpectralFieldIndex;
+
+    // Forward Fourier transform of J and rho
+    field_data.ForwardTransform(*current[0], Idx::Jx, 0);
+    field_data.ForwardTransform(*current[1], Idx::Jy, 0);
+    field_data.ForwardTransform(*current[2], Idx::Jz, 0);
+    field_data.ForwardTransform(*rho, Idx::rho_old, 0);
+    field_data.ForwardTransform(*rho, Idx::rho_new, 1);
+
+    // Loop over boxes
+    for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
+
+        const amrex::Box& bx = field_data.fields[mfi].box();
+
+        // Extract arrays for the fields to be updated
+        amrex::Array4<Complex> fields = field_data.fields[mfi].array();
+
+        // Extract pointers for the k vectors
+        const amrex::Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr();
+        const amrex::Real* const kx_arr = kx_vec[mfi].dataPtr();
+#if (AMREX_SPACEDIM==3)
+        const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+        const amrex::Real* const ky_arr = ky_vec[mfi].dataPtr();
+#endif
+        const amrex::Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr();
+        const amrex::Real* const kz_arr = kz_vec[mfi].dataPtr();
+
+        // Local copy of member variables before GPU loop
+        const amrex::Real dt = m_dt;
+
+        // Comoving velocity
+        const amrex::Real vx = m_v_comoving[0];
+        const amrex::Real vy = m_v_comoving[1];
+        const amrex::Real vz = m_v_comoving[2];
+
+        // Loop over indices within one box
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            // Shortcuts for the values of J and rho
+            const Complex Jx = fields(i,j,k,Idx::Jx);
+            const Complex Jy = fields(i,j,k,Idx::Jy);
+            const Complex Jz = fields(i,j,k,Idx::Jz);
+            const Complex rho_old = fields(i,j,k,Idx::rho_old);
+            const Complex rho_new = fields(i,j,k,Idx::rho_new);
+
+            // k vector values, and coefficients
+            const amrex::Real kx_mod = modified_kx_arr[i];
+            const amrex::Real kx = kx_arr[i];
+#if (AMREX_SPACEDIM==3)
+            const amrex::Real ky_mod = modified_ky_arr[j];
+            const amrex::Real kz_mod = modified_kz_arr[k];
+            const amrex::Real ky = ky_arr[j];
+            const amrex::Real kz = kz_arr[k];
+#else
+            constexpr amrex::Real ky_mod = 0._rt;
+            const     amrex::Real kz_mod = modified_kz_arr[j];
+            constexpr amrex::Real ky = 0._rt;
+            const     amrex::Real kz = kz_arr[j];
+#endif
+            constexpr Complex I = Complex{0._rt,1._rt};
+
+            const amrex::Real knorm_mod = std::sqrt(kx_mod * kx_mod + ky_mod * ky_mod + kz_mod * kz_mod);
+
+            // Correct J
+            if (knorm_mod != 0._rt)
+            {
+                const Complex kmod_dot_J = kx_mod * Jx + ky_mod * Jy + kz_mod * Jz;
+                const amrex::Real k_dot_v = kx * vx + ky * vy + kz * vz;
+
+                if ( k_dot_v != 0._rt ) {
+
+                    const Complex theta = amrex::exp(- I * k_dot_v * dt * 0.5_rt);
+                    const Complex den = 1._rt - theta * theta;
+
+                    fields(i,j,k,Idx::Jx) = Jx - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kx_mod / (knorm_mod * knorm_mod);
+                    fields(i,j,k,Idx::Jy) = Jy - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * ky_mod / (knorm_mod * knorm_mod);
+                    fields(i,j,k,Idx::Jz) = Jz - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kz_mod / (knorm_mod * knorm_mod);
+
+                } else {
+
+                    fields(i,j,k,Idx::Jx) = Jx - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kx_mod / (knorm_mod * knorm_mod);
+                    fields(i,j,k,Idx::Jy) = Jy - (kmod_dot_J - I * (rho_new - rho_old) / dt) * ky_mod / (knorm_mod * knorm_mod);
+                    fields(i,j,k,Idx::Jz) = Jz - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kz_mod / (knorm_mod * knorm_mod);
+                }
+            }
+        });
+    }
+
+    // Backward Fourier transform of J
+    field_data.BackwardTransform(*current[0], Idx::Jx, 0);
+    field_data.BackwardTransform(*current[1], Idx::Jy, 0);
+    field_data.BackwardTransform(*current[2], Idx::Jz, 0);
+}
+
+void
+ComovingPsatdAlgorithm::VayDeposition (SpectralFieldData& /*field_data*/,
+                                       std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/)
+{
+    amrex::Abort("Vay deposition not implemented for comoving PSATD");
+}
+
+#endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
index 83ae83d69..cc43f0489 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
@@ -59,6 +59,12 @@ class GalileanAlgorithm : public SpectralBaseAlgorithm
     private:
         SpectralRealCoefficients C_coef, S_ck_coef;
         SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef;
+        // Centered modified finite-order k vectors
+        KVectorComponent modified_kx_vec_centered;
+#if (AMREX_SPACEDIM==3)
+        KVectorComponent modified_ky_vec_centered;
+#endif
+        KVectorComponent modified_kz_vec_centered;
         amrex::Array<amrex::Real,3> m_v_galilean;
         amrex::Real m_dt;
         bool m_update_with_rho;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
index 0757d17cd..ea5aa7ec0 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
@@ -18,6 +18,16 @@ GalileanAlgorithm::GalileanAlgorithm(const SpectralKSpace& spectral_kspace,
                          const bool update_with_rho)
      // Initialize members of base class
      : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal),
+       // Initialize the centered finite-order modified k vectors: these are computed
+       // always with the assumption of centered grids (argument nodal = true),
+       // for both nodal and staggered simulations
+       modified_kx_vec_centered(spectral_kspace.getModifiedKComponent(dm,0,norder_x,true)),
+#if (AMREX_SPACEDIM==3)
+       modified_ky_vec_centered(spectral_kspace.getModifiedKComponent(dm,1,norder_y,true)),
+       modified_kz_vec_centered(spectral_kspace.getModifiedKComponent(dm,2,norder_z,true)),
+#else
+       modified_kz_vec_centered(spectral_kspace.getModifiedKComponent(dm,1,norder_z,true)),
+#endif
        m_v_galilean(v_galilean),
        m_dt(dt),
        m_update_with_rho(update_with_rho)
@@ -91,8 +101,8 @@ GalileanAlgorithm::pushSpectralFields (SpectralFieldData& f) const
             const Real ky = modified_ky_arr[j];
             const Real kz = modified_kz_arr[k];
 #else
-            constexpr Real ky = 0;
-            const Real kz = modified_kz_arr[j];
+            constexpr Real ky = 0._rt;
+            const     Real kz = modified_kz_arr[j];
 #endif
             // Physical constant c**2 and imaginary unit
             constexpr Real c2   = PhysConst::c*PhysConst::c;
@@ -168,11 +178,14 @@ void GalileanAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& sp
         const Box& bx = ba[mfi];
 
         // Extract pointers for the k vectors
-        const Real* modified_kx = modified_kx_vec[mfi].dataPtr();
+        const Real* kx   = modified_kx_vec[mfi].dataPtr();
+        const Real* kx_c = modified_kx_vec_centered[mfi].dataPtr();
 #if (AMREX_SPACEDIM==3)
-        const Real* modified_ky = modified_ky_vec[mfi].dataPtr();
+        const Real* ky   = modified_ky_vec[mfi].dataPtr();
+        const Real* ky_c = modified_ky_vec_centered[mfi].dataPtr();
 #endif
-        const Real* modified_kz = modified_kz_vec[mfi].dataPtr();
+        const Real* kz   = modified_kz_vec[mfi].dataPtr();
+        const Real* kz_c = modified_kz_vec_centered[mfi].dataPtr();
 
         // Extract arrays for the coefficients
         Array4<Real> C = C_coef[mfi].array();
@@ -194,13 +207,22 @@ void GalileanAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& sp
         ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             // Calculate norm of vector
-            const Real k_norm = std::sqrt(
-                std::pow(modified_kx[i], 2) +
+            const Real knorm = std::sqrt(
+                std::pow(kx[i], 2) +
 #if (AMREX_SPACEDIM==3)
-                std::pow(modified_ky[j], 2) +
-                std::pow(modified_kz[k], 2));
+                std::pow(ky[j], 2) +
+                std::pow(kz[k], 2));
 #else
-                std::pow(modified_kz[j], 2));
+                std::pow(kz[j], 2));
+#endif
+            // Calculate norm of vector
+            const Real knorm_c = std::sqrt(
+                std::pow(kx_c[i], 2) +
+#if (AMREX_SPACEDIM==3)
+                std::pow(ky_c[j], 2) +
+                std::pow(kz_c[k], 2));
+#else
+                std::pow(kz_c[j], 2));
 #endif
             // Physical constants c, c**2, and epsilon_0, and imaginary unit
             constexpr Real c    = PhysConst::c;
@@ -213,142 +235,222 @@ void GalileanAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& sp
             const Real dt3 = dt * dt2;
             Complex X2_old, X3_old;
 
-            // Calculate dot product of k vector with Galilean velocity
-            const Real kv = modified_kx[i]*vx +
+            // Calculate dot product of k vector with Galilean velocity:
+            // this has to be computed always with the centered finite-order modified k
+            // vectors, in order to work correctly for both nodal and staggered simulations
+            const Real kv = kx_c[i]*vx +
 #if (AMREX_SPACEDIM==3)
-                modified_ky[j]*vy + modified_kz[k]*vz;
+                ky_c[j]*vy + kz_c[k]*vz;
 #else
-                modified_kz[j]*vz;
+                kz_c[j]*vz;
 #endif
             // The coefficients in the following refer to the ones given in equations
             // (12a)-(12d) of (Lehe et al, PRE 94, 2016), used to update B and E
             // (equations (11a) and (11b) of the same reference, respectively)
 
-            if (k_norm != 0.) {
+            if (knorm != 0. && knorm_c != 0.) {
 
                 // Auxiliary coefficients
-                const Real    k2   = k_norm * k_norm;
-                const Real    ck   = c * k_norm;
-                const Real    ckdt = ck * dt;
-                const Complex tmp1 = amrex::exp(  I * ckdt); // limit of T2 for nu =  1
-                const Complex tmp2 = amrex::exp(- I * ckdt); // limit of T2 for nu = -1
+                const Real om  = c * knorm;
+                const Real om2 = om * om;
+                const Real om3 = om * om2;
+                const Real om_c  = c * knorm_c;
+                const Real om2_c = om_c * om_c;
+                const Real om3_c = om_c * om2_c;
+                const Complex tmp1 = amrex::exp(  I * om * dt);
+                const Complex tmp2 = amrex::exp(- I * om * dt);
 
                 // See equation (12a)
-                C   (i,j,k) = std::cos(ckdt);
-                S_ck(i,j,k) = std::sin(ckdt) / ck;
+                C   (i,j,k) = std::cos(om * dt);
+                S_ck(i,j,k) = std::sin(om * dt) / om;
 
                 // See equation (12b)
-                const Real    nu         = kv / ck;
-                const Complex theta      = amrex::exp(  I * 0.5_rt * kv * dt);
-                const Complex theta_star = amrex::exp(- I * 0.5_rt * kv * dt);
+                const Real    nu = kv / om_c;
+                const Complex theta      = amrex::exp(  I * nu * om_c * dt * 0.5_rt);
+                const Complex theta_star = amrex::exp(- I * nu * om_c * dt * 0.5_rt);
 
                 // This is exp(i*(k \dot v_gal)*dt)
                 T2(i,j,k) = theta * theta;
 
-                if ( (nu != 1.) && (nu != 0.) ) {
+                if ( (nu != om/om_c) && (nu != -om/om_c) && (nu != 0.) ) {
 
                     // x1 is the coefficient chi_1 in equation (12c)
-                    Complex x1 = 1._rt / (1._rt - nu*nu)
-                        * (theta_star - C(i,j,k) * theta + I * kv * S_ck(i,j,k) * theta);
+                    Complex x1 = om2_c / (om2 - nu * nu * om2_c)
+                        * (theta_star - theta * C(i,j,k) + I * nu * om_c * theta * S_ck(i,j,k));
 
                     // X1 multiplies i*(k \times J) in the update equation for B
-                    X1(i,j,k) = theta * x1 / (ep0 * c2 * k2);
+                    X1(i,j,k) = theta * x1 / (ep0 * om2_c);
 
                     if (update_with_rho) {
                         // X2 multiplies rho_new in the update equation for E
                         // X3 multiplies rho_old in the update equation for E
-                        X2(i,j,k) = (x1 - theta * (1._rt - C(i,j,k))) / (theta_star - theta) / (ep0 * k2);
-                        X3(i,j,k) = (x1 - theta_star * (1._rt - C(i,j,k))) / (theta_star - theta) / (ep0 * k2);
+                        X2(i,j,k) = c2 * (x1 * om2 - theta * (1._rt - C(i,j,k)) * om2_c)
+                            / (theta_star - theta) / (ep0 * om2_c * om2);
+                        X3(i,j,k) = c2 * (x1 * om2 - theta_star * (1._rt - C(i,j,k)) * om2_c)
+                            / (theta_star - theta) / (ep0 * om2_c * om2);
                     } else {
                         // X2_old is the coefficient chi_2 in equation (12d)
                         // X3_old is the coefficient chi_3 in equation (12d)
                         // X2 multiplies (k \dot E) in the update equation for E
                         // X3 multiplies (k \dot J) in the update equation for E
-                        X2_old = (x1 - theta * (1._rt - C(i,j,k))) / (theta_star - theta);
-                        X3_old = (x1 - theta_star * (1._rt - C(i,j,k))) / (theta_star - theta);
-                        X2(i,j,k) = T2(i,j,k) * (X2_old - X3_old) / k2;
-                        X3(i,j,k) = I * X2_old * (T2(i,j,k) - 1._rt) / (ep0 * k2 * kv);
+                        X2_old = (x1 * om2 - theta * (1._rt - C(i,j,k)) * om2_c)
+                            / (theta_star - theta);
+                        X3_old = (x1 * om2 - theta_star * (1._rt - C(i,j,k)) * om2_c)
+                            / (theta_star - theta);
+                        X2(i,j,k) = c2 * T2(i,j,k) * (X2_old - X3_old)
+                            / (om2_c * om2);
+                        X3(i,j,k) = I * c2 * X2_old * (T2(i,j,k) - 1._rt)
+                            / (ep0 * nu * om3_c * om2);
                     }
 
                     // X4 multiplies J in the update equation for E
-                    X4(i,j,k) = I * kv * X1(i,j,k) - T2(i,j,k) * S_ck(i,j,k) / ep0;
+                    X4(i,j,k) = I * nu * om_c * X1(i,j,k) - T2(i,j,k) * S_ck(i,j,k) / ep0;
                 }
 
                 // Limits for nu = 0
                 if (nu == 0.) {
 
                     // X1 multiplies i*(k \times J) in the update equation for B
-                    X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * c2 * k2);
+                    X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2);
 
                     if (update_with_rho) {
                         // X2 multiplies rho_new in the update equation for E
                         // X3 multiplies rho_old in the update equation for E
-                        X2(i,j,k) = (1._rt - S_ck(i,j,k) / dt) / (ep0 * k2);
-                        X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * k2);
+                        X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2);
+                        X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2);
                     } else {
                         // X2 multiplies (k \dot E) in the update equation for E
                         // X3 multiplies (k \dot J) in the update equation for E
-                        X2(i,j,k) = (1._rt - C(i,j,k)) / k2;
-                        X3(i,j,k) = (S_ck(i,j,k) / dt - 1._rt) * dt / (ep0 * k2);
+                        X2(i,j,k) = c2 * (1._rt - C(i,j,k)) / om2;
+                        X3(i,j,k) = c2 * (S_ck(i,j,k) / dt - 1._rt) * dt / (ep0 * om2);
                     }
 
                     // Coefficient multiplying J in update equation for E
                     X4(i,j,k) = - S_ck(i,j,k) / ep0;
                 }
 
-                // Limits for nu = 1
-                if (nu == 1.) {
+                // Limits for nu = omega/omega_c
+                if (nu == om/om_c) {
 
                     // X1 multiplies i*(k \times J) in the update equation for B
-                    X1(i,j,k) = (1._rt - tmp1 * tmp1 + 2._rt * I * ckdt) / (4._rt * ep0 * c2 * k2);
+                    X1(i,j,k) = (1._rt - tmp1 * tmp1 + 2._rt * I * om * dt) / (4._rt * ep0 * om2);
 
                     if (update_with_rho) {
                         // X2 multiplies rho_new in the update equation for E
                         // X3 multiplies rho_old in the update equation for E
-                        X2(i,j,k) = (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * ckdt)
-                            / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
-                        X3(i,j,k) = (3._rt - 2._rt / tmp1 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * ckdt)
-                            / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
+                        X2(i,j,k) = c2 * (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * om * dt)
+                            / (4._rt * ep0 * om2 * (tmp1 - 1._rt));
+                        X3(i,j,k) = c2 * (3._rt - 2._rt * tmp2 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * om * dt)
+                            / (4._rt * ep0 * om2 * (tmp1 - 1._rt));
                     } else {
                         // X2 multiplies (k \dot E) in the update equation for E
                         // X3 multiplies (k \dot J) in the update equation for E
-                        X2(i,j,k) = (1._rt - C(i,j,k)) * tmp1 / k2;
-                        X3(i,j,k) = (2._rt * ckdt - I * tmp1 * tmp1 + 4._rt * I * tmp1 - 3._rt * I)
-                            / (4._rt * ep0 * ck * k2);
+                        X2(i,j,k) = c2 * (1._rt - C(i,j,k)) * tmp1 / om2;
+                        X3(i,j,k) = c2 * (2._rt * om * dt - I * tmp1 * tmp1 + 4._rt * I * tmp1 - 3._rt * I)
+                            / (4._rt * ep0 * om3);
                     }
 
                     // Coefficient multiplying J in update equation for E
-                    X4(i,j,k) = (- I + I * tmp1 * tmp1 - 2._rt * ckdt) / (4._rt * ep0 * ck);
+                    X4(i,j,k) = (- I + I * tmp1 * tmp1 - 2._rt * om * dt) / (4._rt * ep0 * om);
                 }
 
-                // Limits for nu = -1
-                if (nu == -1.) {
+                // Limits for nu = -omega/omega_c
+                if (nu == -om/om_c) {
 
                     // X1 multiplies i*(k \times J) in the update equation for B
-                    X1(i,j,k) = (1._rt - tmp2 * tmp2 - 2._rt * I * ckdt) / (4._rt * ep0 * c2 * k2);
+                    X1(i,j,k) = (1._rt - tmp2 * tmp2 - 2._rt * I * om * dt) / (4._rt * ep0 * om2);
 
                     if (update_with_rho) {
                         // X2 multiplies rho_new in the update equation for E
                         // X3 multiplies rho_old in the update equation for E
-                        X2(i,j,k) = (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * ckdt * tmp1)
-                            / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
-                        X3(i,j,k) = (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * ckdt * tmp1)
-                            / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
+                        X2(i,j,k) = c2 * (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * om * dt * tmp1)
+                            / (4._rt * ep0 * om2 * (tmp1 - 1._rt));
+                        X3(i,j,k) = c2 * (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * om * dt * tmp1)
+                            / (4._rt * ep0 * om2 * (tmp1 - 1._rt));
                     } else {
                         // X2 multiplies (k \dot E) in the update equation for E
                         // X3 multiplies (k \dot J) in the update equation for E
-                        X2(i,j,k) = (1._rt - C(i,j,k)) * tmp2 / k2;
-                        X3(i,j,k) = (2._rt * ckdt + I * tmp2 * tmp2 - 4._rt * I * tmp2 + 3._rt * I)
-                            / (4._rt * ep0 * ck * k2);
+                        X2(i,j,k) = c2 * (1._rt - C(i,j,k)) * tmp2 / om2;
+                        X3(i,j,k) = c2 * (2._rt * om * dt + I * tmp2 * tmp2 - 4._rt * I * tmp2 + 3._rt * I)
+                            / (4._rt * ep0 * om3);
                     }
 
                     // Coefficient multiplying J in update equation for E
-                    X4(i,j,k) = (I - I * tmp2 * tmp2 - 2._rt * ckdt) / (4._rt * ep0 * ck);
+                    X4(i,j,k) = (I - I * tmp2 * tmp2 - 2._rt * om * dt) / (4._rt * ep0 * om);
+                }
+            }
+
+            // Limits for omega_c = 0 only
+            else if (knorm != 0. && knorm_c == 0.) {
+
+                const Real om  = c * knorm;
+                const Real om2 = om * om;
+
+                C   (i,j,k) = std::cos(om * dt);
+                S_ck(i,j,k) = std::sin(om * dt) / om;
+                T2(i,j,k) = 1._rt;
+
+                // X1 multiplies i*(k \times J) in the update equation for B
+                X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2);
+
+                if (update_with_rho) {
+                    // X2 multiplies rho_new in the update equation for E
+                    // X3 multiplies rho_old in the update equation for E
+                    X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2);
+                    X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2);
+                } else {
+                    // X2 multiplies (k \dot E) in the update equation for E
+                    // X3 multiplies (k \dot J) in the update equation for E
+                    X2(i,j,k) = c2 * (1._rt - C(i,j,k)) / om2;
+                    X3(i,j,k) = c2 * (S_ck(i,j,k) / dt - 1._rt) * dt / (ep0 * om2);
                 }
+
+                // Coefficient multiplying J in update equation for E
+                X4(i,j,k) = - S_ck(i,j,k) / ep0;
+
+            }
+
+            // Limits for omega = 0 only
+            else if (knorm == 0. && knorm_c != 0.) {
+
+                const Real om_c  = c * knorm_c;
+                const Real om2_c = om_c * om_c;
+                const Real om3_c = om_c * om2_c;
+                const Real nu = kv / om_c;
+                const Complex theta = amrex::exp(I * nu * om_c * dt * 0.5_rt);
+
+                C(i,j,k) = 1._rt;
+                S_ck(i,j,k) = dt;
+                T2(i,j,k) = theta * theta;
+
+                // X1 multiplies i*(k \times J) in the update equation for B
+                X1(i,j,k) = (-1._rt + T2(i,j,k) - I * nu * om_c * dt * T2(i,j,k))
+                    / (ep0 * nu * nu * om2_c);
+
+                if (update_with_rho) {
+                    // X2 multiplies rho_new in the update equation for E
+                    // X3 multiplies rho_old in the update equation for E
+                    X2(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om_c * dt * T2(i,j,k)
+                        + 0.5_rt * nu * nu * om2_c * dt * dt * T2(i,j,k))
+                        / (ep0 * nu * nu * om2_c * (T2(i,j,k) - 1._rt));
+                    X3(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om_c * dt * T2(i,j,k)
+                        + 0.5_rt * nu * nu * om2_c * dt * dt)
+                        / (ep0 * nu * nu * om2_c * (T2(i,j,k) - 1._rt));
+                } else {
+                    // X2 multiplies (k \dot E) in the update equation for E
+                    // X3 multiplies (k \dot J) in the update equation for E
+                    X2(i,j,k) = c2 * dt * dt * T2(i,j,k) * 0.5_rt;
+                    X3(i,j,k) = c2 * (2._rt * I - 2._rt * nu * om_c * dt * T2(i,j,k)
+                        + I * nu * nu * om2_c * dt * dt * T2(i,j,k))
+                        / (2._rt * ep0 * nu * nu * nu * om3_c);
+                }
+
+                // Coefficient multiplying J in update equation for E
+                X4(i,j,k) = I * (T2(i,j,k) - 1._rt) / (ep0 * nu * om_c);
             }
 
-            // Limits for k = 0
-            else {
+            // Limits for omega = 0 and omega_c = 0
+            else if (knorm == 0. && knorm_c == 0.) {
 
                 // Limits of cos(c*k*dt) and sin(c*k*dt)/(c*k)
                 C(i,j,k) = 1._rt;
@@ -404,11 +506,14 @@ GalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
         amrex::Array4<Complex> fields = field_data.fields[mfi].array();
 
         // Extract pointers for the k vectors
-        const amrex::Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr();
+        const amrex::Real* const modified_kx_arr   = modified_kx_vec[mfi].dataPtr();
+        const amrex::Real* const modified_kx_arr_c = modified_kx_vec_centered[mfi].dataPtr();
 #if (AMREX_SPACEDIM==3)
-        const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+        const amrex::Real* const modified_ky_arr   = modified_ky_vec[mfi].dataPtr();
+        const amrex::Real* const modified_ky_arr_c = modified_ky_vec_centered[mfi].dataPtr();
 #endif
-        const amrex::Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr();
+        const amrex::Real* const modified_kz_arr   = modified_kz_vec[mfi].dataPtr();
+        const amrex::Real* const modified_kz_arr_c = modified_kz_vec_centered[mfi].dataPtr();
 
         // Local copy of member variables before GPU loop
         const amrex::Real dt = m_dt;
@@ -429,13 +534,18 @@ GalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
             const Complex rho_new = fields(i,j,k,Idx::rho_new);
 
             // k vector values, and coefficients
-            const amrex::Real kx = modified_kx_arr[i];
+            const amrex::Real kx   = modified_kx_arr[i];
+            const amrex::Real kx_c = modified_kx_arr_c[i];
 #if (AMREX_SPACEDIM==3)
-            const amrex::Real ky = modified_ky_arr[j];
-            const amrex::Real kz = modified_kz_arr[k];
+            const amrex::Real ky   = modified_ky_arr[j];
+            const amrex::Real kz   = modified_kz_arr[k];
+            const amrex::Real ky_c = modified_ky_arr_c[j];
+            const amrex::Real kz_c = modified_kz_arr_c[k];
 #else
-            constexpr amrex::Real ky = 0._rt;
-            const amrex::Real kz = modified_kz_arr[j];
+            constexpr amrex::Real ky   = 0._rt;
+            const     amrex::Real kz   = modified_kz_arr[j];
+            constexpr amrex::Real ky_c = 0._rt;
+            const     amrex::Real kz_c = modified_kz_arr_c[j];
 #endif
             constexpr Complex I = Complex{0._rt,1._rt};
 
@@ -445,7 +555,7 @@ GalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
             if (k_norm != 0._rt)
             {
                 const Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz;
-                const amrex::Real k_dot_vg = kx * vgx + ky * vgy + kz * vgz;
+                const amrex::Real k_dot_vg = kx_c * vgx + ky_c * vgy + kz_c * vgz;
 
                 if ( k_dot_vg != 0._rt ) {
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
index 57f8c5e3f..71b97b8bc 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
@@ -3,6 +3,7 @@ CEXE_sources += PsatdAlgorithm.cpp
 CEXE_sources += GalileanAlgorithm.cpp
 CEXE_sources += AvgGalileanAlgorithm.cpp
 CEXE_sources += PMLPsatdAlgorithm.cpp
+CEXE_sources += ComovingPsatdAlgorithm.cpp
 
 ifeq ($(USE_RZ),TRUE)
   CEXE_sources += SpectralBaseAlgorithmRZ.cpp
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index 6e3338fa8..2626c6982 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -33,6 +33,7 @@ class 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=false,
                         const bool periodic_single_box=false,
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);
+            }
         }
     }