Implement averaged algo on staggered grids & merge spectral classes (#1544)

* Refactor and clean up some spectral classes

* Abort when current correction or Vay deposition are not implemented

* Implement general equations for averaged Galilean

* Allocate averaged MultiFabs also when aux_is_nodal=1 and do_nodal=0

* Allocate +ngextra guard cells also for averaged MultiFabs

* Make alias MultiFabs for averaged aux data

* With averaging, interpolate from avg_fp (not fp) to aux

* Fix some limits of the coefficients

* Fix bug causing NaNs in spectral coefficients

* Add 2D CI test with same analysis as nodal test

* Add 3D CI test with same analysis as nodal test

* Add limit that was not covered (knorm=0 && knorm_c!=0 && nu=0)

* Allocate T2_coef only if Galilean algorithm is used

* Allocate X4_coef only if Galilean algorithm is used

* Remove extra ghost cell from 'avg_fp' MultiFabs

7 files changed, 922 insertions, 1173 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp
deleted file mode 100644
index ea381235e..000000000
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp
+++ /dev/null
@@ -1,376 +0,0 @@
-#include "FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H"
-#include "Utils/WarpXConst.H"
-#include <cmath>
-
-using namespace amrex;
-
-/* \brief Initialize coefficients for the update equation */
-AvgGalileanAlgorithm::AvgGalileanAlgorithm(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_galilean,
-                         const Real dt)
-     // Initialize members of base class
-     : SpectralBaseAlgorithm( spectral_kspace, dm,
-                              norder_x, norder_y, norder_z, nodal )
-{
-    const BoxArray& ba = spectral_kspace.spectralspace_ba;
-
-    // Allocate the arrays of coefficients
-    C_coef = SpectralRealCoefficients(ba, dm, 1, 0);
-    S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0);
-
-    C1_coef = SpectralRealCoefficients(ba, dm, 1, 0);
-    S1_coef = SpectralRealCoefficients(ba, dm, 1, 0);
-    C3_coef = SpectralRealCoefficients(ba, dm, 1, 0);
-    S3_coef = SpectralRealCoefficients(ba, dm, 1, 0);
-
-    Psi1_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
-    Psi2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
-    Psi3_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
-
-    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);
-
-    A1_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
-    A2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
-
-    Rhoold_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
-    Rhonew_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
-    Jcoef_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
-
-    InitializeSpectralCoefficients(spectral_kspace, dm, v_galilean, dt);
-
-}
-
-void AvgGalileanAlgorithm::InitializeSpectralCoefficients(
-         const SpectralKSpace& spectral_kspace,
-         const amrex::DistributionMapping& dm,
-         const Array<Real, 3>& v_galilean,
-         const amrex::Real dt)
-{
-    const BoxArray& ba = spectral_kspace.spectralspace_ba;
-    // Fill them with the right values:
-    // Loop over boxes and allocate the corresponding coefficients
-    // for each box owned by the local MPI proc
-    for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){
-
-        const Box& bx = ba[mfi];
-
-        // Extract pointers for the k vectors
-        const Real* modified_kx = modified_kx_vec[mfi].dataPtr();
-#if (AMREX_SPACEDIM==3)
-        const Real* modified_ky = modified_ky_vec[mfi].dataPtr();
-#endif
-        const Real* modified_kz = modified_kz_vec[mfi].dataPtr();
-        // Extract arrays for the coefficients
-        Array4<Real> C = C_coef[mfi].array();
-        Array4<Real> S_ck = S_ck_coef[mfi].array();
-        Array4<Real> C1 = C1_coef[mfi].array();
-        Array4<Real> S1 = S1_coef[mfi].array();
-        Array4<Real> C3 = C3_coef[mfi].array();
-        Array4<Real> S3 = S3_coef[mfi].array();
-
-        Array4<Complex> Psi1 = Psi1_coef[mfi].array();
-        Array4<Complex> Psi2 = Psi2_coef[mfi].array();
-        Array4<Complex> Psi3 = Psi3_coef[mfi].array();
-        Array4<Complex> X1 = X1_coef[mfi].array();
-        Array4<Complex> X2 = X2_coef[mfi].array();
-        Array4<Complex> X3 = X3_coef[mfi].array();
-        Array4<Complex> X4 = X4_coef[mfi].array();
-        Array4<Complex> Theta2 = Theta2_coef[mfi].array();
-        Array4<Complex> A1 = A1_coef[mfi].array();
-        Array4<Complex> A2 = A2_coef[mfi].array();
-
-        Array4<Complex> CRhoold = Rhoold_coef[mfi].array();
-        Array4<Complex> CRhonew = Rhonew_coef[mfi].array();
-        Array4<Complex> Jcoef   = Jcoef_coef[mfi].array();
-        // Extract reals (for portability on GPU)
-        Real vx = v_galilean[0];
-#if (AMREX_SPACEDIM==3)
-        Real vy = v_galilean[1];
-#endif
-        Real vz = v_galilean[2];
-
-        // Loop over indices within one box
-        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) +
-#if (AMREX_SPACEDIM==3)
-                std::pow(modified_ky[j], 2) +
-                std::pow(modified_kz[k], 2));
-#else
-                std::pow(modified_kz[j], 2));
-#endif
-
-            // Calculate coefficients
-            constexpr Real c = PhysConst::c;
-            constexpr Real c2 = PhysConst::c*PhysConst::c;
-            constexpr Real ep0 = PhysConst::ep0;
-            const Complex I{0.,1.};
-            if (k_norm != 0){
-
-                C(i,j,k) = std::cos(c*k_norm*dt);
-                S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm);
-
-                C1(i,j,k) = std::cos(0.5_rt*c*k_norm*dt);
-                S1(i,j,k) = std::sin(0.5_rt*c*k_norm*dt);
-                C3(i,j,k) = std::cos(1.5_rt*c*k_norm*dt);
-                S3(i,j,k) = std::sin(1.5_rt*c*k_norm*dt);
-
-                // Calculate dot product with galilean velocity
-                const Real kv = modified_kx[i]*vx +
-#if (AMREX_SPACEDIM==3)
-                                 modified_ky[j]*vy +
-                                 modified_kz[k]*vz;
-#else
-                                 modified_kz[j]*vz;
-#endif
-
-                const Real nu = kv/(k_norm*c);
-                const Complex theta = amrex::exp( 0.5_rt*I*kv*dt );
-                const Complex theta_star = amrex::exp( -0.5_rt*I*kv*dt );
-                const Complex e_theta = amrex::exp( I*c*k_norm*dt );
-
-                Theta2(i,j,k) = theta*theta;
-
-                if ( (nu != 1.) && (nu != 0) ) {
-
-                    // Note: the coefficients X1, X2, X3 do not correspond
-                    // exactly to the original Galilean paper, but the
-                    // update equation have been modified accordingly so that
-                    // the expressions/ below (with the update equations)
-                    // are mathematically equivalent to those of the paper.
-                    Complex x1 = 1._rt/(1._rt-nu*nu) *
-                        (theta_star - C(i,j,k)*theta + I*kv*S_ck(i,j,k)*theta);
-
-                    Complex C_rho = I* c2 /( (1._rt-theta*theta) * ep0);
-
-                    Psi1(i,j,k) = theta * ((S1(i,j,k) + I*nu*C1(i,j,k))
-                                  - Theta2(i,j,k) * (S3(i,j,k) + I*nu*C3(i,j,k))) /(c*k_norm*dt * (nu*nu - 1._rt));
-                    Psi2(i,j,k) = theta * ((C1(i,j,k) - I*nu*S1(i,j,k))
-                                  - Theta2(i,j,k) * (C3(i,j,k) - I*nu*S3(i,j,k))) /(c2*k_norm*k_norm*dt * (nu*nu - 1._rt));
-                    Psi3(i,j,k) = I * theta * (1._rt - theta*theta) /(c*k_norm*dt*nu);
-
-                    A1(i,j,k) = (Psi1(i,j,k)  - 1._rt + I * kv*Psi2(i,j,k)    )/ (c2* k_norm*k_norm * (nu*nu - 1._rt));
-                    A2(i,j,k) = (Psi3(i,j,k) - Psi1(i,j,k)) / (c2*k_norm*k_norm);
-
-                    CRhoold(i,j,k) = C_rho * (theta*theta * A1(i,j,k) - A2(i,j,k));
-                    CRhonew(i,j,k) = C_rho * (A2(i,j,k) - A1(i,j,k));
-                    Jcoef(i,j,k) = (I*kv*A1(i,j,k) + Psi2(i,j,k))/ep0;
-                    // x1, above, is identical to the original paper
-                    X1(i,j,k) = theta*x1/(ep0*c*c*k_norm*k_norm);
-                    // The difference betwen X2 and X3 below, and those
-                    // from the original paper is the factor ep0*k_norm*k_norm
-                    X2(i,j,k) = (x1 - theta*(1._rt - C(i,j,k)))
-                                /(theta_star-theta)/(ep0*k_norm*k_norm);
-                    X3(i,j,k) = (x1 - theta_star*(1._rt - C(i,j,k)))
-                                /(theta_star-theta)/(ep0*k_norm*k_norm);
-                    X4(i,j,k) = I*kv*X1(i,j,k) - theta*theta*S_ck(i,j,k)/ep0;
-                }
-                if ( nu == 0) {
-                    X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0*c*c*k_norm*k_norm);
-                    X2(i,j,k) = (1._rt - S_ck(i,j,k)/dt) / (ep0*k_norm*k_norm);
-                    X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k)/dt) / (ep0*k_norm*k_norm);
-                    X4(i,j,k) = -S_ck(i,j,k)/ep0;
-
-                    Psi1(i,j,k) = (-S1(i,j,k) + S3(i,j,k)) / (c*k_norm*dt);
-                    Psi2(i,j,k) = (-C1(i,j,k) + C3(i,j,k)) / (c2*k_norm*k_norm*dt);
-                    Psi3(i,j,k) = 1._rt;
-                    A1(i,j,k) = (c*k_norm*dt + S1(i,j,k) - S3(i,j,k)) / (c*c2 * k_norm*k_norm*k_norm * dt);
-                    A2(i,j,k) =  (c*k_norm*dt + S1(i,j,k) - S3(i,j,k)) / (c*c2 * k_norm*k_norm*k_norm * dt);
-                    CRhoold(i,j,k) = 2._rt * I * S1(i,j,k)  * ( dt*C(i,j,k) - S_ck(i,j,k))
-                                    / (c*k_norm*k_norm*k_norm*dt*dt*ep0);
-                    CRhonew(i,j,k) =  - I * (c2* k_norm*k_norm * dt*dt - C1(i,j,k) + C3(i,j,k))
-                                    / (c2 * k_norm*k_norm*k_norm*k_norm * ep0 * dt*dt);
-                    Jcoef(i,j,k) = (-C1(i,j,k) + C3(i,j,k)) / (c2*ep0*k_norm*k_norm*dt);
-                }
-                if ( nu == 1.) {
-                    X1(i,j,k) = (1._rt - e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*c*c*ep0*k_norm*k_norm);
-                    X2(i,j,k) = (3._rt - 4._rt*e_theta + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*k_norm*k_norm*(1._rt- e_theta));
-                    X3(i,j,k) = (3._rt - 2._rt/e_theta - 2._rt*e_theta + e_theta*e_theta - 2._rt*I*c*k_norm*dt) / (4._rt*ep0*(e_theta - 1._rt)*k_norm*k_norm);
-                    X4(i,j,k) = I*(-1._rt + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*c*k_norm);
-                }
-
-            } else { // Handle k_norm = 0, by using the analytical limit
-              C(i,j,k) = 1._rt;
-              S_ck(i,j,k) = dt;
-              C1(i,j,k) = 1._rt;
-              S1(i,j,k) =  0._rt;
-              C3(i,j,k) = 1._rt;
-              S3(i,j,k) = 0._rt;
-
-              X1(i,j,k) = dt*dt/(2._rt * ep0);
-              X2(i,j,k) = c2*dt*dt/(6._rt * ep0);
-              X3(i,j,k) = - c2*dt*dt/(3._rt * ep0);
-              X4(i,j,k) = -dt/ep0;
-              Theta2(i,j,k) = 1._rt;
-
-              Psi1(i,j,k) = 1._rt;
-              Psi2(i,j,k) = -dt;
-              Psi3(i,j,k) = 1._rt;
-              A1(i,j,k) = 13._rt * dt*dt /24._rt;
-              A2(i,j,k) = 13._rt * dt*dt /24._rt;
-              CRhoold(i,j,k) = -I*c2 * dt*dt / (3._rt * ep0);
-              CRhonew(i,j,k) = -5._rt*I*c2 * dt*dt / (24._rt * ep0);
-              Jcoef(i,j,k) = -dt/ep0;
-            }
-
-        });
-    }
-}
-
-/* Advance the E and B field in spectral space (stored in `f`)
- * over one time step */
-void
-AvgGalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
-
-    // Loop over boxes
-    for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){
-
-        const Box& bx = f.fields[mfi].box();
-
-        // Extract arrays for the fields to be updated
-        Array4<Complex> fields = f.fields[mfi].array();
-        // Extract arrays for the coefficients
-        Array4<const Real> C_arr = C_coef[mfi].array();
-        Array4<const Real> S_ck_arr = S_ck_coef[mfi].array();
-        Array4<const Complex> X1_arr = X1_coef[mfi].array();
-        Array4<const Complex> X2_arr = X2_coef[mfi].array();
-        Array4<const Complex> X3_arr = X3_coef[mfi].array();
-        Array4<const Complex> X4_arr = X4_coef[mfi].array();
-        Array4<const Complex> Theta2_arr = Theta2_coef[mfi].array();
-        Array4<const Complex> Psi1_arr = Psi1_coef[mfi].array();
-        Array4<const Complex> Psi2_arr = Psi2_coef[mfi].array();
-
-        Array4<const Complex> A1_arr = A1_coef[mfi].array();
-        Array4<const Complex> Rhonew_arr = Rhonew_coef[mfi].array();
-        Array4<const Complex> Rhoold_arr = Rhoold_coef[mfi].array();
-        Array4<const Complex> Jcoef_arr =Jcoef_coef[mfi].array();
-        // Extract pointers for the k vectors
-        const Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr();
-#if (AMREX_SPACEDIM==3)
-        const Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr();
-#endif
-        const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr();
-
-        // Loop over indices within one box
-        ParallelFor(bx,
-        [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-        {
-            // Record old values of the fields to be updated
-            using Idx = SpectralAvgFieldIndex;
-
-            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);
-
-            // Shortcut 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 Real kx = modified_kx_arr[i];
-
-#if (AMREX_SPACEDIM==3)
-            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];
-#endif
-            constexpr Real c2 = PhysConst::c*PhysConst::c;
-            constexpr Real inv_ep0 = 1._rt/PhysConst::ep0;
-            constexpr Complex I = Complex{0,1};
-
-            const Real C = C_arr(i,j,k);
-            const 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);
-            const Complex T2 = Theta2_arr(i,j,k);
-
-            const Complex Psi1 = Psi1_arr(i,j,k);
-            const Complex Psi2 = Psi2_arr(i,j,k);
-            const Complex A1 = A1_arr(i,j,k);
-            const Complex CRhoold= Rhoold_arr(i,j,k);
-            const Complex CRhonew= Rhonew_arr(i,j,k);
-            const Complex Jcoef = Jcoef_arr(i,j,k);
-
-            //Update E (see the original Galilean article)
-            fields(i,j,k,Idx::Ex) = T2*C*Ex_old
-                        + T2*S_ck*c2*I*(ky*Bz_old - kz*By_old)
-                        + X4*Jx - I*(X2*rho_new - T2*X3*rho_old)*kx;
-            fields(i,j,k,Idx::Ey) = T2*C*Ey_old
-                        + T2*S_ck*c2*I*(kz*Bx_old - kx*Bz_old)
-                        + X4*Jy - I*(X2*rho_new - T2*X3*rho_old)*ky;
-            fields(i,j,k,Idx::Ez) = T2*C*Ez_old
-                        + T2*S_ck*c2*I*(kx*By_old - ky*Bx_old)
-                        + X4*Jz - I*(X2*rho_new - T2*X3*rho_old)*kz;
-
-            // Update B (see the original Galilean article)
-            // Note: here X1 is T2*x1/(ep0*c*c*k_norm*k_norm), where
-            // x1 has the same definition as in the original paper
-            fields(i,j,k,Idx::Bx) = T2*C*Bx_old
-                        - T2*S_ck*I*(ky*Ez_old - kz*Ey_old)
-                        +      X1*I*(ky*Jz     - kz*Jy);
-            fields(i,j,k,Idx::By) = T2*C*By_old
-                        - T2*S_ck*I*(kz*Ex_old - kx*Ez_old)
-                        +      X1*I*(kz*Jx     - kx*Jz);
-            fields(i,j,k,Idx::Bz) = T2*C*Bz_old
-                        - T2*S_ck*I*(kx*Ey_old - ky*Ex_old)
-                        +      X1*I*(kx*Jy     - ky*Jx);
-
-            //Update the averaged E,B fields in time on the interval [(n-1/2)dx, (n+1/2)dx]
-            fields(i,j,k,Idx::Ex_avg) = Psi1*Ex_old
-                          - Psi2*c2*I*(ky*Bz_old - kz*By_old)
-                          + Jcoef*Jx + ( CRhonew * rho_new +  CRhoold*rho_old )*kx;
-            fields(i,j,k,Idx::Ey_avg) =  Psi1*Ey_old
-                          - Psi2*c2*I*(kz*Bx_old - kx*Bz_old)
-                          + Jcoef*Jy +( CRhonew * rho_new +  CRhoold*rho_old )*ky;
-            fields(i,j,k,Idx::Ez_avg) =  Psi1*Ez_old
-                          - Psi2*c2*I*(kx*By_old - ky*Bx_old)
-                          + Jcoef*Jz     + ( CRhonew * rho_new +  CRhoold*rho_old )*kz;
-
-            fields(i,j,k,Idx::Bx_avg) =  Psi1*Bx_old
-                          + I*Psi2*(ky*Ez_old - kz*Ey_old)
-                          + A1*I*(ky*Jz     - kz*Jy)*inv_ep0;
-            fields(i,j,k,Idx::By_avg) =  Psi1*By_old
-                          + I*Psi2*(kz*Ex_old - kx*Ez_old)
-                          + A1*I*(kz*Jx     - kx*Jz)*inv_ep0;
-            fields(i,j,k,Idx::Bz_avg) =  Psi1*Bz_old
-                          + I*Psi2*(kx*Ey_old - ky*Ex_old)
-                          + A1*I*(kx*Jy     - ky*Jx)*inv_ep0;
-                        });
-    }
-}
-
-void
-AvgGalileanAlgorithm::CurrentCorrection (SpectralFieldData& /*field_data*/,
-                                         std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/,
-                                         const std::unique_ptr<amrex::MultiFab>& /*rho*/)
-{
-    amrex::Abort("Current correction not implemented for averaged Galilean PSATD");
-}
-
-void
-AvgGalileanAlgorithm::VayDeposition (SpectralFieldData& /*field_data*/,
-                                     std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/)
-{
-    amrex::Abort("Vay deposition not implemented for averaged Galilean PSATD");
-}
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
index ab8e439ed..f46fd9eb9 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
@@ -1,10 +1,8 @@
 target_sources(WarpX
   PRIVATE
-    GalileanAlgorithm.cpp
-    PMLPsatdAlgorithm.cpp
     PsatdAlgorithm.cpp
+    PMLPsatdAlgorithm.cpp
     SpectralBaseAlgorithm.cpp
-    AvgGalileanAlgorithm.cpp
     ComovingPsatdAlgorithm.cpp
 )
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
deleted file mode 100644
index ea5aa7ec0..000000000
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
+++ /dev/null
@@ -1,592 +0,0 @@
-#include "GalileanAlgorithm.H"
-#include "Utils/WarpXConst.H"
-
-#include <cmath>
-
-
-#if WARPX_USE_PSATD
-
-using namespace amrex;
-
-/* \brief Initialize coefficients for the update equation */
-GalileanAlgorithm::GalileanAlgorithm(const SpectralKSpace& spectral_kspace,
-                         const DistributionMapping& dm,
-                         const int norder_x, const int norder_y,
-                         const int norder_z, const bool nodal,
-                         const Array<Real, 3>& v_galilean,
-                         const Real dt,
-                         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)
-{
-    const BoxArray& ba = spectral_kspace.spectralspace_ba;
-
-    // Allocate the arrays of coefficients
-    C_coef = SpectralRealCoefficients(ba, dm, 1, 0);
-    S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0);
-    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);
-
-    InitializeSpectralCoefficients(spectral_kspace, dm, dt);
-}
-
-/* Advance the E and B field in spectral space (stored in `f`) over one time step */
-void
-GalileanAlgorithm::pushSpectralFields (SpectralFieldData& f) const
-{
-    const bool update_with_rho = m_update_with_rho;
-
-    // Loop over boxes
-    for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){
-
-        const Box& bx = f.fields[mfi].box();
-
-        // Extract arrays for the fields to be updated
-        Array4<Complex> fields = f.fields[mfi].array();
-
-        // Extract arrays for the coefficients
-        Array4<const Real> C_arr = C_coef[mfi].array();
-        Array4<const Real> S_ck_arr = S_ck_coef[mfi].array();
-        Array4<const Complex> X1_arr = X1_coef[mfi].array();
-        Array4<const Complex> X2_arr = X2_coef[mfi].array();
-        Array4<const Complex> X3_arr = X3_coef[mfi].array();
-        Array4<const Complex> X4_arr = X4_coef[mfi].array();
-        Array4<const Complex> Theta2_arr = Theta2_coef[mfi].array();
-
-        // Extract pointers for the k vectors
-        const Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr();
-#if (AMREX_SPACEDIM==3)
-        const Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr();
-#endif
-        const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr();
-
-        // Loop over indices within one box
-        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 Real kx = modified_kx_arr[i];
-#if (AMREX_SPACEDIM==3)
-            const Real ky = modified_ky_arr[j];
-            const Real kz = modified_kz_arr[k];
-#else
-            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;
-            constexpr Complex I = Complex{0._rt,1._rt};
-
-            // The definition of these coefficients is explained in more detail
-            // in the function InitializeSpectralCoefficients below
-            const Real C = C_arr(i,j,k);
-            const 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);
-            const Complex T2 = Theta2_arr(i,j,k);
-
-            // The equations in the following are the update equations for B and E,
-            // equations (11a) and (11b) of (Lehe et al, PRE 94, 2016), respectively,
-            // (or their rho-free formulation)
-
-            // Update E (equation (11b) or its rho-free formulation):
-            if (update_with_rho) {
-
-                // Ex
-                fields(i,j,k,Idx::Ex) = T2*C*Ex_old
-                            + T2*S_ck*c2*I*(ky*Bz_old - kz*By_old)
-                            + X4*Jx - I*(X2*rho_new - T2*X3*rho_old)*kx;
-                // Ey
-                fields(i,j,k,Idx::Ey) = T2*C*Ey_old
-                            + T2*S_ck*c2*I*(kz*Bx_old - kx*Bz_old)
-                            + X4*Jy - I*(X2*rho_new - T2*X3*rho_old)*ky;
-                // Ez
-                fields(i,j,k,Idx::Ez) = T2*C*Ez_old
-                            + T2*S_ck*c2*I*(kx*By_old - ky*Bx_old)
-                            + X4*Jz - I*(X2*rho_new - T2*X3*rho_old)*kz;
-            } else {
-
-                Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz;
-                Complex k_dot_E = kx * Ex_old + ky * Ey_old + kz * Ez_old;
-
-                // Ex
-                fields(i,j,k,Idx::Ex) = T2 * C * Ex_old + I * T2 * S_ck * c2 * (ky * Bz_old - kz * By_old)
-                    + X4 * Jx + X2 * k_dot_E * kx + X3 * k_dot_J * kx;
-                // Ey
-                fields(i,j,k,Idx::Ey) = T2 * C * Ey_old + I * T2 * S_ck * c2 * (kz * Bx_old - kx * Bz_old)
-                    + X4 * Jy + X2 * k_dot_E * ky + X3 * k_dot_J * ky;
-                // Ez
-                fields(i,j,k,Idx::Ez) = T2 * C * Ez_old + I * T2 * S_ck * c2 * (kx * By_old - ky * Bx_old)
-                    + X4 * Jz + X2 * k_dot_E * kz + X3 * k_dot_J * kz;
-            }
-
-            // Update B (equation (11a) with X1 rescaled by theta/(epsilon_0*c**2*k**2)):
-            // Bx
-            fields(i,j,k,Idx::Bx) = T2*C*Bx_old - T2*S_ck*I*(ky*Ez_old - kz*Ey_old) + X1*I*(ky*Jz - kz*Jy);
-            // By
-            fields(i,j,k,Idx::By) = T2*C*By_old - T2*S_ck*I*(kz*Ex_old - kx*Ez_old) + X1*I*(kz*Jx - kx*Jz);
-            // Bz
-            fields(i,j,k,Idx::Bz) = T2*C*Bz_old - T2*S_ck*I*(kx*Ey_old - ky*Ex_old) + X1*I*(kx*Jy - ky*Jx);
-        });
-    }
-}
-
-void GalileanAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& spectral_kspace,
-                                                        const amrex::DistributionMapping& dm,
-                                                        const amrex::Real dt)
-{
-    const bool update_with_rho = m_update_with_rho;
-
-    const BoxArray& ba = spectral_kspace.spectralspace_ba;
-
-    // Loop over boxes and allocate the corresponding coefficients for each box
-    for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi) {
-
-        const Box& bx = ba[mfi];
-
-        // Extract pointers for the k vectors
-        const Real* kx   = modified_kx_vec[mfi].dataPtr();
-        const Real* kx_c = modified_kx_vec_centered[mfi].dataPtr();
-#if (AMREX_SPACEDIM==3)
-        const Real* ky   = modified_ky_vec[mfi].dataPtr();
-        const Real* ky_c = modified_ky_vec_centered[mfi].dataPtr();
-#endif
-        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();
-        Array4<Real> S_ck = S_ck_coef[mfi].array();
-        Array4<Complex> X1 = X1_coef[mfi].array();
-        Array4<Complex> X2 = X2_coef[mfi].array();
-        Array4<Complex> X3 = X3_coef[mfi].array();
-        Array4<Complex> X4 = X4_coef[mfi].array();
-        Array4<Complex> T2 = Theta2_coef[mfi].array();
-
-        // Extract Galilean velocity
-        Real vx = m_v_galilean[0];
-#if (AMREX_SPACEDIM==3)
-        Real vy = m_v_galilean[1];
-#endif
-        Real vz = m_v_galilean[2];
-
-        // Loop over indices within one box
-        ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-        {
-            // Calculate norm of vector
-            const 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
-            // 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;
-            constexpr Real c2   = c*c;
-            constexpr Real ep0  = PhysConst::ep0;
-            constexpr Complex I = Complex{0._rt,1._rt};
-
-            // Auxiliary coefficients used when update_with_rho=false
-            const Real dt2 = dt * dt;
-            const Real dt3 = dt * dt2;
-            Complex X2_old, X3_old;
-
-            // 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)
-                ky_c[j]*vy + kz_c[k]*vz;
-#else
-                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 (knorm != 0. && knorm_c != 0.) {
-
-                // Auxiliary coefficients
-                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(om * dt);
-                S_ck(i,j,k) = std::sin(om * dt) / om;
-
-                // See equation (12b)
-                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 != om/om_c) && (nu != -om/om_c) && (nu != 0.) ) {
-
-                    // x1 is the coefficient chi_1 in equation (12c)
-                    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 * 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 * (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 * 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 * 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 * 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 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 * 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) = 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) = 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 * om * dt) / (4._rt * ep0 * om);
-                }
-
-                // 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 * 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) = 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) = 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 * 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 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;
-                S_ck(i,j,k) = dt;
-
-                // X1 multiplies i*(k \times J) in the update equation for B
-                X1(i,j,k) = dt2 / (2._rt * ep0);
-
-                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 * dt2 / (6._rt * ep0);
-                    X3(i,j,k) = - c2 * dt2 / (3._rt * ep0);
-                } 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 * dt2 * 0.5_rt;
-                    X3(i,j,k) = - c2 * dt3 / (6._rt * ep0);
-                }
-
-                // Coefficient multiplying J in update equation for E
-                X4(i,j,k) = -dt / ep0;
-
-                // Limit of exp(I*(k \dot v_gal)*dt)
-                T2(i,j,k) = 1._rt;
-            }
-        });
-    }
-}
-
-void
-GalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
-                                      std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
-                                      const std::unique_ptr<amrex::MultiFab>& rho) {
-    // Profiling
-    BL_PROFILE("GalileanAlgorithm::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 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_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_c = modified_kz_vec_centered[mfi].dataPtr();
-
-        // Local copy of member variables before GPU loop
-        const amrex::Real dt = m_dt;
-
-        // Galilean velocity
-        const amrex::Real vgx = m_v_galilean[0];
-        const amrex::Real vgy = m_v_galilean[1];
-        const amrex::Real vgz = m_v_galilean[2];
-
-        // Loop over indices within one box
-        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   = 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_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_c = 0._rt;
-            const     amrex::Real kz_c = modified_kz_arr_c[j];
-#endif
-            constexpr Complex I = Complex{0._rt,1._rt};
-
-            const amrex::Real k_norm = std::sqrt(kx * kx + ky * ky + kz * kz);
-
-            // Correct J
-            if (k_norm != 0._rt)
-            {
-                const Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz;
-                const amrex::Real k_dot_vg = kx_c * vgx + ky_c * vgy + kz_c * vgz;
-
-                if ( k_dot_vg != 0._rt ) {
-
-                    const Complex rho_old_mod = rho_old * amrex::exp(I * k_dot_vg * dt);
-                    const Complex den = 1._rt - amrex::exp(I * k_dot_vg * dt);
-
-                    fields(i,j,k,Idx::Jx) = Jx - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) * kx / (k_norm * k_norm);
-                    fields(i,j,k,Idx::Jy) = Jy - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) * ky / (k_norm * k_norm);
-                    fields(i,j,k,Idx::Jz) = Jz - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) * kz / (k_norm * k_norm);
-
-                } else {
-
-                    fields(i,j,k,Idx::Jx) = Jx - (k_dot_J - I * (rho_new - rho_old) / dt) * kx / (k_norm * k_norm);
-                    fields(i,j,k,Idx::Jy) = Jy - (k_dot_J - I * (rho_new - rho_old) / dt) * ky / (k_norm * k_norm);
-                    fields(i,j,k,Idx::Jz) = Jz - (k_dot_J - I * (rho_new - rho_old) / dt) * kz / (k_norm * k_norm);
-                }
-            }
-        });
-    }
-
-    // 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
-GalileanAlgorithm::VayDeposition (SpectralFieldData& /*field_data*/,
-                                  std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/)
-{
-    amrex::Abort("Vay deposition not implemented for Galilean PSATD");
-}
-
-#endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
index 71b97b8bc..11a10bb91 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
@@ -1,7 +1,5 @@
 CEXE_sources += SpectralBaseAlgorithm.cpp
 CEXE_sources += PsatdAlgorithm.cpp
