Vay current deposition (#1051)

* Added stub for current correction in RZ spectral solver

* Start implementation of Vay deposition

* Continue implementation of Vay deposition

* Correct deposition of D

* Add phase shift for staggered currents

* Small clean-up

* Fix units in deposition of D

* Implement average of cumulative sum (needs bug fix)

* Start fixing bug in average of cumulative sum

* Still debugging

* Cumulative sums should be correct now

* Subtract averages of cumulative sums:

- current implementation: cumulative sums, inverse Fourier transform, subtraction of averages
- needs to be tested (including units of D after Vay deposition)
- needs to be shortened (too many loops over boxes and ParallelFors)

* [skip CI] Clean up and fix units

* Still fixing units

* [skip CI] Remove temporarily averages of cumulative sums

* [skip CI] Remove distinction between staggered and nodal

* Vay and Esirkepov similar results on periodic single box:

TODO:
- debug (charge not conserved);
- try using compute_shifted_shape_factor as in Esirkepov deposition;
- clean up;
- try on multiple boxes and with correction of mode at 0 frequency.

* [skip CI] Clean up

* Fix bug in 3D deposition

* [skip CI] Clean up

* Fix 2D and 3D implementation:

- simulation results agree between direct and Vay deposition in both 2D and 3D
- Travis CI tests should pass except for check of charge conservation (debug)

* Small clean-up

* Fix bug when compiling in RZ geometry

* Add benchmark json files (will be reset later)

* Do not set zero current at zero frequency

* [skip CI] Revert last commit and clean up

* Fix small bug after reverting commit

* Set nodal test first on Travis

* Fix benchmark for nodal test in 3D

* Fix particle output for nodal test in 3D

* Fix bugs due to staggering

* Rename current nodal Travis tests

* Add Travis tests staggered in 2D and 3D

* Further clean-up after bug fix

* Abort when using Vay deposition with domain decomposition

* Add optional argument of index type to forward FFT

* Fourier shifts can be private members as before

* Small clean-up

* Clean up and improve Doxygen documentation

* Fix small bug in analysis script for 2D tests

* Fix tests (remove E and B fields from particle diags)

* Add option to fill guard cells and docs

* Fix value of last guard cell by enforcing periodicity

* Revert changes merged from #1121

* Clean up style

* Improve docs

* Fix forgotten alignment

* Improve docs

* Make base class functions VayDeposition pure

Co-authored-by: Dave Grote <dpgrote@lbl.gov>

16 files changed, 245 insertions, 12 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H
index f2ae61f33..45a35098b 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H
@@ -26,6 +26,20 @@ class AvgGalileanAlgorithm : public SpectralBaseAlgorithm
            const amrex::Array<amrex::Real, 3>& v_galilean,
            const amrex::Real dt);
 
+        /**
+         * \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
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
+         *
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         */
+        virtual void VayDeposition (SpectralFieldData& field_data,
+                                    std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
+
     private:
         SpectralRealCoefficients C_coef, S_ck_coef, C1_coef, C3_coef, S1_coef,S3_coef;
         SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef, Psi1_coef, Psi2_coef, Psi3_coef, Psi4_coef, A1_coef, A2_coef, Rhoold_coef, Rhonew_coef, Jcoef_coef;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp
index d9db745f1..0ba603031 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp
@@ -369,3 +369,10 @@ AvgGalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
                         });
     }
 };
+
+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/GalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
index 192c453ad..799eca8b1 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
@@ -26,6 +26,20 @@ class GalileanAlgorithm : public SpectralBaseAlgorithm
                                     const amrex::Array<amrex::Real, 3>& v_galilean,
                                     const amrex::Real dt);
 
+        /**
+         * \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
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
+         *
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         */
+        virtual void VayDeposition (SpectralFieldData& field_data,
+                                    std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
+
     private:
         SpectralRealCoefficients C_coef, S_ck_coef;
         SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
index 49e318e3c..e606b3232 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
@@ -243,4 +243,12 @@ void GalileanAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spe
         });
     }
 }
+
+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/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
index ef12ff12e..347012d5a 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
@@ -34,6 +34,20 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm
             return SpectralPMLIndex::n_fields;
         }
 
+        /**
+         * \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
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
+         *
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         */
+        virtual void VayDeposition (SpectralFieldData& field_data,
+                                    std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
+
     private:
         SpectralRealCoefficients C_coef, S_ck_coef;
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
index b3f2aa93d..d2f087706 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
@@ -154,4 +154,12 @@ void PMLPsatdAlgorithm::InitializeSpectralCoefficients (
         });
     }
 };
