Do Not Fill Guard Cells with Inverse FFTs, Unless for Field Damping (#2045)

* Do Not Always Fill Guard Cells with Inverse FFTs

* Query psatd.fill_guards from Inputs

* Clean Up and Reduce Style Changes

* Fix Bug for Periodic Single Box

* Clean Up and Reduce Style Changes

* Fix Bug for RZ PSATD

* Remove Input Parameter, Default 0 Unless Damping

* Fix CI Tests (2D)

* Fix CI Tests (3D)

8 files changed, 35 insertions, 16 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H
index ecbc1578b..577ded61f 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H
@@ -32,6 +32,7 @@ class ComovingPsatdAlgorithm : public SpectralBaseAlgorithm
                                 const int norder_y,
                                 const int norder_z,
                                 const bool nodal,
+                                const amrex::IntVect& fill_guards,
                                 const amrex::Array<amrex::Real,3>& v_comoving,
                                 const amrex::Real dt,
                                 const bool update_with_rho);
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp
index 64920d325..d78ece8f5 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp
@@ -24,11 +24,12 @@ ComovingPsatdAlgorithm::ComovingPsatdAlgorithm (const SpectralKSpace& spectral_k
                                                 const DistributionMapping& dm,
                                                 const int norder_x, const int norder_y,
                                                 const int norder_z, const bool nodal,
+                                                const amrex::IntVect& fill_guards,
                                                 const amrex::Array<amrex::Real, 3>& v_comoving,
                                                 const amrex::Real dt,
                                                 const bool update_with_rho)
      // Members initialization
-     : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal),
+     : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal, fill_guards),
        // Initialize the infinite-order k vectors (the argument n_order = -1 selects
        // the infinite order option, the argument nodal = false is then irrelevant)
        kx_vec(spectral_kspace.getModifiedKComponent(dm, 0, -1, false)),
@@ -424,6 +425,8 @@ ComovingPsatdAlgorithm::CurrentCorrection (const int lev,
     field_data.ForwardTransform(lev, *rho, Idx::rho_old, 0);
     field_data.ForwardTransform(lev, *rho, Idx::rho_new, 1);
 
+    const amrex::IntVect& fill_guards = m_fill_guards;
+
     // Loop over boxes
     for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
 
@@ -504,9 +507,9 @@ ComovingPsatdAlgorithm::CurrentCorrection (const int lev,
     }
 
     // Backward Fourier transform of J
-    field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0);
-    field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0);
-    field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0);
+    field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0, fill_guards);
+    field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0, fill_guards);
+    field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0, fill_guards);
 }
 
 void
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
index 7e17eda07..83c93b25d 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
@@ -31,6 +31,7 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm
                          const amrex::DistributionMapping& dm,
                          const int norder_x, const int norder_y,
                          const int norder_z, const bool nodal,
+                         const amrex::IntVect& fill_guards,
                          const amrex::Real dt,
                          const bool dive_cleaning,
                          const bool divb_cleaning);
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
index d28b5218a..e7459d0f9 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
@@ -32,10 +32,11 @@ using namespace amrex;
 PMLPsatdAlgorithm::PMLPsatdAlgorithm(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 int norder_z, const bool nodal,
+                                     const amrex::IntVect& fill_guards, const Real dt,
                                      const bool dive_cleaning, const bool divb_cleaning)
      // Initialize members of base class
-     : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal),
+     : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal, fill_guards),
        m_dt(dt),
        m_dive_cleaning(dive_cleaning),
        m_divb_cleaning(divb_cleaning)
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
index a6b2e3f7e..a53319327 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
@@ -51,6 +51,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
             const int norder_y,
             const int norder_z,
             const bool nodal,
+            const amrex::IntVect& fill_guards,
             const amrex::Array<amrex::Real,3>& v_galilean,
             const amrex::Real dt,
             const bool update_with_rho,
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index b454d79ba..bc74a36a2 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -33,13 +33,14 @@ PsatdAlgorithm::PsatdAlgorithm(
     const int norder_y,
     const int norder_z,
     const bool nodal,
+    const amrex::IntVect& fill_guards,
     const amrex::Array<amrex::Real,3>& v_galilean,
     const amrex::Real dt,
     const bool update_with_rho,
     const bool time_averaging,
     const bool J_linear_in_time)
     // Initializer list