-CEXE_sources += GalileanAlgorithm.cpp
-CEXE_sources += AvgGalileanAlgorithm.cpp
 CEXE_sources += PMLPsatdAlgorithm.cpp
 CEXE_sources += ComovingPsatdAlgorithm.cpp
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
index 4bfcd9fc9..1ca4aaf31 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
@@ -1,5 +1,4 @@
-/* Copyright 2019 Maxence Thevenet, Remi Lehe, Revathi Jambunathan, Edoardo Zoni
- *
+/* Copyright 2019
  *
  * This file is part of WarpX.
  *
@@ -10,32 +9,45 @@
 
 #include "SpectralBaseAlgorithm.H"
 
-
 #if WARPX_USE_PSATD
-
-/**
- * \brief Class that updates the field in spectral space
+/* \brief Class that updates the field in spectral space
  * and stores the coefficients of the corresponding update equation.
  */
 class PsatdAlgorithm : public SpectralBaseAlgorithm
 {
-
     public:
-        PsatdAlgorithm(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::Real dt,
-                         const bool update_with_rho);
-        // Redefine functions from base class
-        virtual void pushSpectralFields(SpectralFieldData& f) const override final;
-        virtual int getRequiredNumberOfFields() const override final {
-            return SpectralFieldIndex::n_fields;
-        }
 
-        void InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace,
-                                    const amrex::DistributionMapping& dm,
-                                    const amrex::Real dt);
+        // TODO Add Doxygen docs
+        PsatdAlgorithm (
+            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_galilean,
+            const amrex::Real dt,
+            const bool update_with_rho,
+            const bool time_averaging);
+
+        // TODO Add Doxygen docs
+        virtual void pushSpectralFields (SpectralFieldData& f) const override final;
+
+        // TODO Add Doxygen docs
+        virtual int getRequiredNumberOfFields () const override final
+        {
+            if (m_time_averaging) {
+                return SpectralAvgFieldIndex::n_fields;
+            } else {
+                return SpectralFieldIndex::n_fields;
+            }
+        };
+
+        // TODO Add Doxygen docs
+        void InitializeSpectralCoefficients (
+            const SpectralKSpace& spectral_kspace,
+            const amrex::DistributionMapping& dm,
+            const amrex::Real dt);
 
         /**
          * \brief Virtual function for current correction in Fourier space
@@ -49,9 +61,10 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
          *                           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;
+        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
@@ -64,14 +77,34 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
          * \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;
+        virtual void VayDeposition (
+            SpectralFieldData& field_data,
+            std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
 
     private:
-        SpectralRealCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef;
+
+        // These real and complex coefficients are always allocated
+        SpectralRealCoefficients C_coef, S_ck_coef;
+        SpectralComplexCoefficients T2_coef, X1_coef, X2_coef, X3_coef, X4_coef;
+
+        // These real and complex coefficients are allocated only with averaged Galilean PSATD
+        SpectralRealCoefficients C1_coef, C3_coef, S1_coef, S3_coef;
+        SpectralComplexCoefficients Psi1_coef, Psi2_coef, Psi3_coef, Psi4_coef,
+                                    A1_coef, A2_coef, Rhoold_coef, Rhonew_coef, Jcoef_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;
+
+        // Other member variables
+        amrex::Array<amrex::Real,3> m_v_galilean;
         amrex::Real m_dt;
         bool m_update_with_rho;
+        bool m_time_averaging;
+        bool m_is_galilean;
 };
-
 #endif // WARPX_USE_PSATD
 #endif // WARPX_PSATD_ALGORITHM_H_
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index ce0f466aa..0a278a54a 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -1,67 +1,132 @@
-/* Copyright 2019 Remi Lehe, Revathi Jambunathan, Edoardo Zoni
+/* Copyright 2019
  *
  * This file is part of WarpX.
  *
  * License: BSD-3-Clause-LBNL
  */
