Implement WarpXSumGuardCells

3 files changed, 97 insertions, 33 deletions

diff --git a/Source/Parallelization/Make.package b/Source/Parallelization/Make.package
index cbb1b5234..3d1fcf1da 100644
--- a/Source/Parallelization/Make.package
+++ b/Source/Parallelization/Make.package
@@ -1,5 +1,6 @@
 CEXE_sources += WarpXComm.cpp
 CEXE_sources += WarpXRegrid.cpp
+CEXE_headers += WarpXSumGuardCells.H
 
 INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Parallelization
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Parallelization
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index 28971eb0c..348b31c8b 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -1,5 +1,6 @@
 #include <WarpX.H>
 #include <WarpX_f.H>
+#include <WarpXSumGuardCells.H>
 
 #include <AMReX_FillPatchUtil_F.H>
 
@@ -355,7 +356,8 @@ WarpX::SyncCurrent ()
 {
     BL_PROFILE("SyncCurrent()");
 
-    // Restrict fine patch current onto the coarse patch, before fine patch SumBoundary
+    // Restrict fine patch current onto the coarse patch, before
+    // summing the guard cells of the fine patch
     for (int lev = 1; lev <= finest_level; ++lev)
     {
         current_cp[lev][0]->setVal(0.0);
@@ -391,8 +393,9 @@ WarpX::SyncCurrent ()
                 bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]);
                 // Then swap j_fp and current_fp
                 std::swap(j_fp[lev][idim], current_fp[lev][idim]);
-                // At this point, current_fp may have false values close to the 
-                // edges of each FAB. This will be solved with a SumBoundary later.
+                // At this point, current_fp may have false values close to the
+                // edges of each FAB. This will be solved with when summing
+                // the guard cells later.
                 // j_fp contains the exact MultiFab current_fp before this step.
             }
         }
@@ -427,11 +430,11 @@ WarpX::SyncCurrent ()
     {
         const auto& period = Geom(lev).periodicity();
         // When using a bilinear filter with many passes, current_fp has
-        // temporarily more ghost cells here, so that its value inside 
+        // temporarily more ghost cells here, so that its value inside
         // the domain is correct at the end of this stage.
-        current_fp[lev][0]->SumBoundary(period);
-        current_fp[lev][1]->SumBoundary(period);
-        current_fp[lev][2]->SumBoundary(period);
+        WarpXSumGuardCells(*(current_fp[lev][0]),period);
+        WarpXSumGuardCells(*(current_fp[lev][1]),period);
+        WarpXSumGuardCells(*(current_fp[lev][2]),period);
     }
 
     // Add fine level's coarse patch to coarse level's fine patch
@@ -461,9 +464,9 @@ WarpX::SyncCurrent ()
     for (int lev = 1; lev <= finest_level; ++lev)
     {
         const auto& cperiod = Geom(lev-1).periodicity();
-        current_cp[lev][0]->SumBoundary(cperiod);
-        current_cp[lev][1]->SumBoundary(cperiod);
-        current_cp[lev][2]->SumBoundary(cperiod);
+        WarpXSumGuardCells(*(current_cp[lev][0]),cperiod);
+        WarpXSumGuardCells(*(current_cp[lev][1]),cperiod);
+        WarpXSumGuardCells(*(current_cp[lev][2]),cperiod);
     }
 
     if (WarpX::use_filter) {
@@ -476,7 +479,7 @@ WarpX::SyncCurrent ()
                 std::swap(j_fp[lev][idim], current_fp[lev][idim]);
                 // Then copy the interior of j_fp to current_fp.
                 MultiFab::Copy(*current_fp[lev][idim], *j_fp[lev][idim], 0, 0, 1, 0);
-                // current_fp has right number of ghost cells and 
+                // current_fp has right number of ghost cells and
                 // correct filtered values here.
             }
         }
@@ -556,7 +559,8 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof,
 {
     if (!rhof[0]) return;
 
-    // Restrict fine patch onto the coarse patch, before fine patch SumBoundary
+    // Restrict fine patch onto the coarse patch,
+    // before summing the guard cells of the fine patch
     for (int lev = 1; lev <= finest_level; ++lev)
     {
         rhoc[lev]->setVal(0.0);
@@ -607,7 +611,7 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof,
     for (int lev = 0; lev <= finest_level; ++lev)
     {
         const auto& period = Geom(lev).periodicity();
-        rhof[lev]->SumBoundary(period);
+        WarpXSumGuardCells( *(rhof[lev]), period, 0, rhof[lev]->nComp() );
     }
 
     // Add fine level's coarse patch to coarse level's fine patch
@@ -631,7 +635,7 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof,
     for (int lev = 1; lev <= finest_level; ++lev)
     {
         const auto& cperiod = Geom(lev-1).periodicity();
-        rhoc[lev]->SumBoundary(cperiod);
+        WarpXSumGuardCells( *(rhoc[lev]), cperiod, 0, rhoc[lev]->nComp() );
     }
 
     if (WarpX::use_filter) {
@@ -729,10 +733,9 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type)
             ng += bilinear_filter.stencil_length_each_dir-1;
             MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), 1, ng);
             bilinear_filter.ApplyStencil(jf, *j[idim]);
