Merge branch 'dev' into RZgeometry

1 files changed, 264 insertions, 82 deletions

diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp
index 21d348482..8f8a2608e 100644
--- a/Source/BoundaryConditions/PML.cpp
+++ b/Source/BoundaryConditions/PML.cpp
@@ -1,4 +1,3 @@
-
 #include <PML.H>
 #include <WarpX.H>
 #include <WarpXConst.H>
@@ -16,54 +15,70 @@ using namespace amrex;
 
 namespace
 {
-    static void FillLo (int idim, Sigma& sigma, Sigma& sigma_star,
+    static void FillLo (int idim, Sigma& sigma, Sigma& sigma_cumsum,
+                        Sigma& sigma_star, Sigma& sigma_star_cumsum,
                         const Box& overlap, const Box& grid, Real fac)
     {
         int glo = grid.smallEnd(idim);
         int olo = overlap.smallEnd(idim);
         int ohi = overlap.bigEnd(idim);
         int slo = sigma.m_lo;
+        int shi = sigma.m_hi;
         int sslo = sigma_star.m_lo;
+
         for (int i = olo; i <= ohi+1; ++i)
         {
             Real offset = static_cast<Real>(glo-i);
             sigma[i-slo] = fac*(offset*offset);
+            // sigma_cumsum is the analytical integral of sigma function at same points than sigma
+            sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3.)/PhysConst::c;
         }
+
         for (int i = olo; i <= ohi; ++i)
         {
             Real offset = static_cast<Real>(glo-i) - 0.5;
             sigma_star[i-sslo] = fac*(offset*offset);
+            // sigma_star_cumsum is the analytical integral of sigma function at same points than sigma_star
+            sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3.)/PhysConst::c;
         }
     }
 
-    static void FillHi (int idim, Sigma& sigma, Sigma& sigma_star,
+    static void FillHi (int idim, Sigma& sigma, Sigma& sigma_cumsum,
+                        Sigma& sigma_star, Sigma& sigma_star_cumsum,
                         const Box& overlap, const Box& grid, Real fac)
     {
         int ghi = grid.bigEnd(idim);
         int olo = overlap.smallEnd(idim);
         int ohi = overlap.bigEnd(idim);
         int slo = sigma.m_lo;
+        int shi = sigma.m_hi;
         int sslo = sigma_star.m_lo;
         for (int i = olo; i <= ohi+1; ++i)
         {
             Real offset = static_cast<Real>(i-ghi-1);
             sigma[i-slo] = fac*(offset*offset);
+            sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3.)/PhysConst::c;
         }
         for (int i = olo; i <= ohi; ++i)
         {
             Real offset = static_cast<Real>(i-ghi) - 0.5;
             sigma_star[i-sslo] = fac*(offset*offset);
+            sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3.)/PhysConst::c;
         }
     }
 
-    static void FillZero (int idim, Sigma& sigma, Sigma& sigma_star, const Box& overlap)
+    static void FillZero (int idim, Sigma& sigma, Sigma& sigma_cumsum,
+                          Sigma& sigma_star, Sigma& sigma_star_cumsum,
+                          const Box& overlap)
     {
         int olo = overlap.smallEnd(idim);
         int ohi = overlap.bigEnd(idim);
         int slo = sigma.m_lo;
         int sslo = sigma_star.m_lo;
         std::fill(sigma.begin()+(olo-slo), sigma.begin()+(ohi+2-slo), 0.0);
+        std::fill(sigma_cumsum.begin()+(olo-slo), sigma_cumsum.begin()+(ohi+2-slo), 0.0);
         std::fill(sigma_star.begin()+(olo-sslo), sigma_star.begin()+(ohi+1-sslo), 0.0);
+        std::fill(sigma_star_cumsum.begin()+(olo-sslo), sigma_star_cumsum.begin()+(ohi+1-sslo), 0.0);
     }
 }
 