-#include "WarpX.H"
 #include "PsatdAlgorithm.H"
 #include "Utils/WarpXConst.H"
-#include "Utils/WarpXProfilerWrapper.H"
 
 #include <cmath>
 
-
 #if WARPX_USE_PSATD
+
 using namespace amrex;
 
-/**
- * \brief Constructor
- */
-PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace,
-                         const DistributionMapping& dm,
-                         const int norder_x, const int norder_y,
-                         const int norder_z, const bool nodal, const Real dt,
-                         const bool update_with_rho)
-    // Initialize members of base class
+PsatdAlgorithm::PsatdAlgorithm(
+    const SpectralKSpace& spectral_kspace,
+    const DistributionMapping& dm,
+    const int norder_x,
+    const int norder_y,
+    const int norder_z,
+    const bool nodal,
+    const Array<Real,3>& v_galilean,
+    const Real dt,
+    const bool update_with_rho,
+    const bool time_averaging)
+    // Initializer list
     : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal),
-      m_dt(dt),
-      m_update_with_rho(update_with_rho)
+    // 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),
+    m_time_averaging(time_averaging)
 {
     const BoxArray& ba = spectral_kspace.spectralspace_ba;
 
-    // Allocate the arrays of coefficients
+    m_is_galilean = (v_galilean[0] != 0.) || (v_galilean[1] != 0.) || (v_galilean[2] != 0.);
+
+    // Always allocate these coefficients
     C_coef = SpectralRealCoefficients(ba, dm, 1, 0);
     S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0);
