fixing merge conflicts.

1 files changed, 296 insertions, 292 deletions

diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp
index ecbe3c117..a05ec26c1 100644
--- a/Source/WarpXEvolve.cpp
+++ b/Source/WarpXEvolve.cpp
@@ -93,7 +93,9 @@ WarpX::EvolveEM (int numsteps)
         } else {
            // Beyond one step, we have E^{n} and B^{n}.
            // Particles have p^{n-1/2} and x^{n}.
-           UpdateAuxilaryData();
+	  FillBoundaryE();
+	  FillBoundaryB();
+	  UpdateAuxilaryData();
        }
 
         // Push particle from x^{n} to x^{n+1}
@@ -262,104 +264,103 @@ void
 WarpX::EvolveB (int lev, Real dt)
 {
     BL_PROFILE("WarpX::EvolveB()");
-
-    int npatches = (lev == 0) ? 1 : 2;
-
-    for (int ipatch = 0; ipatch < npatches; ++ipatch)
+    EvolveB(lev, PatchType::fine, dt);
+    if (lev > 0)
     {
-        int patch_level = (ipatch == 0) ? lev : lev-1;
-        const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
-        const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
+        EvolveB(lev, PatchType::coarse, dt);
+    }
+}
 
-        MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz;
-        if (ipatch == 0)
-        {
-            Ex = Efield_fp[lev][0].get();
-            Ey = Efield_fp[lev][1].get();
-            Ez = Efield_fp[lev][2].get();
-            Bx = Bfield_fp[lev][0].get();
-            By = Bfield_fp[lev][1].get();
-            Bz = Bfield_fp[lev][2].get();
-        }
-        else
-        {
-            Ex = Efield_cp[lev][0].get();
-            Ey = Efield_cp[lev][1].get();
-            Ez = Efield_cp[lev][2].get();
-            Bx = Bfield_cp[lev][0].get();
-            By = Bfield_cp[lev][1].get();
-            Bz = Bfield_cp[lev][2].get();
-        }
+void
+WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt)
+{
+    const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
+    const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
+    const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
+    
+    MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz;
+    if (patch_type == PatchType::fine)
+    {
+        Ex = Efield_fp[lev][0].get();
+        Ey = Efield_fp[lev][1].get();
+        Ez = Efield_fp[lev][2].get();
+        Bx = Bfield_fp[lev][0].get();
+        By = Bfield_fp[lev][1].get();
+        Bz = Bfield_fp[lev][2].get();
+    }
+    else
+    {
+        Ex = Efield_cp[lev][0].get();
+        Ey = Efield_cp[lev][1].get();
+        Ez = Efield_cp[lev][2].get();
+        Bx = Bfield_cp[lev][0].get();
+        By = Bfield_cp[lev][1].get();
+        Bz = Bfield_cp[lev][2].get();
+    }
 
-        MultiFab* cost = costs[lev].get();
-        const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector();
+    MultiFab* cost = costs[lev].get();
+    const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector();
 
-        // Loop through the grids, and over the tiles within each grid
+    // Loop through the grids, and over the tiles within each grid
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
-        for ( MFIter mfi(*Bx,true); mfi.isValid(); ++mfi )
-        {
-            Real wt = amrex::second();
-
-            const Box& tbx  = mfi.tilebox(Bx_nodal_flag);
-            const Box& tby  = mfi.tilebox(By_nodal_flag);
-            const Box& tbz  = mfi.tilebox(Bz_nodal_flag);
-
-            // Call picsar routine for each tile
-            warpx_push_bvec(
-                tbx.loVect(), tbx.hiVect(),
-                tby.loVect(), tby.hiVect(),
-                tbz.loVect(), tbz.hiVect(),
-                BL_TO_FORTRAN_3D((*Ex)[mfi]),
-                BL_TO_FORTRAN_3D((*Ey)[mfi]),
-                BL_TO_FORTRAN_3D((*Ez)[mfi]),
-                BL_TO_FORTRAN_3D((*Bx)[mfi]),
-                BL_TO_FORTRAN_3D((*By)[mfi]),
-                BL_TO_FORTRAN_3D((*Bz)[mfi]),
-	        &dtsdx[0], &dtsdx[1], &dtsdx[2],
-	        &WarpX::maxwell_fdtd_solver_id);
-
-            if (cost) {
-                Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)});
-                if (ipatch == 1) cbx.refine(rr);
-                wt = (amrex::second() - wt) / cbx.d_numPts();
-                (*cost)[mfi].plus(wt, cbx);
-            }
+    for ( MFIter mfi(*Bx,true); mfi.isValid(); ++mfi )
+    {
+        Real wt = amrex::second();
+        
+        const Box& tbx  = mfi.tilebox(Bx_nodal_flag);
+        const Box& tby  = mfi.tilebox(By_nodal_flag);
+        const Box& tbz  = mfi.tilebox(Bz_nodal_flag);
+        
+        // Call picsar routine for each tile
+        warpx_push_bvec(
+		      tbx.loVect(), tbx.hiVect(),
+		      tby.loVect(), tby.hiVect(),
+		      tbz.loVect(), tbz.hiVect(),
+		      BL_TO_FORTRAN_3D((*Ex)[mfi]),
+		      BL_TO_FORTRAN_3D((*Ey)[mfi]),
+		      BL_TO_FORTRAN_3D((*Ez)[mfi]),
+		      BL_TO_FORTRAN_3D((*Bx)[mfi]),
+		      BL_TO_FORTRAN_3D((*By)[mfi]),
+		      BL_TO_FORTRAN_3D((*Bz)[mfi]),
+		      &dtsdx[0], &dtsdx[1], &dtsdx[2],
+		      &WarpX::maxwell_fdtd_solver_id);
+        
+        if (cost) {
+            Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)});
+            if (patch_type == PatchType::coarse) cbx.refine(rr);
+            wt = (amrex::second() - wt) / cbx.d_numPts();
+            (*cost)[mfi].plus(wt, cbx);
         }
     }
 
     if (do_pml && pml[lev]->ok())
     {
-        for (int ipatch = 0; ipatch < npatches; ++ipatch)
-        {
-            const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
-            const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
-	    int patch_level = (ipatch == 0) ? lev : lev-1;
-	    const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
-	    const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
+        const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
+        const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
+
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
-            for ( MFIter mfi(*pml_B[0],true); mfi.isValid(); ++mfi )
-            {
-                const Box& tbx  = mfi.tilebox(Bx_nodal_flag);
-                const Box& tby  = mfi.tilebox(By_nodal_flag);
-                const Box& tbz  = mfi.tilebox(Bz_nodal_flag);
-
-                WRPX_PUSH_PML_BVEC(
-                    tbx.loVect(), tbx.hiVect(),
-                    tby.loVect(), tby.hiVect(),
-                    tbz.loVect(), tbz.hiVect(),
-                    BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
-                    BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
-                    BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
-                    BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
-                    BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
-                    BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
-		    &dtsdx[0], &dtsdx[1], &dtsdx[2],
-		    &WarpX::maxwell_fdtd_solver_id);
-            }
+        for ( MFIter mfi(*pml_B[0],true); mfi.isValid(); ++mfi )
+        {
+            const Box& tbx  = mfi.tilebox(Bx_nodal_flag);
+            const Box& tby  = mfi.tilebox(By_nodal_flag);
+            const Box& tbz  = mfi.tilebox(Bz_nodal_flag);
+
+            WRPX_PUSH_PML_BVEC(
+			     tbx.loVect(), tbx.hiVect(),
+			     tby.loVect(), tby.hiVect(),
+			     tbz.loVect(), tbz.hiVect(),
+			     BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
+			     BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
+			     BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
+			     BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
+			     BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
+			     BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
+			     &dtsdx[0], &dtsdx[1], &dtsdx[2],
+			     &WarpX::maxwell_fdtd_solver_id);
         }
     }
 }
