1 files changed, 78 insertions, 8 deletions

diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 82ca52c4c..a9267e4cd 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -55,8 +55,7 @@ using namespace amrex;
 #ifdef WARPX_USE_PSATD
 namespace {
 
-    void
-    ForwardTransformVect (
+    void ForwardTransformVect (
         const int lev,
 #ifdef WARPX_DIM_RZ
         SpectralSolverRZ& solver,
@@ -64,15 +63,35 @@ namespace {
         SpectralSolver& solver,
 #endif
         std::array<std::unique_ptr<amrex::MultiFab>,3>& vector_field,
-        const int compx, const int compy, const int compz)
+        const int compx, const int compy, const int compz,
+        const amrex::IntVect& stag_x, const amrex::IntVect& stag_y, const amrex::IntVect& stag_z)
     {
 #ifdef WARPX_DIM_RZ
+        amrex::ignore_unused(stag_x, stag_y, stag_z);
         solver.ForwardTransform(lev, *vector_field[0], compx, *vector_field[1], compy);
+        solver.ForwardTransform(lev, *vector_field[2], compz);
 #else
-        solver.ForwardTransform(lev, *vector_field[0], compx);
-        solver.ForwardTransform(lev, *vector_field[1], compy);
+        solver.ForwardTransform(lev, *vector_field[0], compx, stag_x);
+        solver.ForwardTransform(lev, *vector_field[1], compy, stag_y);
+        solver.ForwardTransform(lev, *vector_field[2], compz, stag_z);
 #endif
-        solver.ForwardTransform(lev, *vector_field[2], compz);
+    }
+
+    AMREX_FORCE_INLINE
+    void ForwardTransformVect (
+        const int lev,
+#ifdef WARPX_DIM_RZ
+        SpectralSolverRZ& solver,
+#else
+        SpectralSolver& solver,
+#endif
+        std::array<std::unique_ptr<amrex::MultiFab>,3>& vector_field,
+        const int compx, const int compy, const int compz)
+    {
+        ForwardTransformVect(lev, solver, vector_field, compx, compy, compz,
+                             (*vector_field[0]).ixType().toIntVect(),
+                             (*vector_field[1]).ixType().toIntVect(),
+                             (*vector_field[2]).ixType().toIntVect());
     }
 
     void
@@ -243,11 +262,41 @@ WarpX::PSATDForwardTransformJ ()
 
     for (int lev = 0; lev <= finest_level; ++lev)
     {
-        ForwardTransformVect(lev, *spectral_solver_fp[lev], current_fp[lev], idx_jx, idx_jy, idx_jz);
+        // With Vay's deposition, J stores a temporary current (D) that is later modified
+        // in Fourier space: its staggering matches that of rho and not J
+        amrex::IntVect jx_stag =
+            (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ?
+            WarpX::m_rho_nodal_flag : current_fp[lev][0]->ixType().toIntVect();
+
+        amrex::IntVect jy_stag =
+            (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ?
+            WarpX::m_rho_nodal_flag : current_fp[lev][1]->ixType().toIntVect();
+
+        amrex::IntVect jz_stag =
+            (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ?
+            WarpX::m_rho_nodal_flag : current_fp[lev][2]->ixType().toIntVect();
+
+        ForwardTransformVect(lev, *spectral_solver_fp[lev], current_fp[lev],
+                             idx_jx, idx_jy, idx_jz, jx_stag, jy_stag, jz_stag);
 
         if (spectral_solver_cp[lev])
         {
-            ForwardTransformVect(lev, *spectral_solver_cp[lev], current_cp[lev], idx_jx, idx_jy, idx_jz);
+            // With Vay's deposition, J stores a temporary current (D) that is later modified
+            // in Fourier space: its staggering matches that of rho and not J
+            jx_stag =
+                (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ?
+                WarpX::m_rho_nodal_flag : current_cp[lev][0]->ixType().toIntVect();
+
+            jy_stag =
+                (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ?
+                WarpX::m_rho_nodal_flag : current_cp[lev][1]->ixType().toIntVect();
+
+            jz_stag =
+                (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ?
+                WarpX::m_rho_nodal_flag : current_cp[lev][2]->ixType().toIntVect();
+
+            ForwardTransformVect(lev, *spectral_solver_cp[lev], current_cp[lev],
+                                 idx_jx, idx_jy, idx_jz, jx_stag, jy_stag, jz_stag);
         }
     }
 
@@ -348,6 +397,19 @@ void WarpX::PSATDCurrentCorrection ()
     }
 }
 
+void WarpX::PSATDVayDeposition ()
+{
+    for (int lev = 0; lev <= finest_level; ++lev)
+    {
+        spectral_solver_fp[lev]->VayDeposition();
+
+        if (spectral_solver_cp[lev])
+        {
+            spectral_solver_cp[lev]->VayDeposition();
+        }
+    }
+}
+
 void
 WarpX::PSATDPushSpectralFields ()
 {
@@ -476,6 +538,14 @@ WarpX::PushPSATD ()
         PSATDBackwardTransformJ();
     }
 
+    // Compute the current in Fourier space according to the Vay deposition scheme, and
+    // transform back to real space so that the Vay deposition is reflected in the diagnostics
+    if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay)
+    {
+        PSATDVayDeposition();
+        PSATDBackwardTransformJ();
+    }
+
 #ifdef WARPX_DIM_RZ
     if (pml_rz[0]) pml_rz[0]->PushPSATD(0);
 #endif