@@ -77,19 +92,31 @@ 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_star    [idim].resize(sz[idim]  );
-        sigma_fac     [idim].resize(sz[idim]+1);
-        sigma_star_fac[idim].resize(sz[idim]  );
-
-        sigma         [idim].m_lo = lo[idim];
-        sigma         [idim].m_hi = hi[idim]+1;
-        sigma_star    [idim].m_lo = lo[idim];
-        sigma_star    [idim].m_hi = hi[idim];
-        sigma_fac     [idim].m_lo = lo[idim];
-        sigma_fac     [idim].m_hi = hi[idim]+1;
-        sigma_star_fac[idim].m_lo = lo[idim];
-        sigma_star_fac[idim].m_hi = hi[idim];
+        sigma                [idim].resize(sz[idim]+1);
+        sigma_cumsum         [idim].resize(sz[idim]+1);
+        sigma_star           [idim].resize(sz[idim]);
+        sigma_star_cumsum    [idim].resize(sz[idim]);
+        sigma_fac            [idim].resize(sz[idim]+1);
+        sigma_cumsum_fac     [idim].resize(sz[idim]+1);
+        sigma_star_fac       [idim].resize(sz[idim]);
+        sigma_star_cumsum_fac[idim].resize(sz[idim]);
+
+        sigma                [idim].m_lo = lo[idim];
+        sigma                [idim].m_hi = hi[idim]+1;
+        sigma_cumsum         [idim].m_lo = lo[idim];
+        sigma_cumsum         [idim].m_hi = hi[idim]+1;
+        sigma_star           [idim].m_lo = lo[idim];
+        sigma_star           [idim].m_hi = hi[idim];
+        sigma_star_cumsum    [idim].m_lo = lo[idim];
+        sigma_star_cumsum    [idim].m_hi = hi[idim];
+        sigma_fac            [idim].m_lo = lo[idim];
+        sigma_fac            [idim].m_hi = hi[idim]+1;
+        sigma_cumsum_fac     [idim].m_lo = lo[idim];
+        sigma_cumsum_fac     [idim].m_hi = hi[idim]+1;
+        sigma_star_fac       [idim].m_lo = lo[idim];
+        sigma_star_fac       [idim].m_hi = hi[idim];
+        sigma_star_cumsum_fac[idim].m_lo = lo[idim];
+        sigma_star_cumsum_fac[idim].m_hi = hi[idim];
     }
 
     Array<Real,AMREX_SPACEDIM> fac;
@@ -157,7 +184,9 @@ 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_star[idim], looverlap, grid_box, fac[idim]);
+                FillLo(idim, sigma[idim], sigma_cumsum[idim],
+                       sigma_star[idim], sigma_star_cumsum[idim],
+                       looverlap, grid_box, fac[idim]);
             }
 
             Box hibox = amrex::adjCellHi(grid_box, idim, ncell);
@@ -167,7 +196,9 @@ 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_star[idim], hioverlap, grid_box, fac[idim]);
+                FillHi(idim, sigma[idim], sigma_cumsum[idim],
+                       sigma_star[idim],  sigma_star_cumsum[idim],
+                       hioverlap, grid_box, fac[idim]);
             }
 
             if (!looverlap.ok() && !hioverlap.ok()) {
@@ -181,8 +212,10 @@ 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_star[idim], overlap);
-            } else {
+                FillZero(idim, sigma[idim], sigma_cumsum[idim],
+                        sigma_star[idim], sigma_star_cumsum[idim], overlap);
+            }
+            else {
                 amrex::Abort("SigmaBox::SigmaBox(): side_side_edges, how did this happen?\n");
             }
         }
@@ -194,13 +227,17 @@ 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_star[idim], looverlap, grid_box, fac[idim]);
+                FillLo(idim, sigma[idim], sigma_cumsum[idim],
+                      sigma_star[idim],  sigma_star_cumsum[idim],
+                      looverlap, grid_box, 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_star[idim], hioverlap, grid_box, fac[idim]);
+                FillHi(idim, sigma[idim], sigma_cumsum[idim],
+                      sigma_star[idim],  sigma_star_cumsum[idim],
+                      hioverlap, grid_box, fac[idim]);
             }
 
             if (!looverlap.ok() && !hioverlap.ok()) {
@@ -218,7 +255,8 @@ 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_star[idim], overlap);
+                FillZero(idim, sigma[idim], sigma_cumsum[idim],
+                        sigma_star[idim], sigma_star_cumsum[idim], overlap);
             } else {
                 amrex::Abort("SigmaBox::SigmaBox(): side_faces, how did this happen?\n");
             }