-    : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal),
+    : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal, fill_guards),
     // 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
@@ -859,6 +860,8 @@ PsatdAlgorithm::CurrentCorrection (
     field_data.ForwardTransform(lev, *rho, Idx::rho_old, 0);
     field_data.ForwardTransform(lev, *rho, Idx::rho_new, 1);
 
+    const amrex::IntVect& fill_guards = m_fill_guards;
+
     // Loop over boxes
     for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
 
@@ -951,9 +954,9 @@ PsatdAlgorithm::CurrentCorrection (
     }
 
     // Backward Fourier transform of J
-    field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0);
-    field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0);
-    field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0);
+    field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0, fill_guards);
+    field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0, fill_guards);
+    field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0, fill_guards);
 }
 
 void
@@ -974,6 +977,8 @@ PsatdAlgorithm::VayDeposition (
     field_data.ForwardTransform(lev, *current[1], Idx::Jy, 0, IntVect(1));
     field_data.ForwardTransform(lev, *current[2], Idx::Jz, 0, IntVect(1));
 
+    const amrex::IntVect& fill_guards = m_fill_guards;
+
     // Loop over boxes
     for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi)
     {
@@ -1028,9 +1033,9 @@ PsatdAlgorithm::VayDeposition (
     }
 
     // Backward Fourier transform of J
-    field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0);
-    field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0);
-    field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0);
+    field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0, fill_guards);
+    field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0, fill_guards);
+    field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0, fill_guards);
 }
 
 #endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
index 624d7870c..f412231b7 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
@@ -80,6 +80,8 @@ class SpectralBaseAlgorithm
 
     protected: // Meant to be used in the subclasses
 
+        amrex::IntVect m_fill_guards;
+
         using SpectralRealCoefficients = \
             amrex::FabArray< amrex::BaseFab <amrex::Real> >;
         using SpectralComplexCoefficients = \
@@ -91,7 +93,8 @@ class SpectralBaseAlgorithm
         SpectralBaseAlgorithm(const SpectralKSpace& spectral_kspace,
                               const amrex::DistributionMapping& dm,
                               const int norder_x, const int norder_y,
-                              const int norder_z, const bool nodal);
+                              const int norder_z, const bool nodal,
+                              const amrex::IntVect& fill_guards);
 
         // Modified finite-order vectors
         KVectorComponent modified_kx_vec;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp
index e57302cc4..4445705cb 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp
@@ -30,8 +30,10 @@ using namespace amrex;
 SpectralBaseAlgorithm::SpectralBaseAlgorithm(const SpectralKSpace& spectral_kspace,
     const amrex::DistributionMapping& dm,
     const int norder_x, const int norder_y,
-    const int norder_z, const bool nodal):
+    const int norder_z, const bool nodal,
+    const amrex::IntVect& fill_guards):
     // Compute and assign the modified k vectors
+        m_fill_guards(fill_guards),
         modified_kx_vec(spectral_kspace.getModifiedKComponent(dm,0,norder_x,nodal)),
 #if (AMREX_SPACEDIM==3)
         modified_ky_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_y,nodal)),
@@ -62,6 +64,8 @@ SpectralBaseAlgorithm::ComputeSpectralDivE (
     field_data.ForwardTransform(lev, *Efield[1], Idx::Ey, 0 );
     field_data.ForwardTransform(lev, *Efield[2], Idx::Ez, 0 );
 
+    const amrex::IntVect& fill_guards = m_fill_guards;
+
     // Loop over boxes
     for (MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
 
@@ -101,5 +105,5 @@ SpectralBaseAlgorithm::ComputeSpectralDivE (
     }
 
     // Backward Fourier transform
-    field_data.BackwardTransform(lev, divE, Idx::divE, 0 );
+    field_data.BackwardTransform(lev, divE, Idx::divE, 0, fill_guards);
 }