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)

20 files changed, 172 insertions, 109 deletions

diff --git a/Regression/Checksum/benchmarks_json/comoving_2d_psatd_hybrid.json b/Regression/Checksum/benchmarks_json/comoving_2d_psatd_hybrid.json
index fa936c3c7..ec447bc56 100644
--- a/Regression/Checksum/benchmarks_json/comoving_2d_psatd_hybrid.json
+++ b/Regression/Checksum/benchmarks_json/comoving_2d_psatd_hybrid.json
@@ -2,43 +2,43 @@
   "beam": {
     "particle_cpu": 0.0,
     "particle_id": 500500.0,
-    "particle_momentum_x": 7.587450161877821e-19,
-    "particle_momentum_y": 4.515786832153623e-19,
-    "particle_momentum_z": 9.228375405679396e-18,
-    "particle_position_x": 0.005388318043043357,
-    "particle_position_y": 0.16843460833541193,
+    "particle_momentum_x": 7.587449921798846e-19,
+    "particle_momentum_y": 4.515786832176223e-19,
+    "particle_momentum_z": 9.228375409541606e-18,
+    "particle_position_x": 0.005388317800845895,
+    "particle_position_y": 0.16843460846512465,
     "particle_weight": 3120754537230.3823
   },
   "electrons": {
     "particle_cpu": 155246.0,
     "particle_id": 26701298988.0,
-    "particle_momentum_x": 9.19759655761835e-19,
-    "particle_momentum_y": 8.111395741629531e-19,
-    "particle_momentum_z": 5.717012342804088e-16,
-    "particle_position_x": 6.634386766502452,
-    "particle_position_y": 58.08979821925373,
+    "particle_momentum_x": 9.197595900579815e-19,
+    "particle_momentum_y": 8.111395743160647e-19,
+    "particle_momentum_z": 5.717012342689313e-16,
+    "particle_position_x": 6.634386766537618,
+    "particle_position_y": 58.0897982192507,
     "particle_weight": 3.212795210423338e+18
   },
   "ions": {
     "particle_cpu": 155246.0,
     "particle_id": 26830093836.0,
-    "particle_momentum_x": 2.0052369737454595e-18,
-    "particle_momentum_y": 8.210971641417303e-19,
-    "particle_momentum_z": 1.0089821790887407e-12,
-    "particle_position_x": 6.622297130412207,
+    "particle_momentum_x": 2.005236974560188e-18,
+    "particle_momentum_y": 8.21097164218749e-19,
+    "particle_momentum_z": 1.0089821790887435e-12,
+    "particle_position_x": 6.622297130412235,
     "particle_position_y": 58.08877602017003,
     "particle_weight": 3.212794492131793e+18
   },
   "lev=0": {
-    "Bx": 1305297.4828241316,
-    "By": 4009300.510439389,
-    "Bz": 203715.33013685,
-    "Ex": 1203889010242205.5,
-    "Ey": 408608837447717.06,
-    "Ez": 217844734969419.84,
-    "jx": 3.5197389872605124e+16,
-    "jy": 2.314621849516091e+16,
-    "jz": 4.137802990284829e+17,
-    "rho": 1352926282.578562
+    "Bx": 1305297.485181685,
+    "By": 4009300.5547388536,
+    "Bz": 203715.3311601148,
+    "Ex": 1203889033875429.0,
+    "Ey": 408608837197114.75,
+    "Ez": 217844729568617.62,
+    "jx": 3.519739071285825e+16,
+    "jy": 2.3146218523962868e+16,
+    "jz": 4.137803403281171e+17,
+    "rho": 1352926431.87181
   }
 }
\ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/galilean_2d_psatd_hybrid.json b/Regression/Checksum/benchmarks_json/galilean_2d_psatd_hybrid.json
index 36159b1cb..b0021fc06 100644
--- a/Regression/Checksum/benchmarks_json/galilean_2d_psatd_hybrid.json
+++ b/Regression/Checksum/benchmarks_json/galilean_2d_psatd_hybrid.json
@@ -2,43 +2,43 @@
   "beam": {
     "particle_cpu": 0.0,
     "particle_id": 500500.0,
-    "particle_momentum_x": 7.781818973078852e-19,
-    "particle_momentum_y": 4.515677165595182e-19,
-    "particle_momentum_z": 8.375527032802002e-18,
-    "particle_position_x": 0.006743036325434187,
-    "particle_position_y": 0.1716091883553772,
+    "particle_momentum_x": 7.781818973095813e-19,
+    "particle_momentum_y": 4.515677165586177e-19,
+    "particle_momentum_z": 8.375527032801872e-18,
+    "particle_position_x": 0.006743036325498549,
+    "particle_position_y": 0.171609188355419,
     "particle_weight": 3120754537230.3823
   },
   "electrons": {
     "particle_cpu": 155383.0,
     "particle_id": 26721507984.0,
-    "particle_momentum_x": 7.343166622433748e-19,
-    "particle_momentum_y": 5.669478712783679e-19,
-    "particle_momentum_z": 5.697759351248577e-16,
-    "particle_position_x": 6.639863939125528,
-    "particle_position_y": 58.2105156240904,
+    "particle_momentum_x": 7.343166618394246e-19,
+    "particle_momentum_y": 5.669478713153082e-19,
+    "particle_momentum_z": 5.697759351247803e-16,
+    "particle_position_x": 6.639863939125319,
+    "particle_position_y": 58.21051562409041,
     "particle_weight": 3.215903744897337e+18
   },
   "ions": {
     "particle_cpu": 155428.0,
     "particle_id": 26814790002.0,
-    "particle_momentum_x": 1.6709338576817887e-18,
-    "particle_momentum_y": 5.73611349305925e-19,
+    "particle_momentum_x": 1.6709338568912299e-18,
+    "particle_momentum_y": 5.736113493536961e-19,
     "particle_momentum_z": 1.010167756546284e-12,
     "particle_position_x": 6.63004997278591,
-    "particle_position_y": 58.249147619440286,
+    "particle_position_y": 58.2491476194403,
     "particle_weight": 3.21662600232034e+18
   },
   "lev=0": {
-    "Bx": 1225047.3704144938,
-    "By": 3461610.5122063416,
-    "Bz": 194622.52078763963,
-    "Ex": 1039200462527705.4,
-    "Ey": 384410222541170.8,
-    "Ez": 200388928565719.44,
-    "jx": 3.100345182047955e+16,
-    "jy": 2.2227570738047532e+16,
-    "jz": 2.9254388838934234e+17,
-    "rho": 995460563.3392936
+    "Bx": 1225047.370205577,
+    "By": 3461610.509386122,
+    "Bz": 194622.52070961153,
+    "Ex": 1039200462646995.5,
+    "Ey": 384410222605049.8,
+    "Ez": 200388928770113.06,
+    "jx": 3.1003451811671292e+16,
+    "jy": 2.222757073772597e+16,
+    "jz": 2.9254388838862925e+17,
+    "rho": 995460563.3297532
   }
 }
\ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/multi_J_2d_psatd.json b/Regression/Checksum/benchmarks_json/multi_J_2d_psatd.json
index 7788a9d85..3398f7ba6 100644
--- a/Regression/Checksum/benchmarks_json/multi_J_2d_psatd.json
+++ b/Regression/Checksum/benchmarks_json/multi_J_2d_psatd.json
@@ -2,47 +2,47 @@
   "driver": {
     "particle_cpu": 0.0,
     "particle_id": 5000050000.0,
-    "particle_momentum_x": 4.2468004298680323e-16,
+    "particle_momentum_x": 4.246800430279664e-16,
     "particle_momentum_y": 0.0,
-    "particle_momentum_z": 9.82279023786301e-09,
-    "particle_position_x": 0.3979349343731913,
+    "particle_momentum_z": 9.822790237863e-09,
+    "particle_position_x": 0.3979349343731916,
     "particle_position_y": 55.175199387968405,
     "particle_weight": 124830181489215.27
   },
   "lev=0": {
     "Bx": 0.0,
-    "By": 938797.6239373332,
+    "By": 938798.4013860568,
     "Bz": 0.0,
-    "Ex": 275592022698142.3,
+    "Ex": 275591837601629.0,
     "Ey": 0.0,
-    "Ez": 93204573102613.78,
-    "jx": 9857319741191368.0,
+    "Ez": 93197772358526.53,
+    "jx": 9857318930321166.0,
     "jy": 0.0,
-    "jz": 6.4634613927785736e+16,
-    "rho": 221319690.13710034,
+    "jz": 6.463459953111897e+16,
+    "rho": 221319857.83191392,
     "rho_driver": 2562225.1199331162,
     "rho_driver_back": 0.0,
-    "rho_plasma_e": 2698164888.5847692,
-    "rho_plasma_p": 2699287913.2703176
+    "rho_plasma_e": 2698164828.605695,
+    "rho_plasma_p": 2699287913.2711463
   },
   "plasma_e": {
     "particle_cpu": 58858.0,
     "particle_id": 7388900221.0,
-    "particle_momentum_x": 7.296651628257412e-19,
+    "particle_momentum_x": 7.296651098332414e-19,
     "particle_momentum_y": 0.0,
-    "particle_momentum_z": 4.494602499033915e-17,
-    "particle_position_x": 5.338124877733829,
-    "particle_position_y": 26.54016202205885,
+    "particle_momentum_z": 4.4946035352884545e-17,
+    "particle_position_x": 5.338124868002372,
+    "particle_position_y": 26.540162019472735,
     "particle_weight": 1.3186130303279046e+17
   },
   "plasma_p": {
     "particle_cpu": 58880.0,
     "particle_id": 7417282880.0,
-    "particle_momentum_x": 1.4768074009968963e-18,
+    "particle_momentum_x": 1.476807478359322e-18,
     "particle_momentum_y": 0.0,
-    "particle_momentum_z": 7.942647962730444e-14,
-    "particle_position_x": 5.290048066122284,
-    "particle_position_y": 26.56734131597245,
+    "particle_momentum_z": 7.942647961709633e-14,
+    "particle_position_x": 5.290048066126035,
+    "particle_position_y": 26.567341315980233,
     "particle_weight": 1.3191059027779915e+17
   }
 }
\ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json b/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json
index f3c0b0a48..ad0acf8f7 100644
--- a/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json
+++ b/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json
@@ -1,11 +1,11 @@
 {
   "lev=0": {
-    "Bx": 1.55895366473558e-08,
-    "By": 1.5015206116838625e-08,
-    "Bz": 1.516725065646249e-08,
-    "Ex": 19.175719716905533,
-    "Ey": 19.25088294791974,
-    "Ez": 19.244587842285785,
+    "Bx": 1.3610687381384125e-08,
+    "By": 1.361341814453787e-08,
+    "Bz": 1.361623444282223e-08,
+    "Ex": 18.733019975918303,
+    "Ey": 18.73132517981054,
+    "Ez": 18.731454226529113,
     "rho": 0.00017713857500552797
   }
 }
\ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/pml_x_psatd.json b/Regression/Checksum/benchmarks_json/pml_x_psatd.json
index f845c7199..042e9f2cf 100644
--- a/Regression/Checksum/benchmarks_json/pml_x_psatd.json
+++ b/Regression/Checksum/benchmarks_json/pml_x_psatd.json
@@ -1,12 +1,12 @@
 {
   "lev=0": {
-    "Bx": 1.2110998243414779e-08,
-    "By": 1.5545485705380978e-08,
-    "Bz": 9.134211013787592e-09,
-    "Ex": 4.050591183040409,
-    "Ey": 4.248614715015242,
-    "Ez": 3.3401282934200776,
-    "divE": 188920.76108729007,
+    "Bx": 1.159349225951124e-08,
+    "By": 1.5110587450944313e-08,
+    "Bz": 8.985531772997968e-09,
+    "Ex": 3.93868514084609,
+    "Ey": 4.13492225197527,
+    "Ez": 3.3085920402662996,
+    "divE": 188916.91920600203,
     "rho": 1.6719346125882405e-06
   }
 }