@@ -367,7 +368,8 @@ WarpX::EvolveB (int lev, Real dt)
 void
 WarpX::EvolveE (Real dt)
 {
-    for (int lev = 0; lev <= finest_level; ++lev) {
+    for (int lev = 0; lev <= finest_level; ++lev)
+    {
         EvolveE(lev, dt);
     }
 }
@@ -376,148 +378,145 @@ void
 WarpX::EvolveE (int lev, Real dt)
 {
     BL_PROFILE("WarpX::EvolveE()");
+    EvolveE(lev, PatchType::fine, dt);
+    if (lev > 0)
+    {
+        EvolveE(lev, PatchType::coarse, dt);
+    }
+}
 
+void
+WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt)
+{
     const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt;
     const Real c2dt = (PhysConst::c*PhysConst::c) * dt;
 
-    int npatches = (lev == 0) ? 1 : 2;
+    int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
+    const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
+    const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]};
 
-    for (int ipatch = 0; ipatch < npatches; ++ipatch)
+    MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F;
+    if (patch_type == PatchType::fine)
     {
-        int patch_level = (ipatch == 0) ? lev : lev-1;
-        const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
-        const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]};
-
-        MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F;
-        if (ipatch == 0)
-        {
-            Ex = Efield_fp[lev][0].get();
-            Ey = Efield_fp[lev][1].get();
-            Ez = Efield_fp[lev][2].get();
-            Bx = Bfield_fp[lev][0].get();
-            By = Bfield_fp[lev][1].get();
-            Bz = Bfield_fp[lev][2].get();
-            jx = current_fp[lev][0].get();
-            jy = current_fp[lev][1].get();
-            jz = current_fp[lev][2].get();
-            F  = F_fp[lev].get();
-        }
-        else
-        {
-            Ex = Efield_cp[lev][0].get();
-            Ey = Efield_cp[lev][1].get();
-            Ez = Efield_cp[lev][2].get();
-            Bx = Bfield_cp[lev][0].get();
-            By = Bfield_cp[lev][1].get();
-            Bz = Bfield_cp[lev][2].get();
-            jx = current_cp[lev][0].get();
-            jy = current_cp[lev][1].get();
-            jz = current_cp[lev][2].get();
-            F  = F_cp[lev].get();
-        }
+        Ex = Efield_fp[lev][0].get();
+        Ey = Efield_fp[lev][1].get();
+        Ez = Efield_fp[lev][2].get();
+        Bx = Bfield_fp[lev][0].get();
+        By = Bfield_fp[lev][1].get();
+        Bz = Bfield_fp[lev][2].get();
+        jx = current_fp[lev][0].get();
+        jy = current_fp[lev][1].get();
+        jz = current_fp[lev][2].get();
+        F  = F_fp[lev].get();
+    }
+    else if (patch_type == PatchType::coarse)
+    {
+        Ex = Efield_cp[lev][0].get();
+        Ey = Efield_cp[lev][1].get();
+        Ez = Efield_cp[lev][2].get();
+        Bx = Bfield_cp[lev][0].get();
+        By = Bfield_cp[lev][1].get();
+        Bz = Bfield_cp[lev][2].get();
+        jx = current_cp[lev][0].get();
+        jy = current_cp[lev][1].get();
+        jz = current_cp[lev][2].get();
+        F  = F_cp[lev].get();
+    }
 
