More efficient PML BoxArray (#2631)

* More efficient PML BoxArray

If the union of the grids is a single rectangular domain, we can simplify
the process and generate more efficient PML BoxArray.

* Update Source/BoundaryConditions/PML.cpp

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Apply suggestions from code review

* reset Python_wrappers benchmark

* fix the computation of sigmas for the new BoxArray

* Revert "reset Python_wrappers benchmark"

This reverts commit 2999304571d525076ba6c1f7cbbcef6e9acafabb.

* fix warning

* fix 1d

* initialize to quiet NaN

* Reset Benchmark: pml_x_psatd

- maximum relative error: 2.50e-06
- new implementation: 10 PML grids
- old implementation: 24 PML grids

* Reset Benchmark: LaserAccelerationMR

- maximum relative error: 2.73e-04
- new implementation: (18,8,8) PML grids
- old implementation: (48,18,18) PML grids

* Reset Benchmark: LaserOnFine

- maximum relative error: 6.44e-05
- new implementation: (2,6,6) PML grids
- old implementation: (2,12,12) PML grids

* Reset Benchmark: PlasmaAccelerationMR

- maximum relative error: 6.84e-04
- new implementation: (10,6,6) PML grids
- old implementation: (24,12,12) PML grids

* Reset Benchmark: RefinedInjection

- maximum relative error: 2.55e-04
- new implementation: (18,8,8) PML grids
- old implementation: (48,18,18) PML grids

* Reset Benchmark: momentum-conserving-gather

- maximum relative error: 7.43e-04
- new implementation: (10,6,6) PML grids
- old implementation: (24,12,12) PML grids

* Reset Benchmark: subcyclingMR

- maximum relative error: 2.41e-05
- new implementation: (6,6,6) PML grids
- old implementation: (12,12,12) PML grids

* Reset Benchmark: Langmuir_multi_2d_MR

- maximum relative error: 1.32e-01 (B numerical artifact)
- new implementation: (0,20,20) PML grids
- old implementation: (0,52,40) PML grids

* Reset Benchmark: Langmuir_multi_2d_MR_psatd

- maximum relative error: 1.05e-01 (B numerical artifact)
- new implementation: (0,20,20) PML grids
- old implementation: (0,52,40) PML grids

* Reset Benchmark: Python_LaserAccelerationMR

- maximum relative error: 2.73e-04
- new implementation: (18,8,8) PML grids
- old implementation: (48,18,18) PML grids

* Reset Benchmark: Python_wrappers

- maximum relative error: 1.07e-08
- new implementation: 8 PML grids
- old implementation: 16 PML grids

* Reset Benchmark: pml_psatd_dive_divb_cleaning

- maximum relative error: 4.91e-03
- new implementation: 24 PML grids
- old implementation: 98 PML grids

* Remove an assertion.  We will fix it later

* Reset Benchmark: Langmuir_multi_2d_MR_anisotropic

- maximum relative error: 1.07e-01 (B numerical artifact)
- new implementation: (0,16,16) PML grids
- old implementation: (0,40,34) PML grids

* Reset Benchmark: PEC_field_mr

- maximum relative error: 3.98e-02
- new implementation: (0,2,2) PML grids
- old implementation: (0,2,2) PML grids
  (different number of ghost cells on coarse PML patch)

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
Co-authored-by: Edoardo Zoni <ezoni@lbl.gov>

1 files changed, 202 insertions, 91 deletions

diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp
index e6af648cc..352166697 100644
--- a/Source/BoundaryConditions/PML.cpp
+++ b/Source/BoundaryConditions/PML.cpp
@@ -44,6 +44,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <limits>
 #include <memory>
 #include <utility>
 #ifdef AMREX_USE_EB
@@ -54,15 +55,12 @@ using namespace amrex;
 
 namespace
 {
-    static void FillLo (int idim, Sigma& sigma, Sigma& sigma_cumsum,
+    static void FillLo (Sigma& sigma, Sigma& sigma_cumsum,
                         Sigma& sigma_star, Sigma& sigma_star_cumsum,
-                        const Box& overlap, const Box& grid, Real fac)
+                        const int olo, const int ohi, const int glo, Real fac)
     {
-        int glo = grid.smallEnd(idim);
-        int olo = overlap.smallEnd(idim);
-        int ohi = overlap.bigEnd(idim);
-        int slo = sigma.m_lo;
-        int sslo = sigma_star.m_lo;
+        const int slo = sigma.m_lo;
+        const int sslo = sigma_star.m_lo;
 
         const int N = ohi+1-olo+1;
         Real* p_sigma = sigma.data();
@@ -85,15 +83,12 @@ namespace
         });
     }
 
