Galilean PSATD: current correction and rho-free formulation (#1151)

* Introduce option to update E with/without rho

* Clean up

* Implement current correction for Galilean PSATD (needs bug fix)

* Include equations in docs

* Fix EOL whitespaces error

* Small clean-up

* Implement Galilean PSATD update without rho

* Clean up

* Fix bug in current correction

* Fix EOL whitespaces

* Clean up

* Fix unused import

* Remove unused variable

* [skip CI] Improve docs

* Clean up style

* Fix EOL whitespaces

* Fix EOL whitespaces

* Clean up style

* Revert analysis script to old status

* [skip CI] Clean up style

* Make equations more human-readable and improve comments

* 2D test with current correction works

* Temporary build fix as in #1197

* 3D test with current correction works

* Rename th and th_star as theta and theta_star

* Fix a couple of wrong comments

* Add vertical spaces to improve readability

* Improve documentation

* Function CurrentCorrection is now pure

* 2D benchmark fields data are now correct

* Add limits of coefficients for nu=-1

* Change default of update_with_rho for Galilean PSATD

12 files changed, 456 insertions, 162 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H
index 45a35098b..9d3facf3f 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H
@@ -27,6 +27,22 @@ class AvgGalileanAlgorithm : public SpectralBaseAlgorithm
            const amrex::Real dt);
 
         /**
+         * \brief Virtual function for current correction in Fourier space
+         * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>).
+         * This function overrides the virtual function \c CurrentCorrection in the
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
+         *
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
+         */
+        virtual void CurrentCorrection (SpectralFieldData& field_data,
+                                        std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+                                        const std::unique_ptr<amrex::MultiFab>& rho) override final;
+
+        /**
          * \brief Virtual function for Vay current deposition in Fourier space
          * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>).
          * This function overrides the virtual function \c VayDeposition in the
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp
index 0ba603031..e5fbc8261 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp
@@ -371,6 +371,14 @@ AvgGalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
 };
 
 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)
 {
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
index 799eca8b1..83ae83d69 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
@@ -10,21 +10,37 @@
 class GalileanAlgorithm : public SpectralBaseAlgorithm
 {
     public:
-        GalileanAlgorithm(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);
+        GalileanAlgorithm (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);
         // Redefine update equation from base class
-        virtual void pushSpectralFields(SpectralFieldData& f) const override final;
-        virtual int getRequiredNumberOfFields() const override final {
+        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::Array<amrex::Real, 3>& v_galilean,
-                                    const amrex::Real dt);
+        void InitializeSpectralCoefficients (const SpectralKSpace& spectral_kspace,
+                                             const amrex::DistributionMapping& dm,
+                                             const amrex::Real dt);
+
+        /**
+         * \brief Virtual function for current correction in Fourier space
+         * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>).
+         * This function overrides the virtual function \c CurrentCorrection in the
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
+         *
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
+         */
+        virtual void CurrentCorrection (SpectralFieldData& field_data,
+                                        std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+                                        const std::unique_ptr<amrex::MultiFab>& rho) override final;
 
         /**
          * \brief Virtual function for Vay current deposition in Fourier space
@@ -43,6 +59,9 @@ class GalileanAlgorithm : public SpectralBaseAlgorithm
     private:
         SpectralRealCoefficients C_coef, S_ck_coef;
         SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef;
+        amrex::Array<amrex::Real,3> m_v_galilean;
+        amrex::Real m_dt;
+        bool m_update_with_rho;
 };
 #endif // WARPX_USE_PSATD
 #endif // WARPX_GALILEAN_ALGORITHM_H_
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
index e606b3232..e969fd6cf 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
@@ -14,10 +14,13 @@ GalileanAlgorithm::GalileanAlgorithm(const SpectralKSpace& spectral_kspace,
                          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 Real dt,
+                         const bool update_with_rho)
      // Initialize members of base class
-     : SpectralBaseAlgorithm( spectral_kspace, dm,
-                              norder_x, norder_y, norder_z, nodal )
+     : m_v_galilean(v_galilean),
+       m_dt(dt),
+       m_update_with_rho(update_with_rho),
+       SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal)
 {
     const BoxArray& ba = spectral_kspace.spectralspace_ba;
 
@@ -30,14 +33,14 @@ GalileanAlgorithm::GalileanAlgorithm(const SpectralKSpace& spectral_kspace,
     X4_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
     Theta2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
 
-    InitializeSpectralCoefficients(spectral_kspace, dm, v_galilean, dt);
-
+    InitializeSpectralCoefficients(spectral_kspace, dm, dt);
 };
 
-/* Advance the E and B field in spectral space (stored in `f`)
- * over one time step */
+/* Advance the E and B field in spectral space (stored in `f`) over one time step */
 void
-GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
+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){
@@ -46,6 +49,7 @@ GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
 
         // 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();
@@ -63,8 +67,7 @@ GalileanAlgorithm::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
         {
             // Record old values of the fields to be updated
             using Idx = SpectralFieldIndex;
@@ -74,13 +77,15 @@ GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
             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
+
+            // 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];
@@ -89,8 +94,12 @@ GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
             constexpr Real ky = 0;
             const Real kz = modified_kz_arr[j];
 #endif
-            constexpr Real c2 = PhysConst::c*PhysConst::c;
-            constexpr Complex I = Complex{0,1};
+            // 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);
@@ -99,43 +108,62 @@ GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
             const Complex X4 = X4_arr(i,j,k);
             const Complex T2 = Theta2_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);
+            // 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 Array<Real, 3>& v_galilean,
-                                    const amrex::Real dt)
+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;
-    // 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){
+
+    // Loop over boxes and allocate the corresponding coefficients for each box
+    for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi) {
 
         const Box& bx = ba[mfi];
 
@@ -145,6 +173,7 @@ void GalileanAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spe
         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();
@@ -152,17 +181,17 @@ void GalileanAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spe
         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();
-        // Extract reals (for portability on GPU)
-        Real vx = v_galilean[0];
+        Array4<Complex> T2 = Theta2_coef[mfi].array();
+
+        // Extract Galilean velocity
+        Real vx = m_v_galilean[0];
 #if (AMREX_SPACEDIM==3)
-        Real vy = v_galilean[1];
+        Real vy = m_v_galilean[1];
 #endif
-        Real vz = v_galilean[2];
+        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(
@@ -173,75 +202,277 @@ void GalileanAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spe
 #else
                 std::pow(modified_kz[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};
 
-            // Calculate coefficients
-            constexpr Real c = PhysConst::c;
-            constexpr Real ep0 = PhysConst::ep0;
-            const Complex I{0.,1.};
-            if (k_norm != 0){
+            // Auxiliary coefficients used when update_with_rho=false
+            const Real dt2 = dt * dt;
+            const Real dt3 = dt * dt2;
+            Complex X2_old, X3_old;
 
-                C(i,j,k) = std::cos(c*k_norm*dt);
-                S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm);
-
-                // Calculate dot product with galilean velocity
-                const Real kv = modified_kx[i]*vx +
+            // Calculate dot product of k vector with Galilean velocity
+            const Real kv = modified_kx[i]*vx +
 #if (AMREX_SPACEDIM==3)
-                                 modified_ky[j]*vy +
-                                 modified_kz[k]*vz;
+                modified_ky[j]*vy + modified_kz[k]*vz;
 #else
-                                 modified_kz[j]*vz;
+                modified_kz[j]*vz;
 #endif
+            // The coefficients in the following refer to the ones given in equations
+            // (12a)-(12d) of (Lehe et al, PRE 94, 2016), used to update B and E
+            // (equations (11a) and (11b) of the same reference, respectively)
+
+            if (k_norm != 0.) {
+
+                // Auxiliary coefficients
+                const Real    k2   = k_norm * k_norm;
+                const Real    ck   = c * k_norm;
+                const Real    ckdt = ck * dt;
+                const Complex tmp1 = amrex::exp(  I * ckdt); // limit of T2 for nu =  1
+                const Complex tmp2 = amrex::exp(- I * ckdt); // limit of T2 for nu = -1
+
+                // See equation (12a)
+                C   (i,j,k) = std::cos(ckdt);
+                S_ck(i,j,k) = std::sin(ckdt) / ck;
+
+                // See equation (12b)
+                const Real    nu         = kv / ck;
+                const Complex theta      = amrex::exp(  I * 0.5_rt * kv * dt);
+                const Complex theta_star = amrex::exp(- I * 0.5_rt * kv * dt);
 
-                const Real nu = kv/(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);
-                    // 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;
+                // This is exp(i*(k \dot v_gal)*dt)
+                T2(i,j,k) = theta * theta;
+
+                if ( (nu != 1.) && (nu != 0.) ) {
+
+                    // x1 is the coefficient chi_1 in equation (12c)
+                    Complex x1 = 1._rt / (1._rt - nu*nu)
+                        * (theta_star - C(i,j,k) * theta + I * kv * S_ck(i,j,k) * theta);
+
+                    // X1 multiplies i*(k \times J) in the update equation for B
+                    X1(i,j,k) = theta * x1 / (ep0 * c2 * k2);
+
+                    if (update_with_rho) {
+                        // X2 multiplies rho_new in the update equation for E
+                        // X3 multiplies rho_old in the update equation for E
+                        X2(i,j,k) = (x1 - theta * (1._rt - C(i,j,k))) / (theta_star - theta) / (ep0 * k2);
+                        X3(i,j,k) = (x1 - theta_star * (1._rt - C(i,j,k))) / (theta_star - theta) / (ep0 * k2);
+                    } else {
+                        // X2_old is the coefficient chi_2 in equation (12d)
+                        // X3_old is the coefficient chi_3 in equation (12d)
+                        // X2 multiplies (k \dot E) in the update equation for E
+                        // X3 multiplies (k \dot J) in the update equation for E
+                        X2_old = (x1 - theta * (1._rt - C(i,j,k))) / (theta_star - theta);
+                        X3_old = (x1 - theta_star * (1._rt - C(i,j,k))) / (theta_star - theta);
+                        X2(i,j,k) = T2(i,j,k) * (X2_old - X3_old) / k2;
+                        X3(i,j,k) = I * X2_old * (T2(i,j,k) - 1._rt) / (ep0 * k2 * kv);
+                    }
+
+                    // X4 multiplies J in the update equation for E
+                    X4(i,j,k) = I * kv * X1(i,j,k) - T2(i,j,k) * S_ck(i,j,k) / ep0;
                 }
-                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;
+
+                // Limits for nu = 0
+                if (nu == 0.) {
+
+                    // X1 multiplies i*(k \times J) in the update equation for B
+                    X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * c2 * k2);
+
+                    if (update_with_rho) {
+                        // X2 multiplies rho_new in the update equation for E
+                        // X3 multiplies rho_old in the update equation for E
+                        X2(i,j,k) = (1._rt - S_ck(i,j,k) / dt) / (ep0 * k2);
+                        X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * k2);
+                    } else {
+                        // X2 multiplies (k \dot E) in the update equation for E
+                        // X3 multiplies (k \dot J) in the update equation for E
+                        X2(i,j,k) = (1._rt - C(i,j,k)) / k2;
+                        X3(i,j,k) = (S_ck(i,j,k) / dt - 1._rt) * dt / (ep0 * k2);
+                    }
+
+                    // Coefficient multiplying J in update equation for E
+                    X4(i,j,k) = - S_ck(i,j,k) / ep0;
                 }
-                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);
+
+                // Limits for nu = 1
+                if (nu == 1.) {
+
+                    // X1 multiplies i*(k \times J) in the update equation for B
+                    X1(i,j,k) = (1._rt - tmp1 * tmp1 + 2._rt * I * ckdt) / (4._rt * ep0 * c2 * k2);
+
+                    if (update_with_rho) {
+                        // X2 multiplies rho_new in the update equation for E
+                        // X3 multiplies rho_old in the update equation for E
+                        X2(i,j,k) = (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * ckdt)
+                            / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
+                        X3(i,j,k) = (3._rt - 2._rt / tmp1 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * ckdt)
+                            / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
+                    } else {
+                        // X2 multiplies (k \dot E) in the update equation for E
+                        // X3 multiplies (k \dot J) in the update equation for E
+                        X2(i,j,k) = (1._rt - C(i,j,k)) * tmp1 / k2;
+                        X3(i,j,k) = (2._rt * ckdt - I * tmp1 * tmp1 + 4._rt * I * tmp1 - 3._rt * I)
+                            / (4._rt * ep0 * ck * k2);
+                    }
+
+                    // Coefficient multiplying J in update equation for E
+                    X4(i,j,k) = (- I + I * tmp1 * tmp1 - 2._rt * ckdt) / (4._rt * ep0 * ck);
                 }
 
-            } else { // Handle k_norm = 0, by using the analytical limit
+                // Limits for nu = -1
+                if (nu == -1.) {
+
+                    // X1 multiplies i*(k \times J) in the update equation for B
+                    X1(i,j,k) = (1._rt - tmp2 * tmp2 - 2._rt * I * ckdt) / (4._rt * ep0 * c2 * k2);
+
+                    if (update_with_rho) {
+                        // X2 multiplies rho_new in the update equation for E
+                        // X3 multiplies rho_old in the update equation for E
+                        X2(i,j,k) = (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * ckdt * tmp1)
+                            / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
+                        X3(i,j,k) = (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * ckdt * tmp1)
+                            / (4._rt * ep0 * k2 * (tmp1 - 1._rt));
+                    } else {
+                        // X2 multiplies (k \dot E) in the update equation for E
+                        // X3 multiplies (k \dot J) in the update equation for E
+                        X2(i,j,k) = (1._rt - C(i,j,k)) * tmp2 / k2;
+                        X3(i,j,k) = (2._rt * ckdt + I * tmp2 * tmp2 - 4._rt * I * tmp2 + 3._rt * I)
+                            / (4._rt * ep0 * ck * k2);
+                    }
+
+                    // Coefficient multiplying J in update equation for E
+                    X4(i,j,k) = (I - I * tmp2 * tmp2 - 2._rt * ckdt) / (4._rt * ep0 * ck);
+                }
+            }
+
+            // Limits for k = 0
+            else {
+
+                // 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(i,j,k) = dt*dt/(2._rt * ep0);
-                X2(i,j,k) = c*c*dt*dt/(6._rt * ep0);
-                X3(i,j,k) = - c*c*dt*dt/(3._rt * ep0);
-                X4(i,j,k) = -dt/ep0;
-                Theta2(i,j,k) = 1._rt;
+
+                // X1 multiplies i*(k \times J) in the update equation for B
+                X1(i,j,k) = dt2 / (2._rt * ep0);
+
+                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();
+#if (AMREX_SPACEDIM==3)
+        const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+#endif
+        const amrex::Real* const modified_kz_arr = modified_kz_vec[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
+        {
+            // Record old values of the fields to be updated
+            using Idx = SpectralFieldIndex;
+
+            // 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];
+#if (AMREX_SPACEDIM==3)
+            const amrex::Real ky = modified_ky_arr[j];
+            const amrex::Real kz = modified_kz_arr[k];
+#else
+            constexpr amrex::Real ky = 0._rt;
+            const amrex::Real kz = modified_kz_arr[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 * vgx + ky * vgy + kz * 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
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
index 347012d5a..348282dce 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
@@ -35,6 +35,22 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm
         }
 
         /**
+         * \brief Virtual function for current correction in Fourier space
+         * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>).
+         * This function overrides the virtual function \c CurrentCorrection in the
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
+         *
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
+         */
+        virtual void CurrentCorrection (SpectralFieldData& field_data,
+                                        std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+                                        const std::unique_ptr<amrex::MultiFab>& rho) override final;
+
+        /**
          * \brief Virtual function for Vay current deposition in Fourier space
          * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>).
          * This function overrides the virtual function \c VayDeposition in the
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
index d2f087706..4bc147cd8 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
@@ -156,6 +156,14 @@ void PMLPsatdAlgorithm::InitializeSpectralCoefficients (
 };
 
 void
+PMLPsatdAlgorithm::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 PML PSATD");
+}
+
+void
 PMLPsatdAlgorithm::VayDeposition (SpectralFieldData& field_data,
                                   std::array<std::unique_ptr<amrex::MultiFab>,3>& current)
 {
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
index bfa9283e6..4bfcd9fc9 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
@@ -39,19 +39,19 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
 
         /**
          * \brief Virtual function for current correction in Fourier space
-         * (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010).
+         * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>).
          * This function overrides the virtual function \c CurrentCorrection in the
-         * base class \c SpectralBaseAlgorithm (qualifier \c override) and cannot be
-         * overridden by further derived classes (qualifier \c final).
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
          *
-         * \param[in,out] field_data all fields in Fourier space
-         * \param[in,out] current    three-dimensional array of unique pointers to MultiFab
-         *                           storing the three components of the current density
-         * \param[in]     rho        unique pointer to MultiFab storing the charge density
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
          */
-        virtual void CurrentCorrection ( SpectralFieldData& field_data,
-                                         std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
-                                         const std::unique_ptr<amrex::MultiFab>& rho ) override final;
+        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
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index 7f9fd3edb..89138641f 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -220,9 +220,9 @@ void PsatdAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spectr
 }
 
 void
-PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data,
+PsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
                                    std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
-                                   const std::unique_ptr<amrex::MultiFab>& rho ) {
+                                   const std::unique_ptr<amrex::MultiFab>& rho) {
     // Profiling
     WARPX_PROFILE( "PsatdAlgorithm::CurrentCorrection" );
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
index fc03f95f2..b4c2597a8 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
@@ -30,19 +30,19 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ
 
         /**
          * \brief Virtual function for current correction in Fourier space
-         * (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010).
+         * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>).
          * This function overrides the virtual function \c CurrentCorrection in the
-         * base class \c SpectralBaseAlgorithmRZ (qualifier \c override) and cannot be
-         * overridden by further derived classes (qualifier \c final).
+         * base class \c SpectralBaseAlgorithmRZ and cannot be overridden by further
+         * derived classes.
          *
-         * \param[in,out] field_data all fields in Fourier space
-         * \param[in,out] current    two-dimensional array of unique pointers to MultiFab
-         *                           storing the three components of the current density
-         * \param[in]     rho        unique pointer to MultiFab storing the charge density
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
          */
-        virtual void CurrentCorrection ( SpectralFieldDataRZ& field_data,
-                                         std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
-                                         const std::unique_ptr<amrex::MultiFab>& rho ) override final;
+        virtual void CurrentCorrection (SpectralFieldDataRZ& 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
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
index 35a478f2c..349f266d1 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
@@ -200,9 +200,9 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const
 void
 PsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data,
                                      std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
-                                     const std::unique_ptr<amrex::MultiFab>& rho )
+                                     const std::unique_ptr<amrex::MultiFab>& rho)
 {
-    amrex::Abort("Current correction not implemented in RZ");
+    amrex::Abort("Current correction not implemented in RZ geometry");
 }
 
 void
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
index 2a34d21ab..73ad773f3 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
@@ -33,19 +33,17 @@ class SpectralBaseAlgorithm
 
         /**
          * \brief Virtual function for current correction in Fourier space
-         * (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010).
-         * This virtual function is not pure: it can be overridden in derived
-         * classes (e.g. PsatdAlgorithm, PMLPsatdAlgorithm), but a base-class
-         * implementation is provided (empty implementation in this case).
+         * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>).
+         * This virtual function is pure and must be defined in derived classes.
          *
-         * \param[in,out] field_data all fields in Fourier space
-         * \param[in,out] current    three-dimensional array of unique pointers to MultiFab
-         *                           storing the three components of the current density
-         * \param[in]     rho        unique pointer to MultiFab storing the charge density
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
          */
-        virtual void CurrentCorrection ( SpectralFieldData& field_data,
-                                         std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
-                                         const std::unique_ptr<amrex::MultiFab>& rho ) {};
+        virtual void CurrentCorrection (SpectralFieldData& field_data,
+                                        std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+                                        const std::unique_ptr<amrex::MultiFab>& rho) = 0;
 
         /**
          * \brief Virtual function for Vay current deposition in Fourier space
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H
index 833a61aec..b0a30de50 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H
@@ -30,19 +30,17 @@ class SpectralBaseAlgorithmRZ
 
         /**
          * \brief Virtual function for current correction in Fourier space
-         * (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010).
-         * This virtual function is not pure: it can be overridden in derived
-         * classes (e.g. PsatdAlgorithmRZ), but a base-class
-         * implementation is provided (empty implementation in this case).
+         * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>).
+         * This virtual function is pure and must be defined in derived classes.
          *
-         * \param[in,out] field_data all fields in Fourier space
-         * \param[in,out] current    two-dimensional array of unique pointers to MultiFab
-         *                           storing the three components of the current density
-         * \param[in]     rho        unique pointer to MultiFab storing the charge density
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
          */
         virtual void CurrentCorrection ( SpectralFieldDataRZ& field_data,
                                          std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
-                                         const std::unique_ptr<amrex::MultiFab>& rho ) {};
+                                         const std::unique_ptr<amrex::MultiFab>& rho ) = 0;
 
         /**
          * \brief Compute spectral divergence of E