-        MultiFab* cost = costs[lev].get();
-        const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector();
+    MultiFab* cost = costs[lev].get();
+    const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector();
 
-        // Loop through the grids, and over the tiles within each grid
+    // Loop through the grids, and over the tiles within each grid
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
-        for ( MFIter mfi(*Ex,true); mfi.isValid(); ++mfi )
+    for ( MFIter mfi(*Ex,true); mfi.isValid(); ++mfi )
+    {
+        Real wt = amrex::second();
+
+        const Box& tex  = mfi.tilebox(Ex_nodal_flag);
+        const Box& tey  = mfi.tilebox(Ey_nodal_flag);
+        const Box& tez  = mfi.tilebox(Ez_nodal_flag);
+
+        // Call picsar routine for each tile
+        warpx_push_evec(
+		      tex.loVect(), tex.hiVect(),
+		      tey.loVect(), tey.hiVect(),
+		      tez.loVect(), tez.hiVect(),
+		      BL_TO_FORTRAN_3D((*Ex)[mfi]),
+		      BL_TO_FORTRAN_3D((*Ey)[mfi]),
+		      BL_TO_FORTRAN_3D((*Ez)[mfi]),
+		      BL_TO_FORTRAN_3D((*Bx)[mfi]),
+		      BL_TO_FORTRAN_3D((*By)[mfi]),
+		      BL_TO_FORTRAN_3D((*Bz)[mfi]),
+		      BL_TO_FORTRAN_3D((*jx)[mfi]),
+		      BL_TO_FORTRAN_3D((*jy)[mfi]),
+		      BL_TO_FORTRAN_3D((*jz)[mfi]),
+		      &mu_c2_dt,
+		      &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]);
+
+        if (F)
         {
-            Real wt = amrex::second();
-
-            const Box& tex  = mfi.tilebox(Ex_nodal_flag);
-            const Box& tey  = mfi.tilebox(Ey_nodal_flag);
-            const Box& tez  = mfi.tilebox(Ez_nodal_flag);
-
-            // Call picsar routine for each tile
-            warpx_push_evec(
-                tex.loVect(), tex.hiVect(),
-                tey.loVect(), tey.hiVect(),
-                tez.loVect(), tez.hiVect(),
-                BL_TO_FORTRAN_3D((*Ex)[mfi]),
-                BL_TO_FORTRAN_3D((*Ey)[mfi]),
-                BL_TO_FORTRAN_3D((*Ez)[mfi]),
-                BL_TO_FORTRAN_3D((*Bx)[mfi]),
-                BL_TO_FORTRAN_3D((*By)[mfi]),
-                BL_TO_FORTRAN_3D((*Bz)[mfi]),
-                BL_TO_FORTRAN_3D((*jx)[mfi]),
-                BL_TO_FORTRAN_3D((*jy)[mfi]),
-                BL_TO_FORTRAN_3D((*jz)[mfi]),
-                &mu_c2_dt,
-                &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]);
-
-            if (F) {
-
-                // Call picsar routine for each tile
-                warpx_push_evec_f(
-                    tex.loVect(), tex.hiVect(),
-                    tey.loVect(), tey.hiVect(),
-                    tez.loVect(), tez.hiVect(),
-                    BL_TO_FORTRAN_3D((*Ex)[mfi]),
-                    BL_TO_FORTRAN_3D((*Ey)[mfi]),
-                    BL_TO_FORTRAN_3D((*Ez)[mfi]),
-                    BL_TO_FORTRAN_3D((*F)[mfi]),
-                    &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2],
-                    &WarpX::maxwell_fdtd_solver_id);
-
-            }
+            warpx_push_evec_f(
+			  tex.loVect(), tex.hiVect(),
+			  tey.loVect(), tey.hiVect(),
+			  tez.loVect(), tez.hiVect(),
+			  BL_TO_FORTRAN_3D((*Ex)[mfi]),
+			  BL_TO_FORTRAN_3D((*Ey)[mfi]),
+			  BL_TO_FORTRAN_3D((*Ez)[mfi]),
+			  BL_TO_FORTRAN_3D((*F)[mfi]),
+			  &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2],
+			  &WarpX::maxwell_fdtd_solver_id);
+        }
 