-    static void FillHi (int idim, Sigma& sigma, Sigma& sigma_cumsum,
+    static void FillHi (Sigma& sigma, Sigma& sigma_cumsum,
                         Sigma& sigma_star, Sigma& sigma_star_cumsum,
-                        const Box& overlap, const Box& grid, Real fac)
+                        const int olo, const int ohi, const int ghi, Real fac)
     {
-        int ghi = grid.bigEnd(idim);
-        int olo = overlap.smallEnd(idim);
-        int ohi = overlap.bigEnd(idim);
-        int slo = sigma.m_lo;
-        int sslo = sigma_star.m_lo;
+        const int slo = sigma.m_lo;
+        const int sslo = sigma_star.m_lo;
 
         const int N = ohi+1-olo+1;
         Real* p_sigma = sigma.data();
@@ -115,14 +110,12 @@ namespace
     }
 
 #if (AMREX_SPACEDIM != 1)
-    static void FillZero (int idim, Sigma& sigma, Sigma& sigma_cumsum,
+    static void FillZero (Sigma& sigma, Sigma& sigma_cumsum,
                           Sigma& sigma_star, Sigma& sigma_star_cumsum,
-                          const Box& overlap)
+                          const int olo, const int ohi)
     {
-        int olo = overlap.smallEnd(idim);
-        int ohi = overlap.bigEnd(idim);
-        int slo = sigma.m_lo;
-        int sslo = sigma_star.m_lo;
+        const int slo = sigma.m_lo;
+        const int sslo = sigma_star.m_lo;
 
         const int N = ohi+1-olo+1;
         Real* p_sigma = sigma.data();
@@ -144,7 +137,8 @@ namespace
 }
 
 
-SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int ncell, int delta)
+SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int ncell, int delta,
+                    const amrex::Box& regdomain)
 {
     BL_ASSERT(box.cellCentered());
 
@@ -154,14 +148,14 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
 
     for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
     {
-        sigma                [idim].resize(sz[idim]+1);
-        sigma_cumsum         [idim].resize(sz[idim]+1);
-        sigma_star           [idim].resize(sz[idim]+1);
-        sigma_star_cumsum    [idim].resize(sz[idim]+1);
-        sigma_fac            [idim].resize(sz[idim]+1);
-        sigma_cumsum_fac     [idim].resize(sz[idim]+1);
-        sigma_star_fac       [idim].resize(sz[idim]+1);
-        sigma_star_cumsum_fac[idim].resize(sz[idim]+1);
+        sigma                [idim].resize(sz[idim]+1,std::numeric_limits<Real>::quiet_NaN());
+        sigma_cumsum         [idim].resize(sz[idim]+1,std::numeric_limits<Real>::quiet_NaN());
+        sigma_star           [idim].resize(sz[idim]+1,std::numeric_limits<Real>::quiet_NaN());
+        sigma_star_cumsum    [idim].resize(sz[idim]+1,std::numeric_limits<Real>::quiet_NaN());
+        sigma_fac            [idim].resize(sz[idim]+1,std::numeric_limits<Real>::quiet_NaN());
+        sigma_cumsum_fac     [idim].resize(sz[idim]+1,std::numeric_limits<Real>::quiet_NaN());
+        sigma_star_fac       [idim].resize(sz[idim]+1,std::numeric_limits<Real>::quiet_NaN());
+        sigma_star_cumsum_fac[idim].resize(sz[idim]+1,std::numeric_limits<Real>::quiet_NaN());
 
         sigma                [idim].m_lo = lo[idim];
         sigma                [idim].m_hi = hi[idim]+1;
@@ -186,6 +180,58 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
         fac[idim] = 4.0_rt*PhysConst::c/(dx[idim]*static_cast<Real>(delta*delta));
     }
 
+    if (regdomain.ok()) { // The union of the regular grids is a single box
+        define_single(regdomain, ncell, fac);
+    } else {
+        define_multiple(box, grids, ncell, fac);
+    }
+}
+
+void SigmaBox::define_single (const Box& regdomain, int ncell,
+                              const Array<Real,AMREX_SPACEDIM>& fac)
+{
+    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+        const int slo = sigma[idim].lo();
+        const int shi = sigma[idim].hi()-1;
+        const int dlo = regdomain.smallEnd(idim);
+        const int dhi = regdomain.bigEnd(idim);
+
+        // Lo
+        int olo = std::max(slo, dlo-ncell);
+        int ohi = std::min(shi, dlo-1);
+        if (ohi >= olo) {
+            FillLo(sigma[idim], sigma_cumsum[idim],
+                   sigma_star[idim], sigma_star_cumsum[idim],
+                   olo, ohi, dlo, fac[idim]);
+        }
+
+#if (AMREX_SPACEDIM != 1)
+        // Mid
+        olo = std::max(slo, dlo);
+        ohi = std::min(shi, dhi);
+        if (ohi >= olo) {
+            FillZero(sigma[idim], sigma_cumsum[idim],
+                     sigma_star[idim], sigma_star_cumsum[idim],
+                     olo, ohi);
+        }
+#endif
+
+        // Hi
+        olo = std::max(slo, dhi+1);
+        ohi = std::min(shi, dhi+ncell);
+        if (ohi >= olo) {
+            FillHi(sigma[idim], sigma_cumsum[idim],
+                   sigma_star[idim], sigma_star_cumsum[idim],
+                   olo, ohi, dhi, fac[idim]);
+        }
+    }
+
+    amrex::Gpu::streamSynchronize();
+}
+
+void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, int ncell,
+                                const Array<Real,AMREX_SPACEDIM>& fac)
+{
     const std::vector<std::pair<int,Box> >& isects = grids.intersections(box, false, ncell);
 
     for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
@@ -251,9 +297,10 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
 #endif
             Box looverlap = lobox & box;
             if (looverlap.ok()) {
-                FillLo(idim, sigma[idim], sigma_cumsum[idim],
+                FillLo(sigma[idim], sigma_cumsum[idim],
                        sigma_star[idim], sigma_star_cumsum[idim],
-                       looverlap, grid_box, fac[idim]);
+                       looverlap.smallEnd(idim), looverlap.bigEnd(idim),
+                       grid_box.smallEnd(idim), fac[idim]);
             }
 
             Box hibox = amrex::adjCellHi(grid_box, idim, ncell);
@@ -263,9 +310,10 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
 #endif
             Box hioverlap = hibox & box;
             if (hioverlap.ok()) {
-                FillHi(idim, sigma[idim], sigma_cumsum[idim],
+                FillHi(sigma[idim], sigma_cumsum[idim],
                        sigma_star[idim],  sigma_star_cumsum[idim],
-                       hioverlap, grid_box, fac[idim]);
+                       hioverlap.smallEnd(idim), hioverlap.bigEnd(idim),
+                       grid_box.bigEnd(idim), fac[idim]);
             }
 
             if (!looverlap.ok() && !hioverlap.ok()) {
@@ -280,8 +328,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
             const Box& grid_box = grids[gid];
             const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box;
             if (overlap.ok()) {
-                FillZero(idim, sigma[idim], sigma_cumsum[idim],
-                        sigma_star[idim], sigma_star_cumsum[idim], overlap);
+                FillZero(sigma[idim], sigma_cumsum[idim],
+                         sigma_star[idim], sigma_star_cumsum[idim],
+                         overlap.smallEnd(idim), overlap.bigEnd(idim));
             }
             else {
                 amrex::Abort("SigmaBox::SigmaBox(): side_side_edges, how did this happen?\n");
@@ -295,17 +344,19 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
             Box lobox = amrex::adjCellLo(grid_box, idim, ncell);
             Box looverlap = lobox.grow(jdim,ncell).grow(kdim,ncell) & box;
             if (looverlap.ok()) {
-                FillLo(idim, sigma[idim], sigma_cumsum[idim],
-                      sigma_star[idim],  sigma_star_cumsum[idim],
-                      looverlap, grid_box, fac[idim]);
+                FillLo(sigma[idim], sigma_cumsum[idim],
+                       sigma_star[idim],  sigma_star_cumsum[idim],
+                       looverlap.smallEnd(idim), looverlap.bigEnd(idim),
+                       grid_box.smallEnd(idim), fac[idim]);
             }
 
             Box hibox = amrex::adjCellHi(grid_box, idim, ncell);
             Box hioverlap = hibox.grow(jdim,ncell).grow(kdim,ncell) & box;
             if (hioverlap.ok()) {
-                FillHi(idim, sigma[idim], sigma_cumsum[idim],
-                      sigma_star[idim],  sigma_star_cumsum[idim],
-                      hioverlap, grid_box, fac[idim]);
+                FillHi(sigma[idim], sigma_cumsum[idim],
+                       sigma_star[idim], sigma_star_cumsum[idim],
+                       hioverlap.smallEnd(idim), hioverlap.bigEnd(idim),
+                       grid_box.bigEnd(idim), fac[idim]);
             }
 
             if (!looverlap.ok() && !hioverlap.ok()) {
@@ -324,8 +375,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
             const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box;
 #endif
             if (overlap.ok()) {
-                FillZero(idim, sigma[idim], sigma_cumsum[idim],
-                        sigma_star[idim], sigma_star_cumsum[idim], overlap);
+                FillZero(sigma[idim], sigma_cumsum[idim],
+                         sigma_star[idim], sigma_star_cumsum[idim],
+                         overlap.smallEnd(idim), overlap.bigEnd(idim));
             } else {
                 amrex::Abort("SigmaBox::SigmaBox(): side_faces, how did this happen?\n");
             }
@@ -339,17 +391,19 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
             const Box& lobox = amrex::adjCellLo(grid_box, idim, ncell);
             Box looverlap = lobox & box;
             if (looverlap.ok()) {
-                FillLo(idim, sigma[idim], sigma_cumsum[idim],
-                      sigma_star[idim],  sigma_star_cumsum[idim],
-                      looverlap, grid_box, fac[idim]);
+                FillLo(sigma[idim], sigma_cumsum[idim],
+                       sigma_star[idim], sigma_star_cumsum[idim],
+                       looverlap.smallEnd(idim), looverlap.bigEnd(idim),
+                       grid_box.smallEnd(idim), fac[idim]);
             }
 
             const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell);
             Box hioverlap = hibox & box;
             if (hioverlap.ok()) {
-                FillHi(idim, sigma[idim], sigma_cumsum[idim],
-                      sigma_star[idim],  sigma_star_cumsum[idim],
-                      hioverlap, grid_box, fac[idim]);
+                FillHi(sigma[idim], sigma_cumsum[idim],
+                       sigma_star[idim], sigma_star_cumsum[idim],
+                       hioverlap.smallEnd(idim), hioverlap.bigEnd(idim),
+                       grid_box.bigEnd(idim), fac[idim]);
             }
 
             if (!looverlap.ok() && !hioverlap.ok()) {
@@ -362,7 +416,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
         }
     }
 
-    amrex::Gpu::synchronize();
+    amrex::Gpu::streamSynchronize();
 }
 
 
@@ -435,9 +489,10 @@ SigmaBox::ComputePMLFactorsE (const Real* a_dx, Real dt)
 }
 
 MultiSigmaBox::MultiSigmaBox (const BoxArray& ba, const DistributionMapping& dm,
-                              const BoxArray& grid_ba, const Real* dx, int ncell, int delta)
+                              const BoxArray& grid_ba, const Real* dx, int ncell, int delta,
+                              const amrex::Box& regular_domain)
     : FabArray<SigmaBox>(ba,dm,1,0,MFInfo(),
-                         SigmaBoxFactory(grid_ba,dx,ncell,delta))
+                         SigmaBoxFactory(grid_ba,dx,ncell,delta, regular_domain))
 {}
 
 void