\ No newline at end of file
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp
index 234cfbabb..033e497a6 100644
--- a/Source/BoundaryConditions/PML.cpp
+++ b/Source/BoundaryConditions/PML.cpp
@@ -603,9 +603,9 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& /*g
         Array<Real,3> const v_comoving_zero = {0., 0., 0.};
         realspace_ba.enclosedCells().grow(nge); // cell-centered + guard cells
         spectral_solver_fp = std::make_unique<SpectralSolver>(lev, realspace_ba, dm,
-            nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, dx, dt, in_pml,
-            periodic_single_box, update_with_rho, fft_do_time_averaging,
-            J_linear_in_time, m_dive_cleaning, m_divb_cleaning);
+            nox_fft, noy_fft, noz_fft, do_nodal, WarpX::fill_guards, v_galilean_zero,
+            v_comoving_zero, dx, dt, in_pml, periodic_single_box, update_with_rho,
+            fft_do_time_averaging, J_linear_in_time, m_dive_cleaning, m_divb_cleaning);
 #endif
     }
 
@@ -711,9 +711,9 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& /*g
             Array<Real,3> const v_comoving_zero = {0., 0., 0.};
             realspace_cba.enclosedCells().grow(nge); // cell-centered + guard cells
             spectral_solver_cp = std::make_unique<SpectralSolver>(lev, realspace_cba, cdm,
-                nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, cdx, dt, in_pml,
-                periodic_single_box, update_with_rho, fft_do_time_averaging,
-                J_linear_in_time, m_dive_cleaning, m_divb_cleaning);
+                nox_fft, noy_fft, noz_fft, do_nodal, WarpX::fill_guards, v_galilean_zero,
+                v_comoving_zero, cdx, dt, in_pml, periodic_single_box, update_with_rho,
+                fft_do_time_averaging, J_linear_in_time, m_dive_cleaning, m_divb_cleaning);
 #endif
         }
     }
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
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 2cb5f4483..9b73174ea 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -152,6 +152,9 @@ public:
     // default is false for standard PSATD and true for Galilean PSATD (set in WarpX.cpp)
     bool update_with_rho = false;
 
+    // PSATD: Whether to fill the guard cells with inverse FFTs based on the boundary conditions
+    static amrex::IntVect fill_guards;
+
     // div(E) and div(B) cleaning
     static bool do_dive_cleaning;
     static bool do_divb_cleaning;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 0745b0941..3c5604a8b 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -99,6 +99,8 @@ Real WarpX::moving_window_v = std::numeric_limits<amrex::Real>::max();
 
 bool WarpX::fft_do_time_averaging = false;
 
+amrex::IntVect WarpX::fill_guards = amrex::IntVect(0);
+
 Real WarpX::quantum_xi_c2 = PhysConst::xi_c2;
 Real WarpX::gamma_boost = 1._rt;
 Real WarpX::beta_boost = 0._rt;
@@ -1089,6 +1091,18 @@ WarpX::ReadParameters ()
                 "field boundary in both lo and hi must be set to Damped for PSATD"
             );
         }
+
+        // Whether to fill the guard cells with inverse FFTs:
+        // WarpX::fill_guards = amrex::IntVect(0) by default,
+        // except for non-periodic directions with damping.
+        for (int dir = 0; dir < AMREX_SPACEDIM; dir++)
+        {
+            if (WarpX::field_boundary_lo[dir] == FieldBoundaryType::Damped ||
+                WarpX::field_boundary_hi[dir] == FieldBoundaryType::Damped)
+            {
+                WarpX::fill_guards[dir] = 1;
+            }
+        }
     }
 
     if (maxwell_solver_id != MaxwellSolverAlgo::PSATD ) {
@@ -1846,6 +1860,7 @@ void WarpX::AllocLevelSpectralSolver (amrex::Vector<std::unique_ptr<SpectralSolv
                                                 noy_fft,
                                                 noz_fft,
                                                 do_nodal,
+                                                WarpX::fill_guards,
                                                 m_v_galilean,
                                                 m_v_comoving,
                                                 dx_vect,