-            if (cost) {
-                Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)});
-                if (ipatch == 1) cbx.refine(rr);
-                wt = (amrex::second() - wt) / cbx.d_numPts();
-                (*cost)[mfi].plus(wt, cbx);
-            }
+        if (cost) {
+            Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)});
+            if (patch_type == PatchType::coarse) cbx.refine(rr);
+            wt = (amrex::second() - wt) / cbx.d_numPts();
+            (*cost)[mfi].plus(wt, cbx);
         }
     }
 
     if (do_pml && pml[lev]->ok())
     {
+        if (F) pml[lev]->ExchangeF(patch_type, F);
 
-        for (int ipatch = 0; ipatch < npatches; ++ipatch)
-        {
-            const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
-            const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
-            const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
-            int patch_level = (ipatch == 0) ? lev : lev-1;
-            const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
-            const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]};
+        const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
+        const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
+        const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
-            for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi )
+        for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi )
+        {
+            const Box& tex  = mfi.tilebox(Ex_nodal_flag);
+            const Box& tey  = mfi.tilebox(Ey_nodal_flag);
+            const Box& tez  = mfi.tilebox(Ez_nodal_flag);
+
+            WRPX_PUSH_PML_EVEC(
+			     tex.loVect(), tex.hiVect(),
+			     tey.loVect(), tey.hiVect(),
+			     tez.loVect(), tez.hiVect(),
+			     BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
+			     BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
+			     BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
+			     BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
+			     BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
+			     BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
+			     &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]);
+            
+            if (pml_F)
             {
-                const Box& tex  = mfi.tilebox(Ex_nodal_flag);
-                const Box& tey  = mfi.tilebox(Ey_nodal_flag);
-                const Box& tez  = mfi.tilebox(Ez_nodal_flag);
-
-                WRPX_PUSH_PML_EVEC(
-                    tex.loVect(), tex.hiVect(),
-                    tey.loVect(), tey.hiVect(),
-                    tez.loVect(), tez.hiVect(),
-                    BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
-                    BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
-                    BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
-                    BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
-                    BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
-                    BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
-                    &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]);
-
-                if (pml_F)
-                {
-                    WRPX_PUSH_PML_EVEC_F(
-                        tex.loVect(), tex.hiVect(),
-                        tey.loVect(), tey.hiVect(),
-                        tez.loVect(), tez.hiVect(),
-                        BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
-                        BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
-                        BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
-                        BL_TO_FORTRAN_3D((*pml_F   )[mfi]),
-                        &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2],
-                        &WarpX::maxwell_fdtd_solver_id);
-                }
+                WRPX_PUSH_PML_EVEC_F(
+				   tex.loVect(), tex.hiVect(),
+				   tey.loVect(), tey.hiVect(),
+				   tez.loVect(), tez.hiVect(),
+				   BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
+				   BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
+				   BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
+				   BL_TO_FORTRAN_3D((*pml_F   )[mfi]),
+				   &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2],
+				   &WarpX::maxwell_fdtd_solver_id);
             }
         }
     }