+
+void
+PMLPsatdAlgorithm::VayDeposition (SpectralFieldData& field_data,
+                                  std::array<std::unique_ptr<amrex::MultiFab>,3>& current)
+{
+    amrex::Abort("Vay deposition not implemented for PML PSATD");
+}
+
 #endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
index b492d9a18..bfa9283e6 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
@@ -53,6 +53,20 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
                                          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
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
+         *
+         * \param[in,out] field_data All fields in Fourier space
+         * \param[in,out] current    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         */
+        virtual void VayDeposition (SpectralFieldData& field_data,
+                                    std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
+
     private:
         SpectralRealCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef;
         amrex::Real m_dt;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index b8b64b80c..7f9fd3edb 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -256,7 +256,6 @@ PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data,
         // 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
@@ -297,4 +296,80 @@ PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data,
     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) {
+    // Profiling
+    WARPX_PROFILE("PsatdAlgorithm::VayDeposition");
+
+    using Idx = SpectralFieldIndex;
+
+    // Forward Fourier transform of D (temporarily stored in current):
+    // D is nodal and does not match the staggering of J, therefore we pass the
+    // actual staggering of D (IntVect(1)) to the ForwardTransform function
+    field_data.ForwardTransform(*current[0], Idx::Jx, 0, IntVect(1));
+    field_data.ForwardTransform(*current[1], Idx::Jy, 0, IntVect(1));
+    field_data.ForwardTransform(*current[2], Idx::Jz, 0, IntVect(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 modified 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();
+
+        // Loop over indices within one box
+        ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            using Idx = SpectralFieldIndex;
+
+            // Shortcuts for the values of D
+            const Complex Dx = fields(i,j,k,Idx::Jx);
+            const Complex Dy = fields(i,j,k,Idx::Jy);
+            const Complex Dz = fields(i,j,k,Idx::Jz);
+
+            // Imaginary unit
+            constexpr Complex I = Complex{0._rt, 1._rt};
+
+            // Modified k vector values
+            const amrex::Real kx_mod = modified_kx_arr[i];
+#if (AMREX_SPACEDIM==3)
+            const amrex::Real ky_mod = modified_ky_arr[j];
+            const amrex::Real kz_mod = modified_kz_arr[k];
+#else
+            constexpr amrex::Real ky_mod = 0._rt;
+            const     amrex::Real kz_mod = modified_kz_arr[j];
+#endif
+
+            // Compute Jx
+            if (kx_mod != 0._rt) fields(i,j,k,Idx::Jx) = I * Dx / kx_mod;
+            else                 fields(i,j,k,Idx::Jx) = 0._rt;
+
+#if (AMREX_SPACEDIM==3)
+            // Compute Jy
+            if (ky_mod != 0._rt) fields(i,j,k,Idx::Jy) = I * Dy / ky_mod;
+            else                 fields(i,j,k,Idx::Jy) = 0._rt;
+#endif
+
+            // 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;
+
+        });
+    }
+
+    // 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);
+}
 #endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
index 6a604155b..fc03f95f2 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
@@ -44,6 +44,20 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ
                                          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
+         * 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    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         */
+        virtual void VayDeposition (SpectralFieldDataRZ& field_data,
+                                    std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
+
     private:
         bool coefficients_initialized;
         // Note that dt is saved to use in InitializeSpectralCoefficients
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
index b094843ac..35a478f2c 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
@@ -204,3 +204,10 @@ PsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data,
 {
     amrex::Abort("Current correction not implemented in RZ");
 }
+
+void
+PsatdAlgorithmRZ::VayDeposition (SpectralFieldDataRZ& field_data,
+                                 std::array<std::unique_ptr<amrex::MultiFab>,3>& current)
+{
+    amrex::Abort("Vay deposition not implemented in RZ geometry");
+}
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
index 7143b7aa9..2a34d21ab 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
@@ -48,6 +48,18 @@ class SpectralBaseAlgorithm
                                          const std::unique_ptr<amrex::MultiFab>& rho ) {};
 
         /**
+         * \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 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    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) = 0;
+
+        /**
          * \brief Compute spectral divergence of E
          */
         void ComputeSpectralDivE ( SpectralFieldData& field_data,
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H
index ea50e979a..833a61aec 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H
@@ -51,6 +51,18 @@ class SpectralBaseAlgorithmRZ
                                    const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
                                    amrex::MultiFab& divE );
 