-    X1_coef = SpectralRealCoefficients(ba, dm, 1, 0);
-    X2_coef = SpectralRealCoefficients(ba, dm, 1, 0);
-    X3_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+    X1_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+    X2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+    X3_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+
+    if (m_is_galilean)
+    {
+        X4_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+        T2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+    }
 
-    // Initialize coefficients for update equations
+    // Allocate these coefficients only with averaged Galilean PSATD
+    if (time_averaging)
+    {
+        C1_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+        S1_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+        C3_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+        S3_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+        Psi1_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+        Psi2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+        Psi3_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+        A1_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+        A2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+        Rhoold_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+        Rhonew_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+        Jcoef_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+    }
+
+    // Initialize coefficients
     InitializeSpectralCoefficients(spectral_kspace, dm, dt);
 }
 
-/**
- * \brief Advance E and B fields in spectral space (stored in `f`) over one time step
- */
 void
-PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
-
+PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
+{
     const bool update_with_rho = m_update_with_rho;
+    const bool time_averaging  = m_time_averaging;
+    const bool is_galilean     = m_is_galilean;
 
     // Loop over boxes
-    for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){
-
+    for (MFIter mfi(f.fields); mfi.isValid(); ++mfi)
+    {
         const Box& bx = f.fields[mfi].box();
 
         // Extract arrays for the fields to be updated
         Array4<Complex> fields = f.fields[mfi].array();
-        // Extract arrays for the coefficients
+
+        // These coefficients are always allocated
         Array4<const Real> C_arr = C_coef[mfi].array();
         Array4<const Real> S_ck_arr = S_ck_coef[mfi].array();
-        Array4<const Real> X1_arr = X1_coef[mfi].array();
-        Array4<const Real> X2_arr = X2_coef[mfi].array();
-        Array4<const Real> X3_arr = X3_coef[mfi].array();
+        Array4<const Complex> X1_arr = X1_coef[mfi].array();
+        Array4<const Complex> X2_arr = X2_coef[mfi].array();
+        Array4<const Complex> X3_arr = X3_coef[mfi].array();
+
+        Array4<const Complex> X4_arr;
+        Array4<const Complex> T2_arr;
+        if (is_galilean)
+        {
+            X4_arr = X4_coef[mfi].array();
+            T2_arr = T2_coef[mfi].array();
+        }
+
+        // These coefficients are allocated only with averaged Galilean PSATD
+        Array4<const Complex> Psi1_arr;
+        Array4<const Complex> Psi2_arr;
+        Array4<const Complex> A1_arr;
+        Array4<const Complex> Rhonew_arr;
+        Array4<const Complex> Rhoold_arr;
+        Array4<const Complex> Jcoef_arr;
+
+        if (time_averaging)
+        {
+            Psi1_arr = Psi1_coef[mfi].array();
+            Psi2_arr = Psi2_coef[mfi].array();
+            A1_arr = A1_coef[mfi].array();
+            Rhonew_arr = Rhonew_coef[mfi].array();
+            Rhoold_arr = Rhoold_coef[mfi].array();
+            Jcoef_arr = Jcoef_coef[mfi].array();
+        }
+
         // Extract pointers for the k vectors
         const Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr();
 #if (AMREX_SPACEDIM==3)
@@ -70,10 +135,12 @@ PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
         const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr();
 
         // Loop over indices within one box
-        ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             using Idx = SpectralFieldIndex;
+            using AvgIdx = SpectralAvgFieldIndex;
 
+            // Record old values of the fields to be updated
             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);
@@ -81,180 +148,783 @@ PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
             const Complex By_old = fields(i,j,k,Idx::By);
             const Complex Bz_old = fields(i,j,k,Idx::Bz);
 
-            // Shortcut for the values of J and rho
+            // 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
+            // k vector values
             const Real kx = modified_kx_arr[i];
 #if (AMREX_SPACEDIM==3)
             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
-            constexpr Real c2 = PhysConst::c*PhysConst::c;
-            constexpr Real inv_eps0 = 1.0_rt/PhysConst::ep0;
 
-            const Complex I = Complex{0,1};
+            // Physical constants and imaginary unit
+            constexpr Real c2 = PhysConst::c * PhysConst::c;
+            constexpr Real inv_ep0 = 1._rt / PhysConst::ep0;
+            constexpr Complex I = Complex{0._rt, 1._rt};
 
+            // These coefficients are initialized in the function InitializeSpectralCoefficients
             const Real C = C_arr(i,j,k);
             const Real S_ck = S_ck_arr(i,j,k);
-            const Real X1 = X1_arr(i,j,k);
-            const Real X2 = X2_arr(i,j,k);
-            const Real X3 = X3_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 = (is_galilean) ? X4_arr(i,j,k) : - S_ck * inv_ep0;
+            const Complex T2 = (is_galilean) ? T2_arr(i,j,k) : 1.0_rt;
 
-            // Update E (see WarpX online documentation: theory section)
+            // Update equations for E in the formulation with rho
+            // T2 = 1 always with standard PSATD (zero Galilean velocity)
 
-            if (update_with_rho) {
+            if (update_with_rho)
+            {
+                fields(i,j,k,Idx::Ex) = T2 * C * Ex_old
+                    + I * c2 * T2 * S_ck * (ky * Bz_old - kz * By_old)
+                    + X4 * Jx - I * (X2 * rho_new - T2 * X3 * rho_old) * kx;
+
+                fields(i,j,k,Idx::Ey) = T2 * C * Ey_old
+                    + I * c2 * T2 * S_ck * (kz * Bx_old - kx * Bz_old)
+                    + X4 * Jy - I * (X2 * rho_new - T2 * X3 * rho_old) * ky;
 
-                fields(i,j,k,Idx::Ex) = C*Ex_old + S_ck*(c2*I*(ky*Bz_old-kz*By_old)-inv_eps0*Jx)
-                                        - I*(X2*rho_new-X3*rho_old)*kx;
+                fields(i,j,k,Idx::Ez) = T2 * C * Ez_old
+                    + I * c2 * T2 * S_ck * (kx * By_old - ky * Bx_old)
+                    + X4 * Jz - I * (X2 * rho_new - T2 * X3 * rho_old) * kz;
+            }
 
-                fields(i,j,k,Idx::Ey) = C*Ey_old + S_ck*(c2*I*(kz*Bx_old-kx*Bz_old)-inv_eps0*Jy)
-                                        - I*(X2*rho_new-X3*rho_old)*ky;
+            // Update equations for E in the formulation without rho
+            // T2 = 1 always with standard PSATD (zero Galilean velocity)
 
-                fields(i,j,k,Idx::Ez) = C*Ez_old + S_ck*(c2*I*(kx*By_old-ky*Bx_old)-inv_eps0*Jz)
-                                        - I*(X2*rho_new-X3*rho_old)*kz;
-            } else {
+            else {
 
-                Complex k_dot_J = kx*Jx + ky*Jy + kz*Jz;
-                Complex k_dot_E = kx*Ex_old + ky*Ey_old + kz*Ez_old;
+                Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz;
+                Complex k_dot_E = kx * Ex_old + ky * Ey_old + kz * Ez_old;
 
-                fields(i,j,k,Idx::Ex) = C*Ex_old + S_ck*(c2*I*(ky*Bz_old-kz*By_old)-inv_eps0*Jx)
-                                        + X2*k_dot_E*kx + X3*inv_eps0*k_dot_J*kx;
+                fields(i,j,k,Idx::Ex) = T2 * C * Ex_old
+                    + I * c2 * T2 * S_ck * (ky * Bz_old - kz * By_old)
+                    + X4 * Jx + X2 * k_dot_E * kx + X3 * k_dot_J * kx;
 
-                fields(i,j,k,Idx::Ey) = C*Ey_old + S_ck*(c2*I*(kz*Bx_old-kx*Bz_old)-inv_eps0*Jy)
-                                        + X2*k_dot_E*ky + X3*inv_eps0*k_dot_J*ky;
+                fields(i,j,k,Idx::Ey) = T2 * C * Ey_old
+                    + I * c2 * T2 * S_ck * (kz * Bx_old - kx * Bz_old)
+                    + X4 * Jy + X2 * k_dot_E * ky + X3 * k_dot_J * ky;
 
-                fields(i,j,k,Idx::Ez) = C*Ez_old + S_ck*(c2*I*(kx*By_old-ky*Bx_old)-inv_eps0*Jz)
-                                        + X2*k_dot_E*kz + X3*inv_eps0*k_dot_J*kz;
+                fields(i,j,k,Idx::Ez) = T2 * C * Ez_old
+                    + I * c2 * T2 * S_ck * (kx * By_old - ky * Bx_old)
+                    + X4 * Jz + X2 * k_dot_E * kz + X3 * k_dot_J * kz;
             }
 
-            // Update B (see WarpX online documentation: theory section)
+            // Update equations for B
+            // T2 = 1 always with standard PSATD (zero Galilean velocity)
+
+            fields(i,j,k,Idx::Bx) = T2 * C * Bx_old - I * T2 * S_ck * (ky * Ez_old - kz * Ey_old)
+                + I * X1 * (ky * Jz - kz * Jy);
 
-            fields(i,j,k,Idx::Bx) = C*Bx_old - S_ck*I*(ky*Ez_old-kz*Ey_old) + X1*I*(ky*Jz-kz*Jy);
+            fields(i,j,k,Idx::By) = T2 * C * By_old - I * T2 * S_ck * (kz * Ex_old - kx * Ez_old)
+                + I * X1 * (kz * Jx - kx * Jz);
 
-            fields(i,j,k,Idx::By) = C*By_old - S_ck*I*(kz*Ex_old-kx*Ez_old) + X1*I*(kz*Jx-kx*Jz);
+            fields(i,j,k,Idx::Bz) = T2 * C * Bz_old - I * T2 * S_ck * (kx * Ey_old - ky * Ex_old)
+                + I * X1 * (kx * Jy - ky * Jx);
 
-            fields(i,j,k,Idx::Bz) = C*Bz_old - S_ck*I*(kx*Ey_old-ky*Ex_old) + X1*I*(kx*Jy-ky*Jx);
-        } );
+            // Additional update equations for averaged Galilean algorithm
+
+            if (time_averaging)
+            {
+                // These coefficients are initialized in the function InitializeSpectralCoefficients below
+                const Complex Psi1 = Psi1_arr(i,j,k);
+                const Complex Psi2 = Psi2_arr(i,j,k);
+                const Complex A1 = A1_arr(i,j,k);
+                const Complex CRhoold = Rhoold_arr(i,j,k);
+                const Complex CRhonew = Rhonew_arr(i,j,k);
+                const Complex Jcoef = Jcoef_arr(i,j,k);
+
+                fields(i,j,k,AvgIdx::Ex_avg) = Psi1 * Ex_old
+                    - I * c2 * Psi2 * (ky * Bz_old - kz * By_old)
+                    + Jcoef * Jx + (CRhonew * rho_new + CRhoold * rho_old) * kx;
+
+                fields(i,j,k,AvgIdx::Ey_avg) = Psi1 * Ey_old
+                    - I * c2 * Psi2 * (kz * Bx_old - kx * Bz_old)
+                    + Jcoef * Jy + (CRhonew * rho_new + CRhoold * rho_old) * ky;
+
+                fields(i,j,k,AvgIdx::Ez_avg) = Psi1 * Ez_old
+                    - I * c2 * Psi2 * (kx * By_old - ky * Bx_old)
+                    + Jcoef * Jz + (CRhonew * rho_new + CRhoold * rho_old) * kz;
+
+                fields(i,j,k,AvgIdx::Bx_avg) = Psi1 * Bx_old + I * Psi2 * (ky * Ez_old - kz * Ey_old)
+                    + I * inv_ep0 * A1 * (ky * Jz - kz * Jy);
+
+                fields(i,j,k,AvgIdx::By_avg) = Psi1 * By_old + I * Psi2 * (kz * Ex_old - kx * Ez_old)
+                    + I * inv_ep0 * A1 * (kz * Jx - kx * Jz);
+
+                fields(i,j,k,AvgIdx::Bz_avg) = Psi1 * Bz_old + I * Psi2 * (kx * Ey_old - ky * Ex_old)
+                    + I * inv_ep0 * A1 * (kx * Jy - ky * Jx);
+            }
+        });
     }
 }
 