@@ -528,7 +527,8 @@ WarpX::EvolveF (Real dt, DtType dt_type)
 {
     if (!do_dive_cleaning) return;
 
-    for (int lev = 0; lev <= finest_level; ++lev) {
+    for (int lev = 0; lev <= finest_level; ++lev)
+    {
         EvolveF(lev, dt, dt_type);
     }
 }
@@ -538,61 +538,65 @@ WarpX::EvolveF (int lev, Real dt, DtType dt_type)
 {
     if (!do_dive_cleaning) return;
 
+    EvolveF(lev, PatchType::fine, dt, dt_type);
+    if (lev > 0) EvolveF(lev, PatchType::coarse, dt, dt_type);
+}
+
+void
+WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type)
+{
+    if (!do_dive_cleaning) return;
+
     BL_PROFILE("WarpX::EvolveF()");
 
     static constexpr Real mu_c2 = PhysConst::mu0*PhysConst::c*PhysConst::c;
 
-    int npatches = (lev == 0) ? 1 : 2;
+    int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
+    const auto& dx = WarpX::CellSize(patch_level);
+    const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
 
-    for (int ipatch = 0; ipatch < npatches; ++ipatch)
+    MultiFab *Ex, *Ey, *Ez, *rho, *F;
+    if (patch_type == PatchType::fine)
     {
-        int patch_level = (ipatch == 0) ? lev : lev-1;
-        const auto& dx = WarpX::CellSize(patch_level);
-        const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
-
-        MultiFab *Ex, *Ey, *Ez, *rho, *F;
-        if (ipatch == 0)
-        {
-            Ex = Efield_fp[lev][0].get();
-            Ey = Efield_fp[lev][1].get();
-            Ez = Efield_fp[lev][2].get();
-            rho = rho_fp[lev].get();
-            F = F_fp[lev].get();
-        }
-        else
-        {
-            Ex = Efield_cp[lev][0].get();
-            Ey = Efield_cp[lev][1].get();
-            Ez = Efield_cp[lev][2].get();
-            rho = rho_cp[lev].get();
-            F = F_cp[lev].get();
-        }
-
-        const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1;
-
-        MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0);
-        ComputeDivE(src, 0, {Ex,Ey,Ez}, dx);
-        MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0);
-        MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0);
-
-        if (do_pml && pml[lev]->ok())
-        {
-            const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
-            const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
+        Ex = Efield_fp[lev][0].get();
+        Ey = Efield_fp[lev][1].get();
+        Ez = Efield_fp[lev][2].get();
+        rho = rho_fp[lev].get();
+        F = F_fp[lev].get();
+    }
+    else
+    {
+        Ex = Efield_cp[lev][0].get();
+        Ey = Efield_cp[lev][1].get();
+        Ez = Efield_cp[lev][2].get();
+        rho = rho_cp[lev].get();
+        F = F_cp[lev].get();
+    }
+    
+    const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1;
+
+    MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0);
+    ComputeDivE(src, 0, {Ex,Ey,Ez}, dx);
+    MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0);
+    MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0);
+    
+    if (do_pml && pml[lev]->ok())
+    {
+        const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
+        const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
 
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
-            for ( MFIter mfi(*pml_F,true); mfi.isValid(); ++mfi )
-            {
-                const Box& bx = mfi.tilebox();
-                WRPX_PUSH_PML_F(bx.loVect(), bx.hiVect(),
-                                BL_TO_FORTRAN_ANYD((*pml_F   )[mfi]),
-                                BL_TO_FORTRAN_ANYD((*pml_E[0])[mfi]),
-                                BL_TO_FORTRAN_ANYD((*pml_E[1])[mfi]),
-                                BL_TO_FORTRAN_ANYD((*pml_E[2])[mfi]),
-                                &dtsdx[0], &dtsdx[1], &dtsdx[2]);
-            }
+        for ( MFIter mfi(*pml_F,true); mfi.isValid(); ++mfi )
+        {
+            const Box& bx = mfi.tilebox();
+            WRPX_PUSH_PML_F(bx.loVect(), bx.hiVect(),
+			  BL_TO_FORTRAN_ANYD((*pml_F   )[mfi]),
+			  BL_TO_FORTRAN_ANYD((*pml_E[0])[mfi]),
+			  BL_TO_FORTRAN_ANYD((*pml_E[1])[mfi]),
+			  BL_TO_FORTRAN_ANYD((*pml_E[2])[mfi]),
+			  &dtsdx[0], &dtsdx[1], &dtsdx[2]);
         }
     }
 }