@@ -231,13 +269,17 @@ 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_star[idim], looverlap, grid_box, fac[idim]);
+                FillLo(idim, sigma[idim], sigma_cumsum[idim],
+                      sigma_star[idim],  sigma_star_cumsum[idim],
+                      looverlap, grid_box, fac[idim]);
             }
 
             const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell);
             Box hioverlap = hibox & box;
             if (hioverlap.ok()) {
-                FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]);
+                FillHi(idim, sigma[idim], sigma_cumsum[idim],
+                      sigma_star[idim],  sigma_star_cumsum[idim],
+                      hioverlap, grid_box, fac[idim]);
             }
 
             if (!looverlap.ok() && !hioverlap.ok()) {
@@ -251,6 +293,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
     }
 }
 
+
 void
 SigmaBox::ComputePMLFactorsB (const Real* dx, Real dt)
 {
@@ -259,6 +302,7 @@ SigmaBox::ComputePMLFactorsB (const Real* dx, Real dt)
         for (int i = 0, N = sigma_star[idim].size(); i < N; ++i)
         {
             sigma_star_fac[idim][i] = std::exp(-sigma_star[idim][i]*dt);
+            sigma_star_cumsum_fac[idim][i] = std::exp(-sigma_star_cumsum[idim][i]*dx[idim]);
         }
     }
 }
@@ -271,6 +315,7 @@ SigmaBox::ComputePMLFactorsE (const Real* dx, Real dt)
         for (int i = 0, N = sigma[idim].size(); i < N; ++i)
         {
             sigma_fac[idim][i] = std::exp(-sigma[idim][i]*dt);
+            sigma_cumsum_fac[idim][i] = std::exp(-sigma_cumsum[idim][i]*dx[idim]);
         }
     }
 }
@@ -320,11 +365,35 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
           Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal,
 #endif
           int do_dive_cleaning, int do_moving_window,
+          int pml_has_particles, int do_pml_in_domain,
           const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi)
     : m_geom(geom),
       m_cgeom(cgeom)
 {
-    const BoxArray& ba = MakeBoxArray(*geom, grid_ba, ncell, do_pml_Lo, do_pml_Hi);
+
+    // When `do_pml_in_domain` is true, the PML overlap with the last `ncell` of the physical domain
+    // (instead of extending `ncell` outside of the physical domain)
+    // 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)
+    Box domain0 = geom->Domain();
+    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+        if ( ! geom->isPeriodic(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& 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.size() == 0) {
         m_ok = false;
         return;
@@ -366,6 +435,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
     pml_B_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::By_nodal_flag), dm, 2, ngb));
     pml_B_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::Bz_nodal_flag), dm, 2, ngb));
 
+
     pml_E_fp[0]->setVal(0.0);
     pml_E_fp[1]->setVal(0.0);
     pml_E_fp[2]->setVal(0.0);
@@ -373,15 +443,30 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
     pml_B_fp[1]->setVal(0.0);
     pml_B_fp[2]->setVal(0.0);
 
+    pml_j_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::jx_nodal_flag), dm, 1, ngb));
+    pml_j_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::jy_nodal_flag), dm, 1, ngb));
+    pml_j_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::jz_nodal_flag), dm, 1, ngb));
+    pml_j_fp[0]->setVal(0.0);
+    pml_j_fp[1]->setVal(0.0);
+    pml_j_fp[2]->setVal(0.0);
+
     if (do_dive_cleaning)
     {
         pml_F_fp.reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()), dm, 3, ngf));
         pml_F_fp->setVal(0.0);
     }
 
-    sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta));
+    if (do_pml_in_domain){
+        sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba_reduced, geom->CellSize(), ncell, delta));
+    }
+    else {
+        sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta));
+    }
+
 
 #ifdef WARPX_USE_PSATD
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE( do_pml_in_domain==false,
+        "PSATD solver cannot be used with `do_pml_in_domain`.");
     const bool in_pml = true; // Tells spectral solver to use split-PML equations
     const RealVect dx{AMREX_D_DECL(geom->CellSize(0), geom->CellSize(1), geom->CellSize(2))};
     // Get the cell-centered box, with guard cells
@@ -400,7 +485,11 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
 
         BoxArray grid_cba = grid_ba;
         grid_cba.coarsen(ref_ratio);
