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)

12 files changed, 69 insertions, 24 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);
 }
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
index 4999a268d..e7764627b 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
@@ -85,7 +85,8 @@ class SpectralFieldData
             ForwardTransform(lev, mf, field_index, i_comp, mf.ixType().toIntVect());
         }
 
-        void BackwardTransform (const int lev, amrex::MultiFab& mf, const int field_index, const int i_comp);
+        void BackwardTransform (const int lev, amrex::MultiFab& mf, const int field_index,
+                                const int i_comp, const amrex::IntVect& fill_guards);
 
         // `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 5d8c52cc0..b26ac4943 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
@@ -223,10 +223,11 @@ SpectralFieldData::ForwardTransform (const int lev,
 /* \brief Transform spectral field specified by `field_index` back to
  * real space, and store it in the component `i_comp` of `mf` */
 void
-SpectralFieldData::BackwardTransform( const int lev,
+SpectralFieldData::BackwardTransform (const int lev,
                                       MultiFab& mf,
                                       const int field_index,
-                                      const int i_comp )
+                                      const int i_comp,
+                                      const amrex::IntVect& fill_guards)
 {
     amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
 
@@ -248,6 +249,9 @@ SpectralFieldData::BackwardTransform( const int lev,
     const int sk = (is_nodal_z) ? 1 : 0;
 #endif
 
+    // Numbers of guard cells
+    const amrex::IntVect& mf_ng = mf.nGrowVect();
+
     // Loop over boxes
     // Note: we do NOT OpenMP parallelize here, since we use OpenMP threads for
     //       the iFFTs on each box!
@@ -295,7 +299,7 @@ SpectralFieldData::BackwardTransform( const int lev,
         // Copy the temporary field tmpRealField to the real-space field mf and
         // normalize, dividing by N, since (FFT + inverse FFT) results in a factor N
         {
-            amrex::Box const& mf_box = (m_periodic_single_box) ? mfi.validbox() : mfi.fabbox();
+            amrex::Box mf_box = (m_periodic_single_box) ? mfi.validbox() : mfi.fabbox();
             amrex::Array4<amrex::Real> mf_arr = mf[mfi].array();
             amrex::Array4<const amrex::Real> tmp_arr = tmpRealField[mfi].array();
 
@@ -317,6 +321,16 @@ SpectralFieldData::BackwardTransform( const int lev,
 #elif (AMREX_SPACEDIM == 3)
             const int lo_k = amrex::lbound(mf_box).z;
 #endif
+            // If necessary, do not fill the guard cells
+            // (shrink box by passing negative number of cells)
+            if (m_periodic_single_box == false)
+            {
+                for (int dir = 0; dir < AMREX_SPACEDIM; dir++)
+                {
+                    if (static_cast<bool>(fill_guards[dir]) == false) mf_box.grow(dir, -mf_ng[dir]);
+                }
+            }
+
             // Loop over cells within full box, including ghost cells
             ParallelFor(mf_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
             {
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index e5d84a2af..45ba1a193 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -71,6 +71,7 @@ class SpectralSolver
                         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::Array<amrex::Real,3>& v_galilean,
                         const amrex::Array<amrex::Real,3>& v_comoving,
                         const amrex::RealVect dx,
@@ -184,6 +185,10 @@ class SpectralSolver
             field_data.fields.mult(scale_factor, icomp, 1);
         }
 
+    protected:
+
+        amrex::IntVect m_fill_guards;
+
     private:
 
         void ReadParameters ();
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index 89a7ce1f5..113ea97c3 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -23,6 +23,7 @@ SpectralSolver::SpectralSolver(
                 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::Array<amrex::Real,3>& v_galilean,
                 const amrex::Array<amrex::Real,3>& v_comoving,
                 const amrex::RealVect dx, const amrex::Real dt,
@@ -45,24 +46,29 @@ SpectralSolver::SpectralSolver(
 
     if (pml) {
         algorithm = std::make_unique<PMLPsatdAlgorithm>(
-            k_space, dm, norder_x, norder_y, norder_z, nodal, dt, dive_cleaning, divb_cleaning);
+            k_space, dm, norder_x, norder_y, norder_z, nodal,
+            fill_guards, dt, dive_cleaning, divb_cleaning);
     }
     else {
         // 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);
+                k_space, dm, norder_x, norder_y, norder_z, nodal,
+                fill_guards, v_comoving, dt, update_with_rho);
         }
         // PSATD algorithms: standard, Galilean, or averaged Galilean
         else {
             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, J_linear_in_time);
+                k_space, dm, norder_x, norder_y, norder_z, nodal, fill_guards,
+                v_galilean, dt, update_with_rho, fft_do_time_averaging, J_linear_in_time);
         }
     }
 
     // - Initialize arrays for fields in spectral space + FFT plans
     field_data = SpectralFieldData( lev, realspace_ba, k_space, dm,
                     algorithm->getRequiredNumberOfFields(), periodic_single_box);
+
+    m_fill_guards = fill_guards;
 }
 
 void
@@ -82,7 +88,7 @@ SpectralSolver::BackwardTransform( const int lev,
                                    const int i_comp )
 {
     WARPX_PROFILE("SpectralSolver::BackwardTransform");
-    field_data.BackwardTransform( lev, mf, field_index, i_comp );
+    field_data.BackwardTransform(lev, mf, field_index, i_comp, m_fill_guards);
 }
 
 void