-/**
- * \brief Initialize coefficients for update equations
- */
-void PsatdAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace,
-                                    const amrex::DistributionMapping& dm,
-                                    const amrex::Real dt)
+void PsatdAlgorithm::InitializeSpectralCoefficients (
+    const SpectralKSpace& spectral_kspace,
+    const amrex::DistributionMapping& dm,
+    const amrex::Real dt)
 {
     const bool update_with_rho = m_update_with_rho;
+    const bool time_averaging  = m_time_averaging;
+    const bool is_galilean     = m_is_galilean;
 
     const BoxArray& ba = spectral_kspace.spectralspace_ba;
 
-    // Loop over boxes and allocate the corresponding coefficients
-    // for each box owned by the local MPI proc
-    for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){
-
+    // Loop over boxes and allocate the corresponding coefficients for each box
+    for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi)
+    {
         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
+        // Coefficients always allocated
         Array4<Real> C = C_coef[mfi].array();
         Array4<Real> S_ck = S_ck_coef[mfi].array();
-        Array4<Real> X1 = X1_coef[mfi].array();
-        Array4<Real> X2 = X2_coef[mfi].array();
-        Array4<Real> X3 = X3_coef[mfi].array();
+        Array4<Complex> X1 = X1_coef[mfi].array();
+        Array4<Complex> X2 = X2_coef[mfi].array();
+        Array4<Complex> X3 = X3_coef[mfi].array();
+
+        Array4<Complex> X4;
+        Array4<Complex> T2;
+        if (is_galilean)
+        {
+            X4 = X4_coef[mfi].array();
+            T2 = T2_coef[mfi].array();
+        }
+
+        // Coefficients allocated only with averaged Galilean PSATD
+        Array4<Real> C1;
+        Array4<Real> S1;
+        Array4<Real> C3;
+        Array4<Real> S3;
+        Array4<Complex> Psi1;
+        Array4<Complex> Psi2;
+        Array4<Complex> Psi3;
+        Array4<Complex> A1;
+        Array4<Complex> A2;
+        Array4<Complex> CRhoold;
+        Array4<Complex> CRhonew;
+        Array4<Complex> Jcoef;
+
+        if (time_averaging)
+        {
+            C1 = C1_coef[mfi].array();
+            S1 = S1_coef[mfi].array();
+            C3 = C3_coef[mfi].array();
+            S3 = S3_coef[mfi].array();
+            Psi1 = Psi1_coef[mfi].array();
+            Psi2 = Psi2_coef[mfi].array();
+            Psi3 = Psi3_coef[mfi].array();
+            A1 = A1_coef[mfi].array();
+            A2 = A2_coef[mfi].array();
+            CRhoold = Rhoold_coef[mfi].array();
+            CRhonew = Rhonew_coef[mfi].array();
+            Jcoef = Jcoef_coef[mfi].array();
+        }
+
+        // Extract Galilean velocity
+        Real vx = m_v_galilean[0];
+#if (AMREX_SPACEDIM==3)
+        Real vy = m_v_galilean[1];
+#endif
+        Real vz = m_v_galilean[2];
 
         // Loop over indices within one box
-        ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        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) +
+            // Calculate norm of k vector
+            const 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
+            // Calculate norm of k vector (centered)
+            const Real knorm_c = std::sqrt(
+                std::pow(kx_c[i], 2) +
 #if (AMREX_SPACEDIM==3)
-                std::pow(modified_ky[j],2) + std::pow(modified_kz[k],2));
+                std::pow(ky_c[j], 2) + std::pow(kz_c[k], 2));
 #else