-        const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba, ncell, do_pml_Lo, do_pml_Hi);
+        const BoxArray grid_cba_reduced = BoxArray(grid_cba.boxList().intersect(domain0));
+
+        const BoxArray& cba = (do_pml_in_domain) ?
+            MakeBoxArray(*cgeom, grid_cba_reduced, ncell, 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);
 
         DistributionMapping cdm{cba};
 
@@ -422,9 +511,20 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
         {
             pml_F_cp.reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()), cdm, 3, ngf));
             pml_F_cp->setVal(0.0);
-        }
 
-        sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta));
+        }
+        pml_j_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::jx_nodal_flag), cdm, 1, ngb));
+        pml_j_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::jy_nodal_flag), cdm, 1, ngb));
+        pml_j_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::jz_nodal_flag), cdm, 1, ngb));
+        pml_j_cp[0]->setVal(0.0);
+        pml_j_cp[1]->setVal(0.0);
+        pml_j_cp[2]->setVal(0.0);
+
+        if (do_pml_in_domain){
+            sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba_reduced, cgeom->CellSize(), ncell, delta));
+        } else {
+            sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta));
+        }
 
 #ifdef WARPX_USE_PSATD
         const bool in_pml = true; // Tells spectral solver to use split-PML equations
@@ -439,7 +539,8 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
 }
 
 BoxArray
-PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell,
+PML::MakeBoxArray (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();
@@ -453,14 +554,18 @@ PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba,
             }
         }
     }
-
     BoxList bl;
     for (int i = 0, N = grid_ba.size(); i < N; ++i)
     {
         const Box& grid_bx = grid_ba[i];
         const IntVect& grid_bx_sz = grid_bx.size();
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell,
-                                         "Consider using larger amr.blocking_factor");
+
+        if (do_pml_in_domain == 0) {
+            // Make sure that, in the case of several distinct refinement patches,
+            //  the PML cells surrounding these patches cannot overlap
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell,
+                            "Consider using larger amr.blocking_factor");
+        }
 
         Box bx = grid_bx;
         bx.grow(ncell);
@@ -530,6 +635,12 @@ PML::GetB_fp ()
 }
 
 std::array<MultiFab*,3>
+PML::Getj_fp ()
+{
+    return {pml_j_fp[0].get(), pml_j_fp[1].get(), pml_j_fp[2].get()};
+}
+
+std::array<MultiFab*,3>
 PML::GetE_cp ()
 {
     return {pml_E_cp[0].get(), pml_E_cp[1].get(), pml_E_cp[2].get()};
@@ -541,6 +652,12 @@ PML::GetB_cp ()
     return {pml_B_cp[0].get(), pml_B_cp[1].get(), pml_B_cp[2].get()};
 }
 