-            jf.SumBoundary(period);
-            MultiFab::Copy(*j[idim], jf, 0, 0, 1, 0);
+            WarpXSumGuardCells(*(j[idim]), jf, period);
         } else {
-            j[idim]->SumBoundary(period);
+            WarpXSumGuardCells(*(j[idim]), period);
         }
     }
 }
@@ -780,8 +783,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev)
             MultiFab::Add(jfb, jfc, 0, 0, 1, ng);
             mf.ParallelAdd(jfb, 0, 0, 1, ng, IntVect::TheZeroVector(), period);
 
-            jfc.SumBoundary(period);
-            MultiFab::Copy(*current_cp[lev+1][idim], jfc, 0, 0, 1, 0);
+            WarpXSumGuardCells(*current_cp[lev+1][idim], jfc, period);
         }
         else if (use_filter) // but no buffer
         {
@@ -792,8 +794,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev)
                         current_cp[lev+1][idim]->DistributionMap(), 1, ng);
             bilinear_filter.ApplyStencil(jf, *current_cp[lev+1][idim]);
             mf.ParallelAdd(jf, 0, 0, 1, ng, IntVect::TheZeroVector(), period);
-            jf.SumBoundary(period);
-            MultiFab::Copy(*current_cp[lev+1][idim], jf, 0, 0, 1, 0);
+            WarpXSumGuardCells(*current_cp[lev+1][idim], jf, period);
         }
         else if (current_buf[lev+1][idim]) // but no filter
         {
@@ -803,14 +804,14 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev)
             mf.ParallelAdd(*current_buf[lev+1][idim], 0, 0, 1,
                            current_buf[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(),
                            period);
-            current_cp[lev+1][idim]->SumBoundary(period);
+            WarpXSumGuardCells(*(current_cp[lev+1][idim]), period);
         }
         else // no filter, no buffer
         {
             mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, 1,
                            current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(),
                            period);
-            current_cp[lev+1][idim]->SumBoundary(period);
+            WarpXSumGuardCells(*(current_cp[lev+1][idim]), period);
         }
         MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, 1, 0);
     }
@@ -840,10 +841,9 @@ WarpX::ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, i
         ng += bilinear_filter.stencil_length_each_dir-1;
         MultiFab rf(r->boxArray(), r->DistributionMap(), ncomp, ng);
         bilinear_filter.ApplyStencil(rf, *r, icomp, 0, ncomp);
-        rf.SumBoundary(period);
-        MultiFab::Copy(*r, rf, 0, icomp, ncomp, 0);
+        WarpXSumGuardCells(*r, rf, period, icomp, ncomp );
     } else {
-        r->SumBoundary(icomp, ncomp, period);
+        WarpXSumGuardCells(*r, period, icomp, ncomp);
     }
 }
 
@@ -885,9 +885,7 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp)
 
             MultiFab::Add(rhofb, rhofc, 0, 0, ncomp, ng);
             mf.ParallelAdd(rhofb, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period);
-
-            rhofc.SumBoundary(period);
-            MultiFab::Copy(*rho_cp[lev+1], rhofc, 0, 0, ncomp, 0);
+            WarpXSumGuardCells( *rho_cp[lev+1], rhofc, period, icomp, ncomp );
         }
         else if (use_filter) // but no buffer
         {
@@ -896,8 +894,7 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp)
             MultiFab rf(rho_cp[lev+1]->boxArray(), rho_cp[lev+1]->DistributionMap(), ncomp, ng);
             bilinear_filter.ApplyStencil(rf, *rho_cp[lev+1], icomp, 0, ncomp);
             mf.ParallelAdd(rf, 0, 0, ncomp, ng, IntVect::TheZeroVector(), period);
-            rf.SumBoundary(0, ncomp, period);
-            MultiFab::Copy(*rho_cp[lev+1], rf, 0, icomp, ncomp, 0);
+            WarpXSumGuardCells( *rho_cp[lev+1], rf, period, icomp, ncomp );
         }
         else if (charge_buf[lev+1]) // but no filter
         {
@@ -908,14 +905,14 @@ WarpX::AddRhoFromFineLevelandSumBoundary(int lev, int icomp, int ncomp)
                            ncomp,
                            charge_buf[lev+1]->nGrowVect(), IntVect::TheZeroVector(),
                            period);
