From d9faf00fe9ee77533454e2394329d5277f5b587d Mon Sep 17 00:00:00 2001
From: Edoardo Zoni <ezoni@lbl.gov>
Date: Thu, 24 Feb 2022 11:16:06 -0800
Subject: [PATCH] Extend `do_pml_Lo/Hi` to MR Levels

---
 Source/Initialization/WarpXInitData.cpp | 23 ++++++++++-------------
 Source/WarpX.H                          |  4 ++--
 Source/WarpX.cpp                        |  8 ++++++--
 3 files changed, 18 insertions(+), 17 deletions(-)

diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 694518b4064..ec321d46ba6 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -228,24 +228,21 @@ WarpX::InitFromScratch ()
 void
 WarpX::InitPML ()
 {
-
     for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
         if (WarpX::field_boundary_lo[idim] == FieldBoundaryType::PML) {
             do_pml = 1;
-            do_pml_Lo[idim] = 1;
+            do_pml_Lo[0][idim] = 1; // on level 0
         }
         if (WarpX::field_boundary_hi[idim] == FieldBoundaryType::PML) {
             do_pml = 1;
-            do_pml_Hi[idim] = 1;
+            do_pml_Hi[0][idim] = 1; // on level 0
         }
     }
     if (finest_level > 0) do_pml = 1;
     if (do_pml)
     {
-        amrex::IntVect do_pml_Lo_corrected = do_pml_Lo;
-
 #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
-        do_pml_Lo_corrected[0] = 0; // no PML at r=0, in cylindrical geometry
+        do_pml_Lo[0][0] = 0; // no PML at r=0, in cylindrical geometry
         pml_rz[0] = std::make_unique<PML_RZ>(0, boxArray(0), DistributionMap(0), &Geom(0), pml_ncell, do_pml_in_domain);
 #else
         pml[0] = std::make_unique<PML>(0, boxArray(0), DistributionMap(0), &Geom(0), nullptr,
@@ -255,28 +252,28 @@ WarpX::InitPML ()
                              do_multi_J,
                              do_pml_dive_cleaning, do_pml_divb_cleaning,
                              guard_cells.ng_FieldSolver.max(),
-                             do_pml_Lo_corrected, do_pml_Hi);
+                             do_pml_Lo[0], do_pml_Hi[0]);
 #endif
 
         for (int lev = 1; lev <= finest_level; ++lev)
         {
-            amrex::IntVect do_pml_Lo_MR = amrex::IntVect::TheUnitVector();
-            amrex::IntVect do_pml_Hi_MR = amrex::IntVect::TheUnitVector();
+            do_pml_Lo[lev] = amrex::IntVect::TheUnitVector();
+            do_pml_Hi[lev] = amrex::IntVect::TheUnitVector();
             // check if fine patch edges co-incide with domain boundary
             amrex::Box levelBox = boxArray(lev).minimalBox();
             // Domain box at level, lev
             amrex::Box DomainBox = Geom(lev).Domain();
             for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                 if (levelBox.smallEnd(idim) == DomainBox.smallEnd(idim))
-                    do_pml_Lo_MR[idim] = do_pml_Lo[idim];
+                    do_pml_Lo[lev][idim] = do_pml_Lo[0][idim];
                 if (levelBox.bigEnd(idim) == DomainBox.bigEnd(idim))
-                    do_pml_Hi_MR[idim] = do_pml_Hi[idim];
+                    do_pml_Hi[lev][idim] = do_pml_Hi[0][idim];
             }
 
 #ifdef WARPX_DIM_RZ
             //In cylindrical geometry, if the edge of the patch is at r=0, do not add PML
             if ((max_level > 0) && (fine_tag_lo[0]==0.)) {
-                do_pml_Lo_MR[0] = 0;
+                do_pml_Lo[lev][0] = 0;
             }
 #endif
             pml[lev] = std::make_unique<PML>(lev, boxArray(lev), DistributionMap(lev),
@@ -286,7 +283,7 @@ WarpX::InitPML ()
                                    do_moving_window, pml_has_particles, do_pml_in_domain,
                                    do_multi_J, do_pml_dive_cleaning, do_pml_divb_cleaning,
                                    guard_cells.ng_FieldSolver.max(),
-                                   do_pml_Lo_MR, do_pml_Hi_MR);
+                                   do_pml_Lo[lev], do_pml_Hi[lev]);
         }
     }
 }
diff --git a/Source/WarpX.H b/Source/WarpX.H
index b74598bf14f..c72d8221eb4 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -1136,8 +1136,8 @@ private:
     static int do_similar_dm_pml;
     bool do_pml_dive_cleaning; // default set in WarpX.cpp
     bool do_pml_divb_cleaning; // default set in WarpX.cpp
-    amrex::IntVect do_pml_Lo = amrex::IntVect::TheZeroVector();
-    amrex::IntVect do_pml_Hi = amrex::IntVect::TheZeroVector();
+    amrex::Vector<amrex::IntVect> do_pml_Lo;
+    amrex::Vector<amrex::IntVect> do_pml_Hi;
     amrex::Vector<std::unique_ptr<PML> > pml;
 #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
     amrex::Vector<std::unique_ptr<PML_RZ> > pml_rz;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 971b1354876..88137fe71e5 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -324,6 +324,10 @@ WarpX::WarpX ()
 #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
     pml_rz.resize(nlevs_max);
 #endif
+
+    do_pml_Lo.resize(nlevs_max);
+    do_pml_Hi.resize(nlevs_max);
+
     costs.resize(nlevs_max);
     load_balance_efficiency.resize(nlevs_max);
 
@@ -2326,8 +2330,8 @@ WarpX::getPMLdirections() const
     {
         for( int i = 0; i < static_cast<int>(dirsWithPML.size()) / 2; ++i )
         {
-            dirsWithPML.at( 2u*i      ) = bool(do_pml_Lo[i]);
-            dirsWithPML.at( 2u*i + 1u ) = bool(do_pml_Hi[i]);
+            dirsWithPML.at( 2u*i      ) = bool(do_pml_Lo[0][i]); // on level 0
+            dirsWithPML.at( 2u*i + 1u ) = bool(do_pml_Hi[0][i]); // on level 0
         }
     }
     return dirsWithPML;