@@ -495,19 +550,23 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
     // Note that this is okay to build pml inside domain for a single patch, or joint patches
     // with same [min,max]. But it does not support multiple disjoint refinement patches.
     Box domain0 = grid_ba.minimalBox();
-    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-        if (do_pml_Lo[idim]){
-            domain0.growLo(idim, -ncell);
-        }
-        if (do_pml_Hi[idim]){
-            domain0.growHi(idim, -ncell);
+    if (do_pml_in_domain) {
+        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+            if (do_pml_Lo[idim]){
+                domain0.growLo(idim, -ncell);
+            }
+            if (do_pml_Hi[idim]){
+                domain0.growHi(idim, -ncell);
+            }
         }
     }
-    const BoxArray grid_ba_reduced = BoxArray(grid_ba.boxList().intersect(domain0));
+    const BoxArray grid_ba_reduced = (do_pml_in_domain) ?
+        BoxArray(grid_ba.boxList().intersect(domain0)) : grid_ba;
+
+    bool is_single_box_domain = domain0.numPts() == grid_ba_reduced.numPts();
+    const BoxArray& ba = MakeBoxArray(is_single_box_domain, domain0, *geom, grid_ba_reduced, ncell,
+                                      do_pml_in_domain, do_pml_Lo, do_pml_Hi);
 