-                std::pow(modified_kz[j],2));
+                std::pow(kz_c[j], 2));
 #endif
-            // Calculate coefficients
+            // Physical constants and imaginary unit
             constexpr Real c = PhysConst::c;
-            constexpr Real eps0 = PhysConst::ep0;
-
-            if (k_norm != 0) {
-                C(i,j,k)    = std::cos(c*k_norm*dt);
-                S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm);
-                X1(i,j,k) = (1.0_rt-C(i,j,k))/(eps0*c*c*k_norm*k_norm);
-                if (update_with_rho) {
-                    X2(i,j,k) = (1.0_rt-S_ck(i,j,k)/dt)/(eps0*k_norm*k_norm);
-                    X3(i,j,k) = (C(i,j,k)-S_ck(i,j,k)/dt)/(eps0*k_norm*k_norm);
-                } else {
-                    X2(i,j,k) = (1.0_rt-C(i,j,k))/(k_norm*k_norm);
-                    X3(i,j,k) = (S_ck(i,j,k)-dt)/(k_norm*k_norm);
+            constexpr Real c2 = c*c;
+            constexpr Real ep0 = PhysConst::ep0;
+            constexpr Complex I = Complex{0._rt, 1._rt};
+
+            // Auxiliary coefficients used when update_with_rho=false
+            const Real dt2 = dt * dt;
+            const Real dt3 = dt * dt2;
+            Complex X2_old, X3_old;
+
+            // Calculate the dot product of the k vector with the Galilean velocity.
+            // This has to be computed always with the centered (that is, nodal) finite-order
+            // modified k vectors, to work correctly for both nodal and staggered simulations.
+            // kv = 0 always with standard PSATD (zero Galilean velocity).
+            const Real kv = kx_c[i]*vx +
+#if (AMREX_SPACEDIM==3)
+                ky_c[j]*vy + kz_c[k]*vz;
+#else
+                kz_c[j]*vz;
+#endif
+
+            // Note that:
+            // - X1 multiplies i*(k \times J) in the update equation for B
+            // - X2 multiplies rho_new if update_with_rho = 1 or (k \dot E)
+            //      if update_with_rho = 0 in the update equation for E
+            // - X3 multiplies rho_old if update_with_rho = 1 or (k \dot J)
+            //      if update_with_rho = 0 in the update equation for E
+            // - X4 multiplies J in the update equation for E
+
+            if (knorm != 0. && knorm_c != 0.)
+            {
+                // Auxiliary coefficients
+                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);
+
+                C   (i,j,k) = std::cos(om * dt);
+                S_ck(i,j,k) = std::sin(om * dt) / om;
+
+                if (time_averaging)
+                {
+                    C1(i,j,k) = std::cos(0.5_rt * om * dt);
+                    S1(i,j,k) = std::sin(0.5_rt * om * dt);
+                    C3(i,j,k) = std::cos(1.5_rt * om * dt);
+                    S3(i,j,k) = std::sin(1.5_rt * om * dt);
+                }
+
+                const Real nu = kv / om_c;
+                const Real nu2 = nu * nu;
+                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);
+
+                // T2 = 1 always with standard PSATD (zero Galilean velocity)
+                if (is_galilean)
+                {
+                    T2(i,j,k) = theta * theta;
+                }
+                const Complex T2_tmp = (is_galilean) ? T2(i,j,k) : 1.0_rt;
+
+                // nu = 0 always with standard PSATD (zero Galilean velocity): skip this block
+                if (nu != om/om_c && nu != -om/om_c && nu != 0.)
+                {
+                    // x1 is the coefficient chi_1 in equation (12c)
+                    Complex x1 = om2_c / (om2 - nu2 * om2_c)
+                        * (theta_star - theta * C(i,j,k) + I * nu * om_c * theta * S_ck(i,j,k));
+
+                    X1(i,j,k) = theta * x1 / (ep0 * om2_c);
+
+                    if (update_with_rho)
+                    {
+                        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 // update_with_rho = 0
+                    {
+                        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_tmp * (X2_old - X3_old) / (om2_c * om2);
+
+                        X3(i,j,k) = I * c2 * X2_old * (T2_tmp - 1._rt) / (ep0 * nu * om3_c * om2);
+                    }
+
+                    if (is_galilean)
+                    {
+                        X4(i,j,k) = I * nu * om_c * X1(i,j,k) - T2_tmp * S_ck(i,j,k) / ep0;
+                    }
+
+                    // Averaged Galilean algorithm
+                    if (time_averaging)
+                    {
+                        Complex C_rho = I * c2 / ((1._rt - T2_tmp) * ep0);
+
+                        Psi1(i,j,k) = theta * ((om * S1(i,j,k) + I * nu * om_c * C1(i,j,k))
+                            - T2_tmp * (om * S3(i,j,k) + I * nu * om_c * C3(i,j,k)))
+                            / (dt * (nu2 * om2_c - om2));
+
+                        Psi2(i,j,k) = theta * ((om * C1(i,j,k) - I * nu * om_c * S1(i,j,k))
+                            - T2_tmp * (om * C3(i,j,k) - I * nu * om_c * S3(i,j,k)))
+                            / (om * dt * (nu2 * om2_c - om2));
+
+                        Psi3(i,j,k) = I * theta * (1._rt - T2_tmp) / (nu * om_c * dt);
+
+                        A1(i,j,k) = (Psi1(i,j,k) - 1._rt + I * nu * om_c * Psi2(i,j,k))
+                            / (nu2 * om2_c - om2);
+
+                        A2(i,j,k) = (Psi3(i,j,k) - Psi1(i,j,k)) / om2;
+
+                        CRhoold(i,j,k) = C_rho * (T2_tmp * A1(i,j,k) - A2(i,j,k));
+
+                        CRhonew(i,j,k) = C_rho * (A2(i,j,k) - A1(i,j,k));
+
+                        Jcoef(i,j,k) = (I * nu * om_c * A1(i,j,k) + Psi2(i,j,k)) / ep0;
+                    }
+                }
+
+                // nu = 0 always with standard PSATD (zero Galilean velocity)
+                if (nu == 0.)
+                {
+                    X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2);
+
+                    if (update_with_rho)
+                    {
+                        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 // update_with_rho = 0
+                    {
+                        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);
+                    }
+
+                    if (is_galilean)
+                    {
+                        X4(i,j,k) = - S_ck(i,j,k) / ep0;
+                    }
+
+                    // Averaged Galilean algorithm
+                    if (time_averaging)
+                    {
+                        Psi1(i,j,k) = (S3(i,j,k) - S1(i,j,k)) / (om * dt);
+
+                        Psi2(i,j,k) = (C3(i,j,k) - C1(i,j,k)) / (om2 * dt);
+
+                        Psi3(i,j,k) = 1._rt;
+
+                        A1(i,j,k) = (om * dt + S1(i,j,k) - S3(i,j,k)) / (om3 * dt);
+
+                        A2(i,j,k) = (om * dt + S1(i,j,k) - S3(i,j,k)) / (om3 * dt);
+
+                        CRhoold(i,j,k) = 2._rt * I * c2 * S1(i,j,k) * (dt * C(i,j,k) - S_ck(i,j,k))
+                            / (om3 * dt2 * ep0);
+
+                        CRhonew(i,j,k) = - I * c2 * (om2 * dt2 - C1(i,j,k) + C3(i,j,k))
+                            / (om2 * om2 * dt2 * ep0);
+
+                        Jcoef(i,j,k) = (I * nu * om_c * A1(i,j,k) + Psi2(i,j,k)) / ep0;
+                    }
+                }
+
+                // nu = 0 always with standard PSATD (zero Galilean velocity): skip this block
+                if (nu == om/om_c)
+                {
+                    X1(i,j,k) = (1._rt - tmp1 * tmp1 + 2._rt * I * om * dt) / (4._rt * ep0 * om2);
+
+                    if (update_with_rho)
+                    {
+                        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 // update_with_rho = 0
+                    {
+                        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);
+                    }
+
+                    if (is_galilean)
+                    {
+                        X4(i,j,k) = (- I + I * tmp1 * tmp1 - 2._rt * om * dt) / (4._rt * ep0 * om);
+                    }
+
+                    // Averaged Galilean algorithm
+                    if (time_averaging)
+                    {
+                        Complex C_rho = I * c2 / ((1._rt - T2_tmp) * ep0);
+
+                        Psi1(i,j,k) = (2._rt * om * dt + I * tmp1 - I * tmp1 * tmp1 * tmp1)
+                            / (4._rt * om * dt);
+
+                        Psi2(i,j,k) = (- 2._rt * I * om * dt - tmp1 + tmp1 * tmp1 * tmp1)
+                            / (4._rt * om2 * dt);
+
+                        Psi3(i,j,k) = I * theta * (1._rt - T2_tmp) / (nu * om_c * dt);
+
+                        A1(i,j,k) = (2._rt * om * dt + I * (4._rt * om2 * dt2 - tmp1 + tmp1 * tmp1 * tmp1))
+                            / (8._rt * om3 * dt);
+
+                        A2(i,j,k) = (Psi3(i,j,k) - Psi1(i,j,k)) / om2;
+
+                        CRhoold(i,j,k) = C_rho * (T2_tmp * A1(i,j,k) - A2(i,j,k));
+
+                        CRhonew(i,j,k) = C_rho * (A2(i,j,k) - A1(i,j,k));
+
+                        Jcoef(i,j,k) = (I * nu * om_c * A1(i,j,k) + Psi2(i,j,k)) / ep0;
+                    }
+                }
+
+                // nu = 0 always with standard PSATD (zero Galilean velocity): skip this block
+                if (nu == -om/om_c)
+                {
+                    X1(i,j,k) = (1._rt - tmp2 * tmp2 - 2._rt * I * om * dt) / (4._rt * ep0 * om2);
+
+                    if (update_with_rho)
+                    {
+                        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 // update_with_rho = 0
+                    {
+                        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);
+                    }
+
+                    if (is_galilean)
+                    {
+                        X4(i,j,k) = (I - I * tmp2 * tmp2 - 2._rt * om * dt) / (4._rt * ep0 * om);
+                    }
+
+                    // Averaged Galilean algorithm
+                    if (time_averaging)
+                    {
+                        Complex C_rho = I * c2 / ((1._rt - T2_tmp) * ep0);
+
+                        Psi1(i,j,k) = (2._rt * om * dt - I * tmp2 + I * tmp2 * tmp2 * tmp2)
+                            / (4._rt * om * dt);
+
+                        Psi2(i,j,k) = (2._rt * I * om * dt - tmp2 + tmp2 * tmp2 * tmp2)
+                            / (4._rt * om2 * dt);
+
+                        Psi3(i,j,k) = I * theta * (1._rt - T2_tmp) / (nu * om_c * dt);
+
+                        A1(i,j,k) = (2._rt * om * dt * (1._rt - 2._rt * I * om * dt)
+                            + I * (tmp2 - tmp2 * tmp2 * tmp2)) / (8._rt * om3 * dt);
+
+                        A2(i,j,k) = (Psi3(i,j,k) - Psi1(i,j,k)) / om2;
+
+                        CRhoold(i,j,k) = C_rho * (T2_tmp * A1(i,j,k) - A2(i,j,k));
+
+                        CRhonew(i,j,k) = C_rho * (A2(i,j,k) - A1(i,j,k));
+
+                        Jcoef(i,j,k) = (I * nu * om_c * A1(i,j,k) + Psi2(i,j,k)) / ep0;
+                    }
+                }
+            }
+
+            else if (knorm != 0. && knorm_c == 0.)
+            {
+                const Real om = c * knorm;
+                const Real om2 = om * om;
+                const Real om3 = om * om2;
+
+                C(i,j,k) = std::cos(om * dt);
+
+                S_ck(i,j,k) = std::sin(om * dt) / om;
+
+                X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2);
+
+                if (update_with_rho)
+                {
+                    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 // update_with_rho = 0
+                {
+                    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);
+                }
+
+                if (is_galilean)
+                {
+                    X4(i,j,k) = - S_ck(i,j,k) / ep0;
+
+                    T2(i,j,k) = 1._rt;
                 }
