Merge branch 'qed_bw_qs_factory_class' into qed_evolve_optical_depth

10 files changed, 709 insertions, 41 deletions

diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index 92f0b4f09..52df3dc25 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -54,6 +54,157 @@ WarpX::UpdateAuxilaryData ()
 {
     BL_PROFILE("UpdateAuxilaryData()");
 
+    if (Bfield_aux[0][0]->ixType() == Bfield_fp[0][0]->ixType()) {
+        UpdateAuxilaryDataSameType();
+    } else {
+        UpdateAuxilaryDataStagToNodal();
+    }
+}
+
+void
+WarpX::UpdateAuxilaryDataStagToNodal ()
+{
+    // For level 0, we only need to do the average.
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+    for (MFIter mfi(*Bfield_aux[0][0]); mfi.isValid(); ++mfi)
+    {
+        Array4<Real> const& bx_aux = Bfield_aux[0][0]->array(mfi);
+        Array4<Real> const& by_aux = Bfield_aux[0][1]->array(mfi);
+        Array4<Real> const& bz_aux = Bfield_aux[0][2]->array(mfi);
+        Array4<Real const> const& bx_fp = Bfield_fp[0][0]->const_array(mfi);
+        Array4<Real const> const& by_fp = Bfield_fp[0][1]->const_array(mfi);
+        Array4<Real const> const& bz_fp = Bfield_fp[0][2]->const_array(mfi);
+
+        Array4<Real> const& ex_aux = Efield_aux[0][0]->array(mfi);
+        Array4<Real> const& ey_aux = Efield_aux[0][1]->array(mfi);
+        Array4<Real> const& ez_aux = Efield_aux[0][2]->array(mfi);
+        Array4<Real const> const& ex_fp = Efield_fp[0][0]->const_array(mfi);
+        Array4<Real const> const& ey_fp = Efield_fp[0][1]->const_array(mfi);
+        Array4<Real const> const& ez_fp = Efield_fp[0][2]->const_array(mfi);
+
+        const Box& bx = mfi.fabbox();
+        amrex::ParallelFor(bx,
+        [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
+        {
+            warpx_interp_nd_bfield_x(j,k,l, bx_aux, bx_fp);
+            warpx_interp_nd_bfield_y(j,k,l, by_aux, by_fp);
+            warpx_interp_nd_bfield_z(j,k,l, bz_aux, bz_fp);
+            warpx_interp_nd_efield_x(j,k,l, ex_aux, ex_fp);
+            warpx_interp_nd_efield_y(j,k,l, ey_aux, ey_fp);
+            warpx_interp_nd_efield_z(j,k,l, ez_aux, ez_fp);
+        });
+    }
+
+    for (int lev = 1; lev <= finest_level; ++lev)
+    {
+        BoxArray const& nba = Bfield_aux[lev][0]->boxArray();
+        BoxArray const& cnba = amrex::coarsen(nba,2);
+        DistributionMapping const& dm = Bfield_aux[lev][0]->DistributionMap();
+        auto const& cperiod = Geom(lev-1).periodicity();
+
+        // Bfield
+        {
+            Array<std::unique_ptr<MultiFab>,3> Btmp;
+            if (Bfield_cax[lev][0]) {
+                for (int i = 0; i < 3; ++i) {
+                    Btmp[i].reset(new MultiFab(*Bfield_cax[lev][i], amrex::make_alias, 0, 1));
+                }
+            } else {
+                IntVect ngtmp = Bfield_aux[lev-1][0]->nGrowVect();
+                for (int i = 0; i < 3; ++i) {
+                    Btmp[i].reset(new MultiFab(cnba, dm, 1, ngtmp));
+                }
+            }
+            // ParallelCopy from coarse level
+            for (int i = 0; i < 3; ++i) {
+                IntVect ng = Btmp[i]->nGrowVect();
+                Btmp[i]->ParallelCopy(*Bfield_aux[lev-1][i], 0, 0, 1, ng, ng, cperiod);
+            }
+
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+            for (MFIter mfi(*Bfield_aux[lev][0]); mfi.isValid(); ++mfi)
+            {
+                Array4<Real> const& bx_aux = Bfield_aux[lev][0]->array(mfi);
+                Array4<Real> const& by_aux = Bfield_aux[lev][1]->array(mfi);
+                Array4<Real> const& bz_aux = Bfield_aux[lev][2]->array(mfi);
+                Array4<Real const> const& bx_fp = Bfield_fp[lev][0]->const_array(mfi);
+                Array4<Real const> const& by_fp = Bfield_fp[lev][1]->const_array(mfi);
+                Array4<Real const> const& bz_fp = Bfield_fp[lev][2]->const_array(mfi);
+                Array4<Real const> const& bx_cp = Bfield_cp[lev][0]->const_array(mfi);
+                Array4<Real const> const& by_cp = Bfield_cp[lev][1]->const_array(mfi);
+                Array4<Real const> const& bz_cp = Bfield_cp[lev][2]->const_array(mfi);
+                Array4<Real const> const& bx_c = Btmp[0]->const_array(mfi);
+                Array4<Real const> const& by_c = Btmp[1]->const_array(mfi);
+                Array4<Real const> const& bz_c = Btmp[2]->const_array(mfi);
+
+                const Box& bx = mfi.fabbox();
+                amrex::ParallelFor(bx,
+                [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
+                {
+                    warpx_interp_nd_bfield_x(j,k,l, bx_aux, bx_fp, bx_cp, bx_c);
+                    warpx_interp_nd_bfield_y(j,k,l, by_aux, by_fp, by_cp, by_c);
+                    warpx_interp_nd_bfield_z(j,k,l, bz_aux, bz_fp, bz_cp, bz_c);
+                });
+            }
+        }
+
+        // Efield
+        {
+            Array<std::unique_ptr<MultiFab>,3> Etmp;
+            if (Efield_cax[lev][0]) {
+                for (int i = 0; i < 3; ++i) {
+                    Etmp[i].reset(new MultiFab(*Efield_cax[lev][i], amrex::make_alias, 0, 1));
+                }
+            } else {
+                IntVect ngtmp = Efield_aux[lev-1][0]->nGrowVect();
+                for (int i = 0; i < 3; ++i) {
+                    Etmp[i].reset(new MultiFab(cnba, dm, 1, ngtmp));
+                }
+            }
+            // ParallelCopy from coarse level
+            for (int i = 0; i < 3; ++i) {
+                IntVect ng = Etmp[i]->nGrowVect();
+                Etmp[i]->ParallelCopy(*Efield_aux[lev-1][i], 0, 0, 1, ng, ng, cperiod);
+            }
+
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+            for (MFIter mfi(*Efield_aux[lev][0]); mfi.isValid(); ++mfi)
+            {
+                Array4<Real> const& ex_aux = Efield_aux[lev][0]->array(mfi);
+                Array4<Real> const& ey_aux = Efield_aux[lev][1]->array(mfi);
+                Array4<Real> const& ez_aux = Efield_aux[lev][2]->array(mfi);
+                Array4<Real const> const& ex_fp = Efield_fp[lev][0]->const_array(mfi);
+                Array4<Real const> const& ey_fp = Efield_fp[lev][1]->const_array(mfi);
+                Array4<Real const> const& ez_fp = Efield_fp[lev][2]->const_array(mfi);
+                Array4<Real const> const& ex_cp = Efield_cp[lev][0]->const_array(mfi);
+                Array4<Real const> const& ey_cp = Efield_cp[lev][1]->const_array(mfi);
+                Array4<Real const> const& ez_cp = Efield_cp[lev][2]->const_array(mfi);
+                Array4<Real const> const& ex_c = Etmp[0]->const_array(mfi);
+                Array4<Real const> const& ey_c = Etmp[1]->const_array(mfi);
+                Array4<Real const> const& ez_c = Etmp[2]->const_array(mfi);
+
+                const Box& bx = mfi.fabbox();
+                amrex::ParallelFor(bx,
+                [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
+                {
+                    warpx_interp_nd_efield_x(j,k,l, ex_aux, ex_fp, ex_cp, ex_c);
+                    warpx_interp_nd_efield_y(j,k,l, ey_aux, ey_fp, ey_cp, ey_c);
+                    warpx_interp_nd_efield_z(j,k,l, ez_aux, ez_fp, ez_cp, ez_c);
+                });
+            }
+        }
+    }
+}
+
+void
+WarpX::UpdateAuxilaryDataSameType ()
+{
     for (int lev = 1; lev <= finest_level; ++lev)
     {
         const auto& crse_period = Geom(lev-1).periodicity();
diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H
index 093323ec3..5da867c9f 100644
--- a/Source/Parallelization/WarpXComm_K.H
+++ b/Source/Parallelization/WarpXComm_K.H
@@ -158,4 +158,489 @@ void warpx_interp_efield_z (int j, int k, int l,
 #endif
 }
 
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_bfield_x (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Bxa,
+                               amrex::Array4<amrex::Real const> const& Bxf,
+                               amrex::Array4<amrex::Real const> const& Bxc,
+                               amrex::Array4<amrex::Real const> const& Bxg)
+{
+    using namespace amrex;
+
+    int jg = amrex::coarsen(j,2);
+    Real wx = (j == jg*2) ? 0.0 : 0.5;
+    Real owx = 1.0-wx;
+
+    int kg = amrex::coarsen(k,2);
+    Real wy = (k == kg*2) ? 0.0 : 0.5;
+    Real owy = 1.0-wy;
+
+#if (AMREX_SPACEDIM == 2)
+
+    // interp from coarse nodal to fine nodal
+    Real bg = owx * owy * Bxg(jg  ,kg  ,0)
+        +     owx *  wy * Bxg(jg  ,kg+1,0)
+        +      wx * owy * Bxg(jg+1,kg  ,0)
+        +      wx *  wy * Bxg(jg+1,kg+1,0);
+
+    // interp from coarse staggered to fine nodal
+    wy = 0.5-wy;  owy = 1.0-wy;
+    Real bc = owx * owy * Bxc(jg  ,kg  ,0)
+        +     owx *  wy * Bxc(jg  ,kg-1,0)
+        +      wx * owy * Bxc(jg+1,kg  ,0)
+        +      wx *  wy * Bxc(jg+1,kg-1,0);
+
+    // interp from fine staggered to fine nodal
+    Real bf = 0.5*(Bxf(j,k-1,0) + Bxf(j,k,0));
+
+#else
+
+    int lg = amrex::coarsen(l,2);
+    Real wz = (l == lg*2) ? 0.0 : 0.5;
+    Real owz = 1.0-wz;
+
+    // interp from coarse nodal to fine nodal
+    Real bg = owx * owy * owz * Bxg(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Bxg(jg+1,kg  ,lg  )
+        +     owx *  wy * owz * Bxg(jg  ,kg+1,lg  )
+        +      wx *  wy * owz * Bxg(jg+1,kg+1,lg  )
+        +     owx * owy *  wz * Bxg(jg  ,kg  ,lg+1)
+        +      wx * owy *  wz * Bxg(jg+1,kg  ,lg+1)
+        +     owx *  wy *  wz * Bxg(jg  ,kg+1,lg+1)
+        +      wx *  wy *  wz * Bxg(jg+1,kg+1,lg+1);
+
+    // interp from coarse staggered to fine nodal
+    wy = 0.5-wy;  owy = 1.0-wy;
+    wz = 0.5-wz;  owz = 1.0-wz;
+    Real bc = owx * owy * owz * Bxc(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Bxc(jg+1,kg  ,lg  )
+        +     owx *  wy * owz * Bxc(jg  ,kg-1,lg  )
+        +      wx *  wy * owz * Bxc(jg+1,kg-1,lg  )
+        +     owx * owy *  wz * Bxc(jg  ,kg  ,lg-1)
+        +      wx * owy *  wz * Bxc(jg+1,kg  ,lg-1)
+        +     owx *  wy *  wz * Bxc(jg  ,kg-1,lg-1)
+        +      wx *  wy *  wz * Bxc(jg+1,kg-1,lg-1);
+
+    // interp from fine stagged to fine nodal
+    Real bf = 0.25*(Bxf(j,k-1,l-1) + Bxf(j,k,l-1) + Bxf(j,k-1,l) + Bxf(j,k,l));
+#endif
+
+    Bxa(j,k,l) = bg + (bf-bc);
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_bfield_y (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Bya,
+                               amrex::Array4<amrex::Real const> const& Byf,
+                               amrex::Array4<amrex::Real const> const& Byc,
+                               amrex::Array4<amrex::Real const> const& Byg)
+{
+    using namespace amrex;
+
+    int jg = amrex::coarsen(j,2);
+    Real wx = (j == jg*2) ? 0.0 : 0.5;
+    Real owx = 1.0-wx;
+
+    int kg = amrex::coarsen(k,2);
+    Real wy = (k == kg*2) ? 0.0 : 0.5;
+    Real owy = 1.0-wy;
+
+#if (AMREX_SPACEDIM == 2)
+
+    // interp from coarse nodal to fine nodal
+    Real bg = owx * owy * Byg(jg  ,kg  ,0)
+        +     owx *  wy * Byg(jg  ,kg+1,0)
+        +      wx * owy * Byg(jg+1,kg  ,0)
+        +      wx *  wy * Byg(jg+1,kg+1,0);
+
+    // interp from coarse stagged (cell-centered for By) to fine nodal
+    wx = 0.5-wx;  owx = 1.0-wx;
+    wy = 0.5-wy;  owy = 1.0-wy;
+    Real bc = owx * owy * Byc(jg  ,kg  ,0)
+        +     owx *  wy * Byc(jg  ,kg-1,0)
+        +      wx * owy * Byc(jg-1,kg  ,0)
+        +      wx *  wy * Byc(jg-1,kg-1,0);
+
+    // interp form fine stagger (cell-centered for By) to fine nodal
+    Real bf = 0.25*(Byf(j,k,0) + Byf(j-1,k,0) + Byf(j,k-1,0) + Byf(j-1,k-1,0));
+
+#else
+
+    int lg = amrex::coarsen(l,2);
+    Real wz = (l == lg*2) ? 0.0 : 0.5;
+    Real owz = 1.0-wz;
+
+    // interp from coarse nodal to fine nodal
+    Real bg = owx * owy * owz * Byg(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Byg(jg+1,kg  ,lg  )
+        +     owx *  wy * owz * Byg(jg  ,kg+1,lg  )
+        +      wx *  wy * owz * Byg(jg+1,kg+1,lg  )
+        +     owx * owy *  wz * Byg(jg  ,kg  ,lg+1)
+        +      wx * owy *  wz * Byg(jg+1,kg  ,lg+1)
+        +     owx *  wy *  wz * Byg(jg  ,kg+1,lg+1)
+        +      wx *  wy *  wz * Byg(jg+1,kg+1,lg+1);
+
+    // interp from coarse staggered to fine nodal
+    wx = 0.5-wx;  owx = 1.0-wx;
+    wz = 0.5-wz;  owz = 1.0-wz;
+    Real bc = owx * owy * owz * Byc(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Byc(jg-1,kg  ,lg  )
+        +     owx *  wy * owz * Byc(jg  ,kg+1,lg  )
+        +      wx *  wy * owz * Byc(jg-1,kg+1,lg  )
+        +     owx * owy *  wz * Byc(jg  ,kg  ,lg-1)
+        +      wx * owy *  wz * Byc(jg-1,kg  ,lg-1)
+        +     owx *  wy *  wz * Byc(jg  ,kg+1,lg-1)
+        +      wx *  wy *  wz * Byc(jg-1,kg+1,lg-1);
+
+    // interp from fine stagged to fine nodal
+    Real bf = 0.25*(Byf(j-1,k,l-1) + Byf(j,k,l-1) + Byf(j-1,k,l) + Byf(j,k,l));
+
+#endif
+
+    Bya(j,k,l) = bg + (bf-bc);
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_bfield_z (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Bza,
+                               amrex::Array4<amrex::Real const> const& Bzf,
+                               amrex::Array4<amrex::Real const> const& Bzc,
+                               amrex::Array4<amrex::Real const> const& Bzg)
+{
+    using namespace amrex;
+
+    int jg = amrex::coarsen(j,2);
+    Real wx = (j == jg*2) ? 0.0 : 0.5;
+    Real owx = 1.0-wx;
+
+    int kg = amrex::coarsen(k,2);
+    Real wy = (k == kg*2) ? 0.0 : 0.5;
+    Real owy = 1.0-wy;
+
+#if (AMREX_SPACEDIM == 2)
+
+    // interp from coarse nodal to fine nodal
+    Real bg = owx * owy * Bzg(jg  ,kg  ,0)
+        +     owx *  wy * Bzg(jg  ,kg+1,0)
+        +      wx * owy * Bzg(jg+1,kg  ,0)
+        +      wx *  wy * Bzg(jg+1,kg+1,0);
+
+    // interp from coarse staggered to fine nodal
+    wx = 0.5-wx;  owx = 1.0-wx;
+    Real bc = owx * owy * Bzc(jg  ,kg  ,0)
+        +     owx *  wy * Bzc(jg  ,kg+1,0)
+        +      wx * owy * Bzc(jg-1,kg  ,0)
+        +      wx *  wy * Bzc(jg-1,kg+1,0);
+
+    // interp from fine staggered to fine nodal
+    Real bf = 0.5*(Bzf(j-1,k,0) + Bzf(j,k,0));
+
+#else
+
+    int lg = amrex::coarsen(l,2);
+    Real wz = (l == lg*2) ? 0.0 : 0.5;
+    Real owz = 1.0-wz;
+
+    // interp from coarse nodal to fine nodal
+    Real bg = owx * owy * owz * Bzg(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Bzg(jg+1,kg  ,lg  )
+        +     owx *  wy * owz * Bzg(jg  ,kg+1,lg  )
+        +      wx *  wy * owz * Bzg(jg+1,kg+1,lg  )
+        +     owx * owy *  wz * Bzg(jg  ,kg  ,lg+1)
+        +      wx * owy *  wz * Bzg(jg+1,kg  ,lg+1)
+        +     owx *  wy *  wz * Bzg(jg  ,kg+1,lg+1)
+        +      wx *  wy *  wz * Bzg(jg+1,kg+1,lg+1);
+
+    // interp from coarse staggered to fine nodal
+    wx = 0.5-wx;  owx = 1.0-wx;
+    wy = 0.5-wy;  owy = 1.0-wy;
+    Real bc = owx * owy * owz * Bzc(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Bzc(jg-1,kg  ,lg  )
+        +     owx *  wy * owz * Bzc(jg  ,kg-1,lg  )
+        +      wx *  wy * owz * Bzc(jg-1,kg-1,lg  )
+        +     owx * owy *  wz * Bzc(jg  ,kg  ,lg+1)
+        +      wx * owy *  wz * Bzc(jg-1,kg  ,lg+1)
+        +     owx *  wy *  wz * Bzc(jg  ,kg-1,lg+1)
+        +      wx *  wy *  wz * Bzc(jg-1,kg-1,lg+1);
+
+    // interp from fine stagged to fine nodal
+    Real bf = 0.25*(Bzf(j-1,k-1,l) + Bzf(j,k-1,l) + Bzf(j-1,k,l) + Bzf(j,k,l));
+
+#endif
+
+    Bza(j,k,l) = bg + (bf-bc);
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_bfield_x (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Bxa,
+                               amrex::Array4<amrex::Real const> const& Bxf)
+{
+#if (AMREX_SPACEDIM == 2)
+    Bxa(j,k,0) = 0.5*(Bxf(j,k-1,0) + Bxf(j,k,0));
+#else
+    Bxa(j,k,l) = 0.25*(Bxf(j,k-1,l-1) + Bxf(j,k,l-1) + Bxf(j,k-1,l) + Bxf(j,k,l));
+#endif
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_bfield_y (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Bya,
+                               amrex::Array4<amrex::Real const> const& Byf)
+{
+#if (AMREX_SPACEDIM == 2)
+    Bya(j,k,0) = 0.25*(Byf(j,k,0) + Byf(j-1,k,0) + Byf(j,k-1,0) + Byf(j-1,k-1,0));
+#else
+    Bya(j,k,l) = 0.25*(Byf(j-1,k,l-1) + Byf(j,k,l-1) + Byf(j-1,k,l) + Byf(j,k,l));
+#endif
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_bfield_z (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Bza,
+                               amrex::Array4<amrex::Real const> const& Bzf)
+{
+#if (AMREX_SPACEDIM == 2)
+    Bza(j,k,0) = 0.5*(Bzf(j-1,k,0) + Bzf(j,k,0));
+#else
+    Bza(j,k,l) = 0.25*(Bzf(j-1,k-1,l) + Bzf(j,k-1,l) + Bzf(j-1,k,l) + Bzf(j,k,l));
+#endif
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_efield_x (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Exa,
+                               amrex::Array4<amrex::Real const> const& Exf,
+                               amrex::Array4<amrex::Real const> const& Exc,
+                               amrex::Array4<amrex::Real const> const& Exg)
+{
+    using namespace amrex;
+
+    int jg = amrex::coarsen(j,2);
+    Real wx = (j == jg*2) ? 0.0 : 0.5;
+    Real owx = 1.0-wx;
+
+    int kg = amrex::coarsen(k,2);
+    Real wy = (k == kg*2) ? 0.0 : 0.5;
+    Real owy = 1.0-wy;
+
+#if (AMREX_SPACEDIM == 2)
+
+    // interp from coarse nodal to fine nodal
+    Real eg = owx * owy * Exg(jg  ,kg  ,0)
+        +     owx *  wy * Exg(jg  ,kg+1,0)
+        +      wx * owy * Exg(jg+1,kg  ,0)
+        +      wx *  wy * Exg(jg+1,kg+1,0);
+
+    // interp from coarse staggered to fine nodal
+    wx = 0.5-wx;  owx = 1.0-wx;
+    Real ec = owx * owy * Exc(jg  ,kg  ,0)
+        +     owx *  wy * Exc(jg  ,kg+1,0)
+        +      wx * owy * Exc(jg-1,kg  ,0)
+        +      wx *  wy * Exc(jg-1,kg+1,0);
+
+    // interp from fine staggered to fine nodal
+    Real ef = 0.5*(Exf(j-1,k,0) + Exf(j,k,0));
+
+#else
+
+    int lg = amrex::coarsen(l,2);
+    Real wz = (l == lg*2) ? 0.0 : 0.5;
+    Real owz = 1.0-wz;
+
+    // interp from coarse nodal to fine nodal
+    Real eg = owx * owy * owz * Exg(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Exg(jg+1,kg  ,lg  )
+        +     owx *  wy * owz * Exg(jg  ,kg+1,lg  )
+        +      wx *  wy * owz * Exg(jg+1,kg+1,lg  )
+        +     owx * owy *  wz * Exg(jg  ,kg  ,lg+1)
+        +      wx * owy *  wz * Exg(jg+1,kg  ,lg+1)
+        +     owx *  wy *  wz * Exg(jg  ,kg+1,lg+1)
+        +      wx *  wy *  wz * Exg(jg+1,kg+1,lg+1);
+
+    // interp from coarse staggered to fine nodal
+    wx = 0.5-wx;  owx = 1.0-wx;
+    Real ec = owx * owy * owz * Exc(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Exc(jg-1,kg  ,lg  )
+        +     owx *  wy * owz * Exc(jg  ,kg+1,lg  )
+        +      wx *  wy * owz * Exc(jg-1,kg+1,lg  )
+        +     owx * owy *  wz * Exc(jg  ,kg  ,lg+1)
+        +      wx * owy *  wz * Exc(jg-1,kg  ,lg+1)
+        +     owx *  wy *  wz * Exc(jg  ,kg+1,lg+1)
+        +      wx *  wy *  wz * Exc(jg-1,kg+1,lg+1);
+
+    // interp from fine staggered to fine nodal
+    Real ef = 0.5*(Exf(j-1,k,l) + Exf(j,k,l));
+
+#endif
+
+    Exa(j,k,l) = eg + (ef-ec);
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_efield_y (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Eya,
+                               amrex::Array4<amrex::Real const> const& Eyf,
+                               amrex::Array4<amrex::Real const> const& Eyc,
+                               amrex::Array4<amrex::Real const> const& Eyg)
+{
+    using namespace amrex;
+
+    int jg = amrex::coarsen(j,2);
+    Real wx = (j == jg*2) ? 0.0 : 0.5;
+    Real owx = 1.0-wx;
+
+    int kg = amrex::coarsen(k,2);
+    Real wy = (k == kg*2) ? 0.0 : 0.5;
+    Real owy = 1.0-wy;
+
+#if (AMREX_SPACEDIM == 2)
+
+    // interp from coarse nodal and coarse staggered to fine nodal
+    Real eg = owx * owy * (Eyg(jg  ,kg  ,0) + Eyc(jg  ,kg  ,0))
+        +     owx *  wy * (Eyg(jg  ,kg+1,0) + Eyc(jg  ,kg+1,0))
+        +      wx * owy * (Eyg(jg+1,kg  ,0) + Eyc(jg+1,kg  ,0))
+        +      wx *  wy * (Eyg(jg+1,kg+1,0) + Eyc(jg+1,kg+1,0));
+    Real ec = 0.0;
+
+    // interp from fine staggered to fine nodal
+    Real ef = Eyf(j,k,0);
+
+#else
+
+    int lg = amrex::coarsen(l,2);
+    Real wz = (l == lg*2) ? 0.0 : 0.5;
+    Real owz = 1.0-wz;
+
+    // interp from coarse nodal to fine nodal
+    Real eg = owx * owy * owz * Eyg(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Eyg(jg+1,kg  ,lg  )
+        +     owx *  wy * owz * Eyg(jg  ,kg+1,lg  )
+        +      wx *  wy * owz * Eyg(jg+1,kg+1,lg  )
+        +     owx * owy *  wz * Eyg(jg  ,kg  ,lg+1)
+        +      wx * owy *  wz * Eyg(jg+1,kg  ,lg+1)
+        +     owx *  wy *  wz * Eyg(jg  ,kg+1,lg+1)
+        +      wx *  wy *  wz * Eyg(jg+1,kg+1,lg+1);
+
+    // interp from coarse staggered to fine nodal
+    wy = 0.5-wy;  owy = 1.0-wy;
+    Real ec = owx * owy * owz * Eyc(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Eyc(jg+1,kg  ,lg  )
+        +     owx *  wy * owz * Eyc(jg  ,kg-1,lg  )
+        +      wx *  wy * owz * Eyc(jg+1,kg-1,lg  )
+        +     owx * owy *  wz * Eyc(jg  ,kg  ,lg+1)
+        +      wx * owy *  wz * Eyc(jg+1,kg  ,lg+1)
+        +     owx *  wy *  wz * Eyc(jg  ,kg-1,lg+1)
+        +      wx *  wy *  wz * Eyc(jg+1,kg-1,lg+1);
+
+    // interp from fine staggered to fine nodal
+    Real ef = 0.5*(Eyf(j,k-1,l) + Eyf(j,k,l));
+
+#endif
+
+    Eya(j,k,l) = eg + (ef-ec);
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_efield_z (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Eza,
+                               amrex::Array4<amrex::Real const> const& Ezf,
+                               amrex::Array4<amrex::Real const> const& Ezc,
+                               amrex::Array4<amrex::Real const> const& Ezg)
+{
+    using namespace amrex;
+
+    int jg = amrex::coarsen(j,2);
+    Real wx = (j == jg*2) ? 0.0 : 0.5;
+    Real owx = 1.0-wx;
+
+    int kg = amrex::coarsen(k,2);
+    Real wy = (k == kg*2) ? 0.0 : 0.5;
+    Real owy = 1.0-wy;
+
+#if (AMREX_SPACEDIM == 2)
+
+    // interp from coarse nodal to fine nodal
+    Real eg = owx * owy * Ezg(jg  ,kg  ,0)
+        +     owx *  wy * Ezg(jg  ,kg+1,0)
+        +      wx * owy * Ezg(jg+1,kg  ,0)
+        +      wx *  wy * Ezg(jg+1,kg+1,0);
+
+    // interp from coarse stagged to fine nodal
+    wy = 0.5-wy;  owy = 1.0-wy;
+    Real ec = owx * owy * Ezc(jg  ,kg  ,0)
+        +     owx *  wy * Ezc(jg  ,kg-1,0)
+        +      wx * owy * Ezc(jg+1,kg  ,0)
+        +      wx *  wy * Ezc(jg+1,kg-1,0);
+
+    // interp from fine staggered to fine nodal
+    Real ef = 0.5*(Ezf(j,k-1,0) + Ezf(j,k,0));
+
+#else
+
+    int lg = amrex::coarsen(l,2);
+    Real wz = (l == lg*2) ? 0.0 : 0.5;
+    Real owz = 1.0-wz;
+
+    // interp from coarse nodal to fine nodal
+    Real eg = owx * owy * owz * Ezg(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Ezg(jg+1,kg  ,lg  )
+        +     owx *  wy * owz * Ezg(jg  ,kg+1,lg  )
+        +      wx *  wy * owz * Ezg(jg+1,kg+1,lg  )
+        +     owx * owy *  wz * Ezg(jg  ,kg  ,lg+1)
+        +      wx * owy *  wz * Ezg(jg+1,kg  ,lg+1)
+        +     owx *  wy *  wz * Ezg(jg  ,kg+1,lg+1)
+        +      wx *  wy *  wz * Ezg(jg+1,kg+1,lg+1);
+
+    // interp from coarse staggered to fine nodal
+    wz = 0.5-wz;  owz = 1.0-wz;
+    Real ec = owx * owy * owz * Ezc(jg  ,kg  ,lg  )
+        +      wx * owy * owz * Ezc(jg+1,kg  ,lg  )
+        +     owx *  wy * owz * Ezc(jg  ,kg+1,lg  )
+        +      wx *  wy * owz * Ezc(jg+1,kg+1,lg  )
+        +     owx * owy *  wz * Ezc(jg  ,kg  ,lg-1)
+        +      wx * owy *  wz * Ezc(jg+1,kg  ,lg-1)
+        +     owx *  wy *  wz * Ezc(jg  ,kg+1,lg-1)
+        +      wx *  wy *  wz * Ezc(jg+1,kg+1,lg-1);
+
+    // interp from fine staggered to fine nodal
+    Real ef = 0.5*(Ezf(j,k,l-1) + Ezf(j,k,l));
+
+#endif
+
+    Eza(j,k,l) = eg + (ef-ec);
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_efield_x (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Exa,
+                               amrex::Array4<amrex::Real const> const& Exf)
+{
+    Exa(j,k,l) = 0.5*(Exf(j-1,k,l) + Exf(j,k,l));
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_efield_y (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Eya,
+                               amrex::Array4<amrex::Real const> const& Eyf)
+{
+#if (AMREX_SPACEDIM == 2)
+    Eya(j,k,0) = Eyf(j,k,0);
+#else
+    Eya(j,k,l) = 0.5*(Eyf(j,k-1,l) + Eyf(j,k,l));
+#endif
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_nd_efield_z (int j, int k, int l,
+                               amrex::Array4<amrex::Real> const& Eza,
+                               amrex::Array4<amrex::Real const> const& Ezf)
+{
+#if (AMREX_SPACEDIM == 2)
+    Eza(j,k,0) = 0.5*(Ezf(j,k-1,0) + Ezf(j,k,0));
+#else
+    Eza(j,k,l) = 0.5*(Ezf(j,k,l-1) + Ezf(j,k,l));
+#endif
+}
+
 #endif
diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp
index 5441755f5..2ae167283 100644
--- a/Source/Parallelization/WarpXRegrid.cpp
+++ b/Source/Parallelization/WarpXRegrid.cpp
@@ -91,7 +91,7 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa
 
         // Aux patch
 
-        if (lev == 0)
+        if (lev == 0 && Bfield_aux[0][0]->ixType() == Bfield_fp[0][0]->ixType())
         {
             for (int idim = 0; idim < 3; ++idim) {
                 Bfield_aux[lev][idim].reset(new MultiFab(*Bfield_fp[lev][idim], amrex::make_alias, 0, Bfield_aux[lev][idim]->nComp()));
diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H
index c19e166b0..75dda58af 100644
--- a/Source/Particles/MultiParticleContainer.H
+++ b/Source/Particles/MultiParticleContainer.H
@@ -218,8 +218,8 @@ protected:
 
 #ifdef WARPX_QED
     // The QED engines
-    BreitWheelerEngine bw_engine;
-    QuantumSynchrotronEngine qs_engine;
+    std::shared_ptr<BreitWheelerEngine> shr_p_bw_engine;
+    std::shared_ptr<QuantumSynchrotronEngine> shr_p_qs_engine;
     //_______________________________
 
     /**
diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp
index 566422a7a..54a4396c1 100644
--- a/Source/Particles/MultiParticleContainer.cpp
+++ b/Source/Particles/MultiParticleContainer.cpp
@@ -744,16 +744,17 @@ MultiParticleContainer::doFieldIonization ()
 #ifdef WARPX_QED
 void MultiParticleContainer::InitQED ()
 {
+    shr_p_qs_engine = std::make_shared<QuantumSynchrotronEngine>();
+    shr_p_bw_engine = std::make_shared<BreitWheelerEngine>();
+
     for (auto& pc : allcontainers) {
         if(pc->has_quantum_sync()){
             pc->set_quantum_sync_engine_ptr
-                (std::make_shared<QuantumSynchrotronEngine>(qs_engine));
-            someone_has_quantum_sync = true;
+                (shr_p_qs_engine);
         }
         if(pc->has_breit_wheeler()){
             pc->set_breit_wheeler_engine_ptr
-                (std::make_shared<BreitWheelerEngine>(bw_engine));
-            someone_has_breit_wheeler = true;
+                (shr_p_bw_engine);
         }
     }
 
@@ -781,8 +782,8 @@ void MultiParticleContainer::InitQuantumSync ()
 
 
     if(generate_table && ParallelDescriptor::IOProcessor()){
-        qs_engine.compute_lookup_tables(ctrl);
-        Vector<char> all_data = qs_engine.export_lookup_tables_data();
+        shr_p_qs_engine->compute_lookup_tables(ctrl);
+        Vector<char> all_data = shr_p_qs_engine->export_lookup_tables_data();
         WarpXUtilIO::WriteBinaryDataOnFile(filename, all_data);
    }
     ParallelDescriptor::Barrier();
@@ -794,10 +795,10 @@ void MultiParticleContainer::InitQuantumSync ()
     //No need to initialize from raw data for the processor that
     //has just generated the table
     if(!generate_table || !ParallelDescriptor::IOProcessor()){
-        qs_engine.init_lookup_tables_from_raw_data(table_data);
+        shr_p_qs_engine->init_lookup_tables_from_raw_data(table_data);
     }
 
-    if(!qs_engine.are_lookup_tables_initialized())
+    if(!shr_p_qs_engine->are_lookup_tables_initialized())
         amrex::Error("Table initialization has failed!\n");
 }
 
@@ -816,8 +817,8 @@ void MultiParticleContainer::InitBreitWheeler ()
     //_________________________________________________
 
     if(generate_table && ParallelDescriptor::IOProcessor()){
-        bw_engine.compute_lookup_tables(ctrl);
-        Vector<char> all_data = bw_engine.export_lookup_tables_data();
+        shr_p_bw_engine->compute_lookup_tables(ctrl);
+        Vector<char> all_data =shr_p_bw_engine->export_lookup_tables_data();
         WarpXUtilIO::WriteBinaryDataOnFile(filename, all_data);
     }
     ParallelDescriptor::Barrier();
@@ -829,10 +830,10 @@ void MultiParticleContainer::InitBreitWheeler ()
     //No need to initialize from raw data for the processor that
     //has just generated the table
     if(!generate_table || !ParallelDescriptor::IOProcessor()){
-        bw_engine.init_lookup_tables_from_raw_data(table_data);
+        shr_p_bw_engine->init_lookup_tables_from_raw_data(table_data);
     }
 
-    if(!bw_engine.are_lookup_tables_initialized())
+    if(!shr_p_bw_engine->are_lookup_tables_initialized())
         amrex::Error("Table initialization has failed!\n");
 }
 
@@ -840,7 +841,7 @@ std::tuple<bool,std::string,PicsarQuantumSynchrotronCtrl>
 MultiParticleContainer::ParseQuantumSyncParams ()
 {
     PicsarQuantumSynchrotronCtrl ctrl =
-        qs_engine.get_default_ctrl();
+        shr_p_qs_engine->get_default_ctrl();
     bool generate_table{false};
     std::string table_name;
 
@@ -898,7 +899,7 @@ std::tuple<bool,std::string,PicsarBreitWheelerCtrl>
 MultiParticleContainer::ParseBreitWheelerParams ()
 {
     PicsarBreitWheelerCtrl ctrl =
-        bw_engine.get_default_ctrl();
+        shr_p_bw_engine->get_default_ctrl();
     bool generate_table{false};
     std::string table_name;
 
diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp
index 5bd7eb372..ea55e04c6 100644
--- a/Source/Particles/PhotonParticleContainer.cpp
+++ b/Source/Particles/PhotonParticleContainer.cpp
@@ -45,9 +45,8 @@ void PhotonParticleContainer::InitData()
 {
     AddParticles(0); // Note - add on level 0
 
-    if (maxLevel() > 0) {
-        Redistribute();  // We then redistribute
-    }
+    Redistribute();  // We then redistribute
+
 }
 
 void
diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H
index 269171a5f..7d26e7af5 100644
--- a/Source/Utils/WarpXAlgorithmSelection.H
+++ b/Source/Utils/WarpXAlgorithmSelection.H
@@ -34,9 +34,9 @@ struct ChargeDepositionAlgo {
 };
 
 struct GatheringAlgo {
-    // Only the Standard algorithm is implemented
     enum {
-         Standard = 0
+         EnergyConserving = 0,
+         MomentumConserving
     };
 };
 
diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp
index be9cdd118..4b66a0809 100644
--- a/Source/Utils/WarpXAlgorithmSelection.cpp
+++ b/Source/Utils/WarpXAlgorithmSelection.cpp
@@ -34,8 +34,9 @@ const std::map<std::string, int> charge_deposition_algo_to_int = {
 };
 
 const std::map<std::string, int> gathering_algo_to_int = {
-    {"standard",   GatheringAlgo::Standard },
-    {"default",    GatheringAlgo::Standard }
+    {"energy-conserving",   GatheringAlgo::EnergyConserving },
+    {"momentum-conserving", GatheringAlgo::MomentumConserving },
+    {"default",             GatheringAlgo::EnergyConserving }
 };
 
 
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 29b89686e..6ecc584a2 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -204,6 +204,8 @@ public:
     // This function does aux(lev) = fp(lev) + I(aux(lev-1)-cp(lev)).
     // Caller must make sure fp and cp have ghost cells filled.
     void UpdateAuxilaryData ();
+    void UpdateAuxilaryDataStagToNodal ();
+    void UpdateAuxilaryDataSameType ();
 
     // Fill boundary cells including coarse/fine boundaries
     void FillBoundaryB ();
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index d94541f17..9ba0741bb 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -541,9 +541,13 @@ WarpX::ReadParameters ()
         ParmParse pp("algo");
         current_deposition_algo = GetAlgorithmInteger(pp, "current_deposition");
         charge_deposition_algo = GetAlgorithmInteger(pp, "charge_deposition");
-        field_gathering_algo = GetAlgorithmInteger(pp, "field_gathering");
         particle_pusher_algo = GetAlgorithmInteger(pp, "particle_pusher");
         maxwell_fdtd_solver_id = GetAlgorithmInteger(pp, "maxwell_fdtd_solver");
+        field_gathering_algo = GetAlgorithmInteger(pp, "field_gathering");
+        if (field_gathering_algo == GatheringAlgo::MomentumConserving) {
+            // Use same shape factors in all directions, for gathering
+            l_lower_order_in_v = false;
+        }
     }
 
 #ifdef WARPX_USE_PSATD
@@ -811,16 +815,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
     ncomps = n_rz_azimuthal_modes*2 - 1;
 #endif
 
+    bool aux_is_nodal = (field_gathering_algo == GatheringAlgo::MomentumConserving);
+    IntVect ngextra(static_cast<int>(aux_is_nodal and !do_nodal));
+
     //
     // The fine patch
     //
-    Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE));
-    Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE));
-    Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE));
+    Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE+ngextra));
+    Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE+ngextra));
+    Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE+ngextra));
 
-    Efield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE));
-    Efield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE));
-    Efield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE));
+    Efield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE+ngextra));
+    Efield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE+ngextra));
+    Efield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE+ngextra));
 
     current_fp[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,ncomps,ngJ));
     current_fp[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ));
@@ -867,7 +874,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
     //
     // The Aux patch (i.e., the full solution)
     //
-    if (lev == 0)
+    if (aux_is_nodal and !do_nodal)
+    {
+        // Create aux multifabs on Nodal Box Array
+        BoxArray const nba = amrex::convert(ba,IntVect::TheNodeVector());
+        Bfield_aux[lev][0].reset( new MultiFab(nba,dm,ncomps,ngE));
+        Bfield_aux[lev][1].reset( new MultiFab(nba,dm,ncomps,ngE));
+        Bfield_aux[lev][2].reset( new MultiFab(nba,dm,ncomps,ngE));
+
+        Efield_aux[lev][0].reset( new MultiFab(nba,dm,ncomps,ngE));
+        Efield_aux[lev][1].reset( new MultiFab(nba,dm,ncomps,ngE));
+        Efield_aux[lev][2].reset( new MultiFab(nba,dm,ncomps,ngE));
+    }
+    else if (lev == 0)
     {
         for (int idir = 0; idir < 3; ++idir) {
             Efield_aux[lev][idir].reset(new MultiFab(*Efield_fp[lev][idir], amrex::make_alias, 0, ncomps));
@@ -948,15 +967,25 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
         cba.coarsen(refRatio(lev-1));
 
         if (n_field_gather_buffer > 0 || mypc->nSpeciesGatherFromMainGrid() > 0) {
-            // Create the MultiFabs for B
-            Bfield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE));
-            Bfield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,ncomps,ngE));
-            Bfield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,ncomps,ngE));
-
-            // Create the MultiFabs for E
-            Efield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,ncomps,ngE));
-            Efield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE));
-            Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE));
+            if (aux_is_nodal) {
+                BoxArray const& cnba = amrex::convert(cba,IntVect::TheNodeVector());
+                Bfield_cax[lev][0].reset( new MultiFab(cnba,dm,ncomps,ngE));
+                Bfield_cax[lev][1].reset( new MultiFab(cnba,dm,ncomps,ngE));
+                Bfield_cax[lev][2].reset( new MultiFab(cnba,dm,ncomps,ngE));
+                Efield_cax[lev][0].reset( new MultiFab(cnba,dm,ncomps,ngE));
+                Efield_cax[lev][1].reset( new MultiFab(cnba,dm,ncomps,ngE));
+                Efield_cax[lev][2].reset( new MultiFab(cnba,dm,ncomps,ngE));
+            } else {
+                // Create the MultiFabs for B
+                Bfield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE));
+                Bfield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,ncomps,ngE));
+                Bfield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,ncomps,ngE));
+
+                // Create the MultiFabs for E
+                Efield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,ncomps,ngE));
+                Efield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE));
+                Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE));
+            }
 
             gather_buffer_masks[lev].reset( new iMultiFab(ba, dm, ncomps, 1) );
             // Gather buffer masks have 1 ghost cell, because of the fact