Adding support for pml in domain for multiple patches

Source/BoundaryConditions/PML.cpp

@@ -565,35 +565,43 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
     // In order to implement this, a reduced domain is created here (decreased by ncells in all direction)
     // and passed to `MakeBoxArray`, which surrounds it by PML boxes
     // (thus creating the PML boxes at the right position, where they overlap with the original domain)
-    // minimalBox provides the bounding box around grid_ba for level, lev.
-    // 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();
+
+    BoxArray grid_ba_reduced = grid_ba;
     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);
+        BoxList bl = grid_ba.simplified_list();
+        for (auto& b : bl) {
+            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+                if (do_pml_Lo[idim]) {
+                    Box const& bb = amrex::adjCellLo(b, idim);
+                    if ( ! grid_ba.intersects(bb) ) {
+                        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(b.length(idim) > ncell, " box length must be greater that pml size");
+                        b.growLo(idim, -ncell);
+                    }
+                }
+                if (do_pml_Hi[idim]) {
+                    Box const& bb = amrex::adjCellHi(b, idim);
+                    if ( ! grid_ba.intersects(bb) ) {
+                        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(b.length(idim) > ncell, " box length must be greater that pml size");
+                        b.growHi(idim, -ncell);
+                    }
+                }
             }
         }
+        grid_ba_reduced = BoxArray(std::move(bl));
     }
-    const BoxArray grid_ba_reduced = (do_pml_in_domain) ?
-        BoxArray(grid_ba.boxList().intersect(domain0)) : grid_ba;
-
+    Box const domain0 = grid_ba_reduced.minimalBox();
     const bool is_single_box_domain = domain0.numPts() == grid_ba_reduced.numPts();
     const BoxArray& ba = MakeBoxArray(is_single_box_domain, domain0, *geom, grid_ba_reduced,
                                       IntVect(ncell), do_pml_in_domain, do_pml_Lo, do_pml_Hi);
 
+
     if (ba.empty()) {
         m_ok = false;
         return;
     } else {
         m_ok = true;
     }
-
-    // Define the number of guard cells in each direction, for E, B, and F
+    // Define the number of guard cells in each di;rection, for E, B, and F
     auto nge = IntVect(AMREX_D_DECL(2, 2, 2));
     auto ngb = IntVect(AMREX_D_DECL(2, 2, 2));
     int ngf_int = 0;
@@ -760,24 +768,28 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
         BoxArray grid_cba = grid_ba;
         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 cdomain = grid_cba.minimalBox();
+        BoxArray grid_cba_reduced = grid_cba;
         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]);
+            BoxList bl = grid_cba.simplified_list();
+            for (auto& b : bl) {
+                for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+                    if (do_pml_Lo[idim]) {
+                        Box const& bb = amrex::adjCellLo(b, idim);
+                        if ( ! grid_cba.intersects(bb) ) {
+                            b.growLo(idim, -ncell/ref_ratio[idim]);
+                        }
+                    }
+                    if (do_pml_Hi[idim]) {
+                        Box const& bb = amrex::adjCellHi(b, idim);
+                        if ( ! grid_cba.intersects(bb) ) {
+                            b.growHi(idim, -ncell/ref_ratio[idim]);
+                        }
+                    }
                 }
             }
+            grid_cba_reduced = BoxArray(std::move(bl));
         }
-        const BoxArray grid_cba_reduced = (do_pml_in_domain) ?
-            BoxArray(grid_cba.boxList().intersect(cdomain)) : grid_cba;
+        Box const cdomain = grid_cba_reduced.minimalBox();
 
         const IntVect cncells = IntVect(ncell)/ref_ratio;
         const IntVect cdelta = IntVect(delta)/ref_ratio;
@@ -792,7 +804,6 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
         } else {
             cdm.define(cba);
         }
-
         const amrex::BoxArray cba_Ex = amrex::convert(cba, WarpX::GetInstance().getEfield_cp(1,0).ixType().toIntVect());
         const amrex::BoxArray cba_Ey = amrex::convert(cba, WarpX::GetInstance().getEfield_cp(1,1).ixType().toIntVect());
         const amrex::BoxArray cba_Ez = amrex::convert(cba, WarpX::GetInstance().getEfield_cp(1,2).ixType().toIntVect());