+std::array<MultiFab*,3>
+PML::Getj_cp ()
+{
+    return {pml_j_cp[0].get(), pml_j_cp[1].get(), pml_j_cp[2].get()};
+}
+
 MultiFab*
 PML::GetF_fp ()
 {
@@ -555,116 +672,181 @@ PML::GetF_cp ()
 
 void
 PML::ExchangeB (const std::array<amrex::MultiFab*,3>& B_fp,
-                const std::array<amrex::MultiFab*,3>& B_cp)
+                const std::array<amrex::MultiFab*,3>& B_cp,
+                int do_pml_in_domain)
 {
-  ExchangeB(PatchType::fine, B_fp);
-  ExchangeB(PatchType::coarse, B_cp);
+  ExchangeB(PatchType::fine, B_fp, do_pml_in_domain);
+  ExchangeB(PatchType::coarse, B_cp, do_pml_in_domain);
 }
 
 void
 PML::ExchangeB (PatchType patch_type,
-                const std::array<amrex::MultiFab*,3>& Bp)
+                const std::array<amrex::MultiFab*,3>& Bp,
+                int do_pml_in_domain)
 {
     if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0])
     {
-        Exchange(*pml_B_fp[0], *Bp[0], *m_geom);
-        Exchange(*pml_B_fp[1], *Bp[1], *m_geom);
-        Exchange(*pml_B_fp[2], *Bp[2], *m_geom);
+        Exchange(*pml_B_fp[0], *Bp[0], *m_geom, do_pml_in_domain);
+        Exchange(*pml_B_fp[1], *Bp[1], *m_geom, do_pml_in_domain);
+        Exchange(*pml_B_fp[2], *Bp[2], *m_geom, do_pml_in_domain);
     }
     else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0])
     {
-        Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom);
-        Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom);
-        Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom);
+        Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom, do_pml_in_domain);
+        Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom, do_pml_in_domain);
+        Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom, do_pml_in_domain);
     }
 }
 
 void
 PML::ExchangeE (const std::array<amrex::MultiFab*,3>& E_fp,
-                const std::array<amrex::MultiFab*,3>& E_cp)
+                const std::array<amrex::MultiFab*,3>& E_cp,
+                int do_pml_in_domain)
 {
-    ExchangeE(PatchType::fine, E_fp);
-    ExchangeE(PatchType::coarse, E_cp);
+    ExchangeE(PatchType::fine, E_fp, do_pml_in_domain);
+    ExchangeE(PatchType::coarse, E_cp, do_pml_in_domain);
 }
 
 void
 PML::ExchangeE (PatchType patch_type,
-                const std::array<amrex::MultiFab*,3>& Ep)
+                const std::array<amrex::MultiFab*,3>& Ep,
+                int do_pml_in_domain)
 {
     if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0])
     {
-        Exchange(*pml_E_fp[0], *Ep[0], *m_geom);
-        Exchange(*pml_E_fp[1], *Ep[1], *m_geom);
-        Exchange(*pml_E_fp[2], *Ep[2], *m_geom);
+        Exchange(*pml_E_fp[0], *Ep[0], *m_geom, do_pml_in_domain);
+        Exchange(*pml_E_fp[1], *Ep[1], *m_geom, do_pml_in_domain);
+        Exchange(*pml_E_fp[2], *Ep[2], *m_geom, do_pml_in_domain);
     }
     else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0])
     {
-        Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom);
-        Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom);
-        Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom);
+        Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom, do_pml_in_domain);
+        Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom, do_pml_in_domain);
+        Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom, do_pml_in_domain);
     }
 }
 
 void
-PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp)
+PML::CopyJtoPMLs (PatchType patch_type,
+                const std::array<amrex::MultiFab*,3>& jp)
+{
+    if (patch_type == PatchType::fine && pml_j_fp[0] && jp[0])
+    {
+        CopyToPML(*pml_j_fp[0], *jp[0], *m_geom);
+        CopyToPML(*pml_j_fp[1], *jp[1], *m_geom);
+        CopyToPML(*pml_j_fp[2], *jp[2], *m_geom);
+    }
+    else if (patch_type == PatchType::coarse && pml_j_cp[0] && jp[0])
+    {
+        CopyToPML(*pml_j_cp[0], *jp[0], *m_cgeom);
+        CopyToPML(*pml_j_cp[1], *jp[1], *m_cgeom);
+        CopyToPML(*pml_j_cp[2], *jp[2], *m_cgeom);
+    }
+}
+
+void
+PML::CopyJtoPMLs (const std::array<amrex::MultiFab*,3>& j_fp,
+                const std::array<amrex::MultiFab*,3>& j_cp)
+{
+    CopyJtoPMLs(PatchType::fine, j_fp);
+    CopyJtoPMLs(PatchType::coarse, j_cp);
+}
+
+
+void
+PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain)
 {
-    ExchangeF(PatchType::fine, F_fp);
-    ExchangeF(PatchType::coarse, F_cp);
+    ExchangeF(PatchType::fine, F_fp, do_pml_in_domain);
+    ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain);
 }
 
 void
-PML::ExchangeF (PatchType patch_type, MultiFab* Fp)
+PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain)
 {
     if (patch_type == PatchType::fine && pml_F_fp && Fp) {
-        Exchange(*pml_F_fp, *Fp, *m_geom);
+        Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain);
     } else if (patch_type == PatchType::coarse && pml_F_cp && Fp) {
-        Exchange(*pml_F_cp, *Fp, *m_cgeom);
+        Exchange(*pml_F_cp, *Fp, *m_cgeom, do_pml_in_domain);
     }
 }
 
+
 void
-PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom)
+PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom,
+                int do_pml_in_domain)
 {
+    BL_PROFILE("PML::Exchange");
+
     const IntVect& ngr = reg.nGrowVect();
     const IntVect& ngp = pml.nGrowVect();
     const int ncp = pml.nComp();
     const auto& period = geom.periodicity();
 
+    // Create temporary MultiFab to copy to and from the PML
     MultiFab tmpregmf(reg.boxArray(), reg.DistributionMap(), ncp, ngr);
 
-    if (ngp.max() > 0)  // Copy from pml to the ghost cells of regular data
-    {
-        MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, 0);
-        MultiFab::LinComb(totpmlmf, 1.0, pml, 0, 1.0, pml, 1, 0, 1, 0);
-        if (ncp == 3) {
-            MultiFab::Add(totpmlmf,pml,2,0,1,0);
-        }
-
-        MultiFab::Copy(tmpregmf, reg, 0, 0, 1, ngr);
-        tmpregmf.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), ngr, period);
+    // Create the sum of the split fields, in the PML
+    MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, 0); // Allocate
+    MultiFab::LinComb(totpmlmf, 1.0, pml, 0, 1.0, pml, 1, 0, 1, 0); // Sum
+    if (ncp == 3) {
+        MultiFab::Add(totpmlmf,pml,2,0,1,0); // Sum the third split component
+    }
 
+    // Copy from the sum of PML split field to valid cells of regular grid
+    if (do_pml_in_domain){
+        // Valid cells of the PML and of the regular grid overlap
+        // Copy from valid cells of the PML to valid cells of the regular grid
+        reg.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), IntVect(0), period);
+    } else {
+        // Valid cells of the PML only overlap with guard cells of regular grid
+        // (and outermost valid cell of the regular grid, for nodal direction)
+        // Copy from valid cells of PML to ghost cells of regular grid
+        // but avoid updating the outermost valid cell
+        if (ngr.max() > 0) {
+            MultiFab::Copy(tmpregmf, reg, 0, 0, 1, ngr);
+            tmpregmf.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), ngr, period);
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
-        for (MFIter mfi(reg); mfi.isValid(); ++mfi)
-        {
-            const FArrayBox& src = tmpregmf[mfi];
-            FArrayBox& dst = reg[mfi];
-            const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox());
-            for (const Box& bx : bl)
+            for (MFIter mfi(reg); mfi.isValid(); ++mfi)
             {
-                dst.copy(src, bx, 0, bx, 0, 1);
+                const FArrayBox& src = tmpregmf[mfi];
+                FArrayBox& dst = reg[mfi];
+                const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox());
+                // boxDiff avoids the outermost valid cell
+                for (const Box& bx : bl) {
+                    dst.copy(src, bx, 0, bx, 0, 1);
+                }
             }
         }
     }
 
-    // Copy from regular data to PML's first component
+    // Copy from valid cells of the regular grid to guard cells of the PML
+    // (and outermost valid cell in the nodal direction)
+    // More specifically, copy from regular data to PML's first component
     // Zero out the second (and third) component
-    MultiFab::Copy(tmpregmf,reg,0,0,1,0);
-    tmpregmf.setVal(0.0, 1, ncp-1, 0);
+    MultiFab::Copy(tmpregmf,reg,0,0,1,0); // Fill first component of tmpregmf
+    tmpregmf.setVal(0.0, 1, ncp-1, 0); // Zero out the second (and third) component
+    if (do_pml_in_domain){
+        // Where valid cells of tmpregmf overlap with PML valid cells,
+        // copy the PML (this is order to avoid overwriting PML valid cells,
+        // in the next `ParallelCopy`)
+        tmpregmf.ParallelCopy(pml,0, 0, ncp, IntVect(0), IntVect(0), period);
+    }
     pml.ParallelCopy(tmpregmf, 0, 0, ncp, IntVect(0), ngp, period);
 }
 
+
+void
+PML::CopyToPML (MultiFab& pml, MultiFab& reg, const Geometry& geom)
+{
+  const IntVect& ngr = reg.nGrowVect();
+  const IntVect& ngp = pml.nGrowVect();
+  const auto& period = geom.periodicity();
+
+  pml.ParallelCopy(reg, 0, 0, 1, ngr, ngp, period);
+}
+
 void
 PML::FillBoundary ()
 {