-            rho_cp[lev+1]->SumBoundary(icomp, ncomp, period);
+            WarpXSumGuardCells(*(rho_cp[lev+1]), period, icomp, ncomp);
         }
         else // no filter, no buffer
         {
             mf.ParallelAdd(*rho_cp[lev+1], icomp, 0, ncomp,
                            rho_cp[lev+1]->nGrowVect(), IntVect::TheZeroVector(),
                            period);
-            rho_cp[lev+1]->SumBoundary(icomp, ncomp, period);
+            WarpXSumGuardCells(*(rho_cp[lev+1]), period, icomp, ncomp);
         }
         ApplyFilterandSumBoundaryRho(lev, PatchType::fine, icomp, ncomp);
         MultiFab::Add(*rho_fp[lev], mf, 0, icomp, ncomp, 0);
diff --git a/Source/Parallelization/WarpXSumGuardCells.H b/Source/Parallelization/WarpXSumGuardCells.H
new file mode 100644
index 000000000..dfb255ab8
--- /dev/null
+++ b/Source/Parallelization/WarpXSumGuardCells.H
@@ -0,0 +1,66 @@
+#ifndef WARPX_SUM_GUARD_CELLS_H_
+#define WARPX_SUM_GUARD_CELLS_H_
+
+#include <AMReX_MultiFab.H>
+
+/* \brief Sum the values of `mf`, where the different boxes overlap
+ * (i.e. in the guard cells)
+ *
+ * This is typically called for the sources of the Maxwell equations (J/rho)
+ * after deposition from the macroparticles.
+ *
+ *  - When WarpX is compiled with a finite-difference scheme: this only
+ *    updates the *valid* cells of `mf`
+ *  - When WarpX is compiled with a spectral scheme (WARPX_USE_PSATD): this
+ *    updates both the *valid* cells and *guard* cells. (This is because a
+ *    spectral solver requires the value of the sources over a large stencil.)
+ */
+void
+WarpXSumGuardCells(amrex::MultiFab& mf, const amrex::Periodicity& period,
+                   const int icomp=0, const int ncomp=1){
+#ifdef WARPX_USE_PSATD
+    // Update both valid cells and guard cells
+    const int n_grow = mf.nGrow();
+    // Create temporary MultiFab, for safe summation
+    amrex::MultiFab tmp(mf.boxArray(), mf.DistributionMap(), ncomp, n_grow);
+    amrex::MultiFab::Copy(tmp, mf, icomp, 0, ncomp, n_grow);
+    mf.setVal(0.0, icomp, ncomp, 0);
+    // Perform summation on overlapping grid points
+    mf.copy(tmp, 0, icomp, ncomp, n_grow, n_grow, period, amrex::FabArrayBase::ADD);
+#else
+    // Update only the valid cells
+    mf.SumBoundary(icomp, ncomp, period);
+#endif
+}
+
+/* \brief Sum the values of `src` where the different boxes overlap
+ * (i.e. in the guard cells) and copy them into `dst`
+ *
+ * This is typically called for the sources of the Maxwell equations (J/rho)
+ * after deposition from the macroparticles + filtering.
+ *
+ *  - When WarpX is compiled with a finite-difference scheme: this only
+ *    updates the *valid* cells of `dst`
+ *  - When WarpX is compiled with a spectral scheme (WARPX_USE_PSATD): this
+ *    updates both the *valid* cells and *guard* cells. (This is because a
+ *    spectral solver requires the value of the sources over a large stencil.)
+ *
+ * Note: `i_comp` is the component where the results will be stored in `dst`;
+ *       The component from which we copy in `src` is always 0.
+ */
+void
+WarpXSumGuardCells(amrex::MultiFab& dst, amrex::MultiFab& src,
+                   const amrex::Periodicity& period,
+                   const int icomp=0, const int ncomp=1){
+#ifdef WARPX_USE_PSATD
+    // Update both valid cells and guard cells
+    dst.ParallelCopy( src, 0, icomp, ncomp,
+            src.nGrow(), dst.nGrow(), period, amrex::FabArrayBase::ADD );
+#else
+    // Update only the valid cells
+    src.SumBoundary(icomp, ncomp, period);
+    amrex::MultiFab::Copy( dst, src, 0, icomp, ncomp, 0 );
+#endif
+}
+
+#endif // WARPX_SUM_GUARD_CELLS_H_