-    const BoxArray& ba = (do_pml_in_domain)?
-          MakeBoxArray(*geom, grid_ba_reduced, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi) :
-          MakeBoxArray(*geom, grid_ba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi);
     if (ba.empty()) {
         m_ok = false;
         return;
@@ -660,12 +719,10 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
         pml_G_fp->setVal(0.0);
     }
 
-    if (do_pml_in_domain){
-        sigba_fp = std::make_unique<MultiSigmaBox>(ba, dm, grid_ba_reduced, geom->CellSize(), ncell, delta);
-    }
-    else {
-        sigba_fp = std::make_unique<MultiSigmaBox>(ba, dm, grid_ba, geom->CellSize(), ncell, delta);
-    }
+    Box single_domain_box = is_single_box_domain ? domain0 : Box();
+    // Empty box (i.e., Box()) means it's not a single box domain.
+    sigba_fp = std::make_unique<MultiSigmaBox>(ba, dm, grid_ba_reduced, geom->CellSize(),
+                                               ncell, delta, single_domain_box);
 
     if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
 #ifndef WARPX_USE_PSATD
@@ -705,26 +762,31 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
         grid_cba.coarsen(ref_ratio);
 
         // assuming that the bounding box around grid_cba is a single patch, and not disjoint patches, similar to fine patch.
-        amrex::Box domain1 = grid_cba.minimalBox();
-        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-            if (do_pml_Lo[idim]){
-                // ncell is divided by refinement ratio to ensure that the
-                // physical width of the PML region is equal in fine and coarse patch
-                domain1.growLo(idim, -ncell/ref_ratio[idim]);
-            }
-            if (do_pml_Hi[idim]){
-                // ncell is divided by refinement ratio to ensure that the
-                // physical width of the PML region is equal in fine and coarse patch
-                domain1.growHi(idim, -ncell/ref_ratio[idim]);
+        amrex::Box cdomain = grid_cba.minimalBox();
+        if (do_pml_in_domain) {
+            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+                if (do_pml_Lo[idim]){
+                    // ncell is divided by refinement ratio to ensure that the
+                    // physical width of the PML region is equal in fine and coarse patch
+                    cdomain.growLo(idim, -ncell/ref_ratio[idim]);
+                }
+                if (do_pml_Hi[idim]){
+                    // ncell is divided by refinement ratio to ensure that the
+                    // physical width of the PML region is equal in fine and coarse patch
+                    cdomain.growHi(idim, -ncell/ref_ratio[idim]);
+                }
             }
         }
-        const BoxArray grid_cba_reduced = BoxArray(grid_cba.boxList().intersect(domain1));
+        const BoxArray grid_cba_reduced = (do_pml_in_domain) ?
+            BoxArray(grid_cba.boxList().intersect(cdomain)) : grid_cba;
 
-        // Assuming that refinement ratio is equal in all dimensions
-        const BoxArray& cba = (do_pml_in_domain) ?
-            MakeBoxArray(*cgeom, grid_cba_reduced, ncell/ref_ratio[0], do_pml_in_domain, do_pml_Lo, do_pml_Hi) :
-            MakeBoxArray(*cgeom, grid_cba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi);
+// This needs to be fixed.        AMREX_ALWAYS_ASSERT(ref_ratio == ref_ratio[0]);
+        const int cncells = ncell/ref_ratio[0];
+        const int cdelta = delta/ref_ratio[0];
 