@@ -600,8 +604,6 @@ WarpX::EvolveF (int lev, Real dt, DtType dt_type)
 void
 WarpX::DampPML ()
 {
-    if (!do_pml) return;
-
     for (int lev = 0; lev <= finest_level; ++lev) {
         DampPML(lev);
     }
@@ -610,54 +612,56 @@ WarpX::DampPML ()
 void
 WarpX::DampPML (int lev)
 {
+    DampPML(lev, PatchType::fine);
+    if (lev > 0) DampPML(lev, PatchType::coarse);
+}
+
+void
+WarpX::DampPML (int lev, PatchType patch_type)
+{
     if (!do_pml) return;
 
     BL_PROFILE("WarpX::DampPML()");
 
-    int npatches = (lev == 0) ? 1 : 2;
-
-    for (int ipatch = 0; ipatch < npatches; ++ipatch)
+    if (pml[lev]->ok())
     {
-        if (pml[lev]->ok())
-        {
-            const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
-            const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
-            const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
-            const auto& sigba = (ipatch==0) ? pml[lev]->GetMultiSigmaBox_fp()
-                                            : pml[lev]->GetMultiSigmaBox_cp();
+        const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
+        const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
+        const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
+        const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp()
+                                                              : pml[lev]->GetMultiSigmaBox_cp();
 
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
-            for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi )
-            {
-                const Box& tex  = mfi.tilebox(Ex_nodal_flag);
-                const Box& tey  = mfi.tilebox(Ey_nodal_flag);
-                const Box& tez  = mfi.tilebox(Ez_nodal_flag);
-                const Box& tbx  = mfi.tilebox(Bx_nodal_flag);
-                const Box& tby  = mfi.tilebox(By_nodal_flag);
-                const Box& tbz  = mfi.tilebox(Bz_nodal_flag);
-
-                WRPX_DAMP_PML(tex.loVect(), tex.hiVect(),
-                              tey.loVect(), tey.hiVect(),
-                              tez.loVect(), tez.hiVect(),
-                              tbx.loVect(), tbx.hiVect(),
-                              tby.loVect(), tby.hiVect(),
-                              tbz.loVect(), tbz.hiVect(),
-                              BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
-                              BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
-                              BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
-                              BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
-                              BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
-                              BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
-                              WRPX_PML_TO_FORTRAN(sigba[mfi]));
-
-                if (pml_F) {
-                    const Box& tnd  = mfi.nodaltilebox();
-                    WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(),
-                                    BL_TO_FORTRAN_3D((*pml_F)[mfi]),
-                                    WRPX_PML_TO_FORTRAN(sigba[mfi]));
-                }
+        for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi )
+        {
+            const Box& tex  = mfi.tilebox(Ex_nodal_flag);
+            const Box& tey  = mfi.tilebox(Ey_nodal_flag);
+            const Box& tez  = mfi.tilebox(Ez_nodal_flag);
+            const Box& tbx  = mfi.tilebox(Bx_nodal_flag);
+            const Box& tby  = mfi.tilebox(By_nodal_flag);
+            const Box& tbz  = mfi.tilebox(Bz_nodal_flag);
+            
+            WRPX_DAMP_PML(tex.loVect(), tex.hiVect(),
+			    tey.loVect(), tey.hiVect(),
+			    tez.loVect(), tez.hiVect(),
+			    tbx.loVect(), tbx.hiVect(),
+    			tby.loVect(), tby.hiVect(),
+	    		tbz.loVect(), tbz.hiVect(),
+		    	BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
+			    BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
+			    BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
+			    BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
+			    BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
+			    BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
+			    WRPX_PML_TO_FORTRAN(sigba[mfi]));
+            
+            if (pml_F) {
+                const Box& tnd  = mfi.nodaltilebox();
+                WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(),
+			        BL_TO_FORTRAN_3D((*pml_F)[mfi]),
+			        WRPX_PML_TO_FORTRAN(sigba[mfi]));
             }
         }
     }
@@ -681,7 +685,7 @@ WarpX::PushParticlesandDepose (int lev, Real cur_time)
                  *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2],
                  *current_fp[lev][0],*current_fp[lev][1],*current_fp[lev][2],
                  current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(),
-                 rho_fp[lev].get(),
+                 rho_fp[lev].get(), charge_buf[lev].get(),
                  Efield_cax[lev][0].get(), Efield_cax[lev][1].get(), Efield_cax[lev][2].get(),
                  Bfield_cax[lev][0].get(), Bfield_cax[lev][1].get(), Bfield_cax[lev][2].get(),
                  cur_time, dt[lev]);