+        /**
+         * \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 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    Array of unique pointers to \c MultiFab storing
+         *                           the three components of the current density
+         */
+        virtual void VayDeposition (SpectralFieldDataRZ& field_data,
+                                    std::array<std::unique_ptr<amrex::MultiFab>,3>& current) = 0;
+
     protected: // Meant to be used in the subclasses
 
         using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
index b445054cc..f48272744 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
@@ -55,10 +55,17 @@ class SpectralFieldData
         SpectralFieldData() = default; // Default constructor
         SpectralFieldData& operator=(SpectralFieldData&& field_data) = default;
         ~SpectralFieldData();
-        void ForwardTransform( const amrex::MultiFab& mf,
-                               const int field_index, const int i_comp);
-        void BackwardTransform( amrex::MultiFab& mf,
-                               const int field_index, const int i_comp);
+
+        void ForwardTransform (const amrex::MultiFab& mf, const int field_index,
+                               const int i_comp, const amrex::IntVect& stag);
+        AMREX_FORCE_INLINE
+        void ForwardTransform (const amrex::MultiFab& mf, const int field_index, const int i_comp)
+        {
+            ForwardTransform(mf, field_index, i_comp, mf.ixType().toIntVect());
+        }
+
+        void BackwardTransform (amrex::MultiFab& mf, const int field_index, const int i_comp);
+
         // `fields` stores fields in spectral space, as multicomponent FabArray
         SpectralField fields;
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
index 0f49e695b..e99e7d57b 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
@@ -94,17 +94,16 @@ SpectralFieldData::~SpectralFieldData()
  *  to spectral space, and store the corresponding result internally
  *  (in the spectral field specified by `field_index`) */
 void
-SpectralFieldData::ForwardTransform( const MultiFab& mf,
-                                     const int field_index,
-                                     const int i_comp )
+SpectralFieldData::ForwardTransform (const MultiFab& mf, const int field_index,
+                                     const int i_comp, const IntVect& stag)
 {
     // Check field index type, in order to apply proper shift in spectral space
-    const bool is_nodal_x = mf.is_nodal(0);
+    const bool is_nodal_x = (stag[0] == amrex::IndexType::NODE) ? true : false;
 #if (AMREX_SPACEDIM == 3)
-    const bool is_nodal_y = mf.is_nodal(1);
-    const bool is_nodal_z = mf.is_nodal(2);
+    const bool is_nodal_y = (stag[1] == amrex::IndexType::NODE) ? true : false;
+    const bool is_nodal_z = (stag[2] == amrex::IndexType::NODE) ? true : false;
 #else
-    const bool is_nodal_z = mf.is_nodal(1);
+    const bool is_nodal_z = (stag[1] == amrex::IndexType::NODE) ? true : false;
 #endif
 
     // Loop over boxes
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index 91cfc84f1..ba805c790 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -84,6 +84,20 @@ class SpectralSolver
              algorithm->CurrentCorrection( field_data, current, rho );
         };
 
+        /**
+         * \brief Public interface to call the virtual function \c VayDeposition,
+         * declared in the base class SpectralBaseAlgorithm and defined in its
+         * derived classes, from objects of class SpectralSolver through the private
+         * unique pointer \c algorithm.
+         *
+         * \param[in,out] current Array of unique pointers to \c MultiFab storing
+         *                        the three components of the current density
+         */
+        void VayDeposition (std::array<std::unique_ptr<amrex::MultiFab>,3>& current)
+        {
+            algorithm->VayDeposition(field_data, current);
+        }
+
     private:
         void ReadParameters ();
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
index 5762e1436..e1179d9cf 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
@@ -104,6 +104,20 @@ class SpectralSolverRZ
              algorithm->CurrentCorrection( field_data, current, rho );
         };
 
+        /**
+         * \brief Public interface to call the virtual function \c VayDeposition,
+         * declared in the base class SpectralBaseAlgorithmRZ and defined in its
+         * derived classes, from objects of class SpectralSolverRZ through the private
+         * unique pointer \c algorithm.
+         *
+         * \param[in,out] current Array of unique pointers to \c MultiFab storing
+         *                        the three components of the current density
+         */
+        void VayDeposition (std::array<std::unique_ptr<amrex::MultiFab>,3>& current)
+        {
+            algorithm->VayDeposition(field_data, current);
+        }
+
     private:
 
         SpectralFieldDataRZ field_data; // Store field in spectral space