+        // Assuming that refinement ratio is equal in all dimensions
+        const BoxArray& cba = MakeBoxArray(is_single_box_domain, cdomain, *cgeom, grid_cba_reduced,
+                                           cncells, do_pml_in_domain, do_pml_Lo, do_pml_Hi);
         DistributionMapping cdm;
         if (WarpX::do_similar_dm_pml) {
             auto ng_sim = amrex::elemwiseMax(amrex::elemwiseMax(nge, ngb), ngf);
@@ -781,12 +843,9 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
         pml_j_cp[1]->setVal(0.0);
         pml_j_cp[2]->setVal(0.0);
 
-        if (do_pml_in_domain){
-            // Note - assuming that the refinement ratio is equal in all dimensions
-            sigba_cp = std::make_unique<MultiSigmaBox>(cba, cdm, grid_cba_reduced, cgeom->CellSize(), ncell/ref_ratio[0], delta/ref_ratio[0]);
-        } else {
-            sigba_cp = std::make_unique<MultiSigmaBox>(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta);
-        }
+        single_domain_box = is_single_box_domain ? cdomain : Box();
+        sigba_cp = std::make_unique<MultiSigmaBox>(cba, cdm, grid_cba_reduced, cgeom->CellSize(),
+                                                   cncells, cdelta, single_domain_box);
 
         if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
 #ifndef WARPX_USE_PSATD
@@ -815,10 +874,62 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
 }
 
 BoxArray
-PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba,
+PML::MakeBoxArray (bool is_single_box_domain, const amrex::Box& regular_domain,
+                   const amrex::Geometry& geom, const amrex::BoxArray& grid_ba,
                    int ncell, int do_pml_in_domain,
                    const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi)
 {
+    if (is_single_box_domain) {
+        return MakeBoxArray_single(regular_domain, grid_ba, ncell, do_pml_Lo, do_pml_Hi);
+    } else { // the union of the regular grids is *not* a single rectangular domain
+        return MakeBoxArray_multiple(geom, grid_ba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi);
+    }
+}
+
+BoxArray
+PML::MakeBoxArray_single (const amrex::Box& regular_domain, const amrex::BoxArray& grid_ba,
+                          int ncell, const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi)
+{
+    BoxList bl;
+    for (int i = 0, N = grid_ba.size(); i < N; ++i) {
+        Box const& b = grid_ba[i];
+        for (OrientationIter oit; oit.isValid(); ++oit) {
+            // In 3d, a Box has 6 faces.  This iterates over the 6 faces.
+            // 3 of them are on the lower side and the others are on the
+            // higher side.
+            Orientation ori = oit();
+            const int idim = ori.coordDir(); // either 0 or 1 or 2 (i.e., x, y, z-direction)
+            bool pml_bndry = false;
+            if (ori.isLow() && do_pml_Lo[idim]) {  // This is one of the lower side faces.
+                pml_bndry = b.smallEnd(idim) == regular_domain.smallEnd(idim);
+            } else if (ori.isHigh() && do_pml_Hi[idim]) { // This is one of the higher side faces.
+                pml_bndry = b.bigEnd(idim) == regular_domain.bigEnd(idim);
+            }
+            if (pml_bndry) {
+                Box bbox = amrex::adjCell(b, ori, ncell);
+                for (int jdim = 0; jdim < idim; ++jdim) {
+                    if (do_pml_Lo[jdim] &&
+                        bbox.smallEnd(jdim) == regular_domain.smallEnd(jdim)) {
+                        bbox.growLo(jdim, ncell);
+                    }
+                    if (do_pml_Hi[jdim] &&
+                        bbox.bigEnd(jdim) == regular_domain.bigEnd(jdim)) {
+                        bbox.growHi(jdim, ncell);
+                    }
+                }
+                bl.push_back(bbox);
+            }
+        }
+    }
+
+    return BoxArray(std::move(bl));
+}
+
+BoxArray
+PML::MakeBoxArray_multiple (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba,
+                            int ncell, int do_pml_in_domain,
+                            const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi)
+{
     Box domain = geom.Domain();
     for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
         if (do_pml_Lo[idim]){