-            } else { // Handle k_norm = 0 with analytical limit
-                C(i,j,k) = 1.0_rt;
+
+                // Averaged Galilean algorithm
+                if (time_averaging)
+                {
+                    C1(i,j,k) = std::cos(0.5_rt * om * dt);
+                    S1(i,j,k) = std::sin(0.5_rt * om * dt);
+                    C3(i,j,k) = std::cos(1.5_rt * om * dt);
+                    S3(i,j,k) = std::sin(1.5_rt * om * dt);
+
+                    Psi1(i,j,k) = (S3(i,j,k) - S1(i,j,k)) / (om * dt);
+
+                    Psi2(i,j,k) = (C3(i,j,k) - C1(i,j,k)) / (om2 * dt);
+
+                    Psi3(i,j,k) = 1._rt;
+
+                    A1(i,j,k) = (om * dt + S1(i,j,k) - S3(i,j,k)) / (om3 * dt);
+
+                    A2(i,j,k) = (om * dt + S1(i,j,k) - S3(i,j,k)) / (om3 * dt);
+
+                    CRhoold(i,j,k) = 2._rt * I * c2 * S1(i,j,k) * (dt * C(i,j,k) - S_ck(i,j,k))
+                        / (om3 * dt2 * ep0);
+
+                    CRhonew(i,j,k) = - I * c2 * (om2 * dt2 - C1(i,j,k) + C3(i,j,k))
+                        / (om2 * om2 * dt2 * ep0);
+
+                    Jcoef(i,j,k) = Psi2(i,j,k) / ep0;
+                }
+            }
+
+            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 Real nu2 = nu * nu;
+                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;
-                X1(i,j,k) = 0.5_rt*dt*dt/eps0;
-                if (update_with_rho) {
-                    X2(i,j,k) = c*c*dt*dt/(6.0_rt*eps0);
-                    X3(i,j,k) = -c*c*dt*dt/(3.0_rt*eps0);
-                } else {
-                    X2(i,j,k) = 0.5_rt*dt*dt*c*c;
-                    X3(i,j,k) = -c*c*dt*dt*dt/6.0_rt;
+
+                if (is_galilean)
+                {
+                    T2(i,j,k) = theta * theta;
+                }
+                const Complex T2_tmp = (is_galilean) ? T2(i,j,k) : 1.0_rt;
+
+                // nu = 0 always with standard PSATD (zero Galilean velocity): skip this block
+                if (nu != 0.)
+                {
+                    X1(i,j,k) = (- 1._rt + T2_tmp - I * nu * om_c * dt * T2_tmp) / (ep0 * nu2 * om2_c);
+
+                    if (update_with_rho)
+                    {
+                        X2(i,j,k) = c2 * (1._rt - T2_tmp + I * nu * om_c * dt * T2_tmp
+                            + 0.5_rt * nu2 * om2_c * dt2 * T2_tmp) / (ep0 * nu2 * om2_c * (T2_tmp - 1._rt));
+
+                        X3(i,j,k) = c2 * (1._rt - T2_tmp + I * nu * om_c * dt * T2_tmp
+                            + 0.5_rt * nu2 * om2_c * dt2) / (ep0 * nu2 * om2_c * (T2_tmp - 1._rt));
+                    }
+
+                    else // update_with_rho = 0
+                    {
+                        X2(i,j,k) = c2 * dt2 * T2_tmp * 0.5_rt;
+
+                        X3(i,j,k) = c2 * (2._rt * I - 2._rt * nu * om_c * dt * T2_tmp
+                            + I * nu2 * om2_c * dt2 * T2_tmp) / (2._rt * ep0 * nu2 * nu * om3_c);
+                    }
+
+                    if (is_galilean)
+                    {
+                        X4(i,j,k) = I * (T2_tmp - 1._rt) / (ep0 * nu * om_c);
+                    }
+
+                    // Averaged Galilean algorithm
+                    if (time_averaging)
+                    {
+                        Complex C_rho = I * c2 / ((1._rt - T2_tmp) * ep0);
+
+                        Psi1(i,j,k) = I * theta * (1._rt - T2_tmp) / (nu * om_c * dt);
+
+                        Psi2(i,j,k) = theta * (2._rt - I * nu * om_c * dt + T2_tmp * (3._rt * I * nu * om_c * dt - 2._rt))
+                            / (2._rt * nu2 * om2_c * dt);
+
+                        Psi3(i,j,k) = I * theta * (1._rt - T2_tmp) / (nu * om_c * dt);
+
+                        A1(i,j,k) = (Psi1(i,j,k) - 1._rt + I * nu * om_c * Psi2(i,j,k)) / (nu2 * om2_c);
+
+                        A2(i,j,k) = theta * (8._rt * I * (T2_tmp - 1._rt) + 4._rt * nu * om_c * dt
+                            * (3._rt * T2_tmp - 1._rt) + I * nu2 * om2_c * dt2 * (1._rt - 9._rt * T2_tmp))
+                            / (8._rt * nu2 * nu * om2_c * om_c * dt);
+
+                        CRhoold(i,j,k) = C_rho * (T2_tmp * A1(i,j,k) - A2(i,j,k));
+
+                        CRhonew(i,j,k) = C_rho * (A2(i,j,k) - A1(i,j,k));
+
+                        Jcoef(i,j,k) = (I * nu * om_c * A1(i,j,k) + Psi2(i,j,k)) / ep0;
+                    }
+                }
+
+                else // nu = 0
+                {
+                    X1(i,j,k) = dt2 / (2._rt * ep0);
+
+                    if (update_with_rho)
+                    {
+                        X2(i,j,k) = c2 * dt2 / (6._rt * ep0);
+
+                        X3(i,j,k) = - c2 * dt2 / (3._rt * ep0);
+                    }
+
+                    else // update_with_rho = 0
+                    {
+                        X2(i,j,k) = c2 * dt2 * 0.5_rt;
+
+                        X3(i,j,k) = - c2 * dt3 / (6._rt * ep0);
+                    }
+
+                    if (is_galilean)
+                    {
+                        X4(i,j,k) = - dt / ep0;
+                    }
+
+                    // Averaged Galilean algorithm
+                    if (time_averaging)
+                    {
+                        Psi1(i,j,k) = 1._rt;
+
+                        Psi2(i,j,k) = - dt;
+
+                        Psi3(i,j,k) = 1._rt;
+
+                        A1(i,j,k) = 13._rt * dt2 / 24._rt;
+
+                        A2(i,j,k) = 13._rt * dt2 / 24._rt;
+
+                        CRhoold(i,j,k) = - I * c2 * dt2 / (3._rt * ep0);
+
+                        CRhonew(i,j,k) = - 5._rt * I * c2 * dt2 / (24._rt * ep0);
+
+                        Jcoef(i,j,k) = - dt / ep0;
+                    }
+                }
+            }
+
+            else if (knorm == 0. && knorm_c == 0.)
+            {
+                C(i,j,k) = 1._rt;
+
+                S_ck(i,j,k) = dt;
+
+                X1(i,j,k) = dt2 / (2._rt * ep0);
+
+                if (update_with_rho)
+                {
+                    X2(i,j,k) = c2 * dt2 / (6._rt * ep0);
+
+                    X3(i,j,k) = - c2 * dt2 / (3._rt * ep0);
+                }
+
+                else // update_with_rho = 0
+                {
+                    X2(i,j,k) = c2 * dt2 * 0.5_rt;
+
+                    X3(i,j,k) = - c2 * dt3 / (6._rt * ep0);
+                }
+
+                // T2 = 1 always with standard PSATD (zero Galilean velocity)
+                if (is_galilean)
+                {
+                    X4(i,j,k) = - dt / ep0;
+
+                    T2(i,j,k) = 1._rt;
+                }
+
+                // Averaged Galilean algorithm
+                if (time_averaging)
+                {
+                    Psi1(i,j,k) = 1._rt;
+
+                    Psi2(i,j,k) = - dt;
+
+                    Psi3(i,j,k) = 1._rt;
+
+                    A1(i,j,k) = 13._rt * dt2 / 24._rt;
+
+                    A2(i,j,k) = 13._rt * dt2 / 24._rt;
+
+                    CRhoold(i,j,k) = - I * c2 * dt2 / (3._rt * ep0);
+
+                    CRhonew(i,j,k) = - 5._rt * I * c2 * dt2 / (24._rt * ep0);
+
+                    Jcoef(i,j,k) = - dt / ep0;
                 }
             }
-        } );
-     }
+        });
+    }
 }
 
 void
-PsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
-                                   std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
-                                   const std::unique_ptr<amrex::MultiFab>& rho) {
+PsatdAlgorithm::CurrentCorrection (
+    SpectralFieldData& field_data,
+    std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+    const std::unique_ptr<amrex::MultiFab>& rho)
+{
     // Profiling
-    WARPX_PROFILE("PsatdAlgorithm::CurrentCorrection()");
+    BL_PROFILE("PsatdAlgorithm::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 );
+    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 (MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
+    for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
 
-        const Box& bx = field_data.fields[mfi].box();
+        const amrex::Box& bx = field_data.fields[mfi].box();
 
         // Extract arrays for the fields to be updated
-        Array4<Complex> fields = field_data.fields[mfi].array();
+        amrex::Array4<Complex> fields = field_data.fields[mfi].array();
 
         // Extract pointers for the k vectors
-        const 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 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 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 Real dt = m_dt;
+        const amrex::Real dt = m_dt;
+
+        // Galilean velocity
+        const amrex::Real vgx = m_v_galilean[0];
+        const amrex::Real vgy = m_v_galilean[1];
+        const amrex::Real vgz = m_v_galilean[2];
 
         // Loop over indices within one box
-        ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        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);
@@ -264,42 +934,73 @@ PsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
             const Complex rho_new = fields(i,j,k,Idx::rho_new);
 
             // k vector values, and coefficients
-            const 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 Real ky = modified_ky_arr[j];
-            const 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 Real ky = 0;
-            const 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
-            const Real k_norm = std::sqrt( kx*kx + ky*ky + kz*kz );
-
-            constexpr Complex I = Complex{0,1};
+            constexpr Complex I = Complex{0._rt, 1._rt};
 
-            // div(J) in Fourier space
-            const Complex k_dot_J = kx*Jx + ky*Jy + kz*Jz;
+            const amrex::Real k_norm = std::sqrt(kx * kx + ky * ky + kz * kz);
 
             // Correct J
-            if ( k_norm != 0 )
+            if (k_norm != 0._rt)
             {
-                fields(i,j,k,Idx::Jx) = Jx - (k_dot_J-I*(rho_new-rho_old)/dt)*kx/(k_norm*k_norm);
-                fields(i,j,k,Idx::Jy) = Jy - (k_dot_J-I*(rho_new-rho_old)/dt)*ky/(k_norm*k_norm);
-                fields(i,j,k,Idx::Jz) = Jz - (k_dot_J-I*(rho_new-rho_old)/dt)*kz/(k_norm*k_norm);
+                const Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz;
+                const amrex::Real k_dot_vg = kx_c * vgx + ky_c * vgy + kz_c * vgz;
+
+                // k_dot_vg = 0 always with standard PSATD (zero Galilean velocity)
+                if ( k_dot_vg != 0._rt )
+                {
+                    const Complex rho_old_mod = rho_old * amrex::exp(I * k_dot_vg * dt);
+                    const Complex den = 1._rt - amrex::exp(I * k_dot_vg * dt);
+
+                    fields(i,j,k,Idx::Jx) = Jx - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den)
+                        * kx / (k_norm * k_norm);
+
+                    fields(i,j,k,Idx::Jy) = Jy - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den)
+                        * ky / (k_norm * k_norm);
+
+                    fields(i,j,k,Idx::Jz) = Jz - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den)
+                        * kz / (k_norm * k_norm);
+                }
+
+                else
+                {
+                    fields(i,j,k,Idx::Jx) = Jx - (k_dot_J - I * (rho_new - rho_old) / dt)
+                        * kx / (k_norm * k_norm);
+
+                    fields(i,j,k,Idx::Jy) = Jy - (k_dot_J - I * (rho_new - rho_old) / dt)
+                        * ky / (k_norm * k_norm);
+
+                    fields(i,j,k,Idx::Jz) = Jz - (k_dot_J - I * (rho_new - rho_old) / dt)
+                        * kz / (k_norm * k_norm);
+                }
             }
-        } );
+        });
     }
 
     // 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 );
+    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
-PsatdAlgorithm::VayDeposition (SpectralFieldData& field_data,
-                               std::array<std::unique_ptr<amrex::MultiFab>,3>& current) {
+PsatdAlgorithm::VayDeposition (
+    SpectralFieldData& field_data,
+    std::array<std::unique_ptr<amrex::MultiFab>,3>& current)
+{
     // Profiling
-    WARPX_PROFILE("PsatdAlgorithm::VayDeposition()");
+    BL_PROFILE("PsatdAlgorithm::VayDeposition()");
 
     using Idx = SpectralFieldIndex;
 
@@ -311,8 +1012,8 @@ PsatdAlgorithm::VayDeposition (SpectralFieldData& field_data,
     field_data.ForwardTransform(*current[2], Idx::Jz, 0, IntVect(1));
 
     // Loop over boxes
-    for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi) {
-
+    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
@@ -360,7 +1061,6 @@ PsatdAlgorithm::VayDeposition (SpectralFieldData& field_data,
             // Compute Jz
             if (kz_mod != 0._rt) fields(i,j,k,Idx::Jz) = I * Dz / kz_mod;
             else                 fields(i,j,k,Idx::Jz) = 0._rt;
-
         });
     }
 
@@ -369,4 +1069,5 @@ PsatdAlgorithm::VayDeposition (SpectralFieldData& field_data,
     field_data.BackwardTransform(*current[1], Idx::Jy, 0);
     field_data.BackwardTransform(*current[2], Idx::Jz, 0);
 }
+
 #endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index 4c9cb6967..725cfcf92 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -7,8 +7,6 @@
 #include "SpectralKSpace.H"
 #include "SpectralSolver.H"
 #include "SpectralAlgorithms/PsatdAlgorithm.H"
-#include "SpectralAlgorithms/GalileanAlgorithm.H"
-#include "SpectralAlgorithms/AvgGalileanAlgorithm.H"
 #include "SpectralAlgorithms/PMLPsatdAlgorithm.H"
 #include "SpectralAlgorithms/ComovingPsatdAlgorithm.H"
 #include "WarpX.H"
@@ -61,26 +59,15 @@ SpectralSolver::SpectralSolver(
             k_space, dm, norder_x, norder_y, norder_z, nodal, dt);
     }
     else {
-        if (fft_do_time_averaging){
-            algorithm = std::make_unique<AvgGalileanAlgorithm>(
-                k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt);
+        // Comoving PSATD algorithm
+        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);
         }
+        // PSATD algorithms: standard, Galilean, or averaged Galilean
         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);
-            }
+            algorithm = std::make_unique<PsatdAlgorithm>(
+                k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt, update_with_rho, fft_do_time_averaging);
         }
     }