Merge branch 'dev' of https://github.com/ECP-WarpX/WarpX into doNotDepositCurrent

7 files changed, 888 insertions, 95 deletions

diff --git a/Source/Parallelization/InterpolateCurrentFineToCoarse.H b/Source/Parallelization/InterpolateCurrentFineToCoarse.H
index 148b725d0..cbbcdfab5 100644
--- a/Source/Parallelization/InterpolateCurrentFineToCoarse.H
+++ b/Source/Parallelization/InterpolateCurrentFineToCoarse.H
@@ -56,6 +56,8 @@ public:
         int const k
     ) const noexcept
     {
+        using namespace amrex;
+
         auto const & fine_unsafe = m_fine;   // out-of-bounds access not secured with zero-values yet
         auto const & coarse      = m_coarse; // out-of-bounds access not secured but will also not occur
 
@@ -71,29 +73,29 @@ public:
         int const kk = k * m_refinement_ratio;
 #if AMREX_SPACEDIM == 2
         if (IDim == 0) {
-            coarse(i, j, k) = 0.25 * (
+            coarse(i, j, k) = 0.25_rt * (
                 fine(ii, jj, kk) + fine(ii + 1, jj, kk) +
-                0.5  * (
+                0.5_rt  * (
                     fine(ii, jj - 1, kk) + fine(ii + 1, jj - 1, kk) +
                     fine(ii, jj + 1, kk) + fine(ii + 1, jj + 1, kk)
                 )
             );
         } else if (IDim == 2) {
-            coarse(i, j, k) = 0.25 * (
+            coarse(i, j, k) = 0.25_rt * (
                 fine(ii, jj, kk) + fine(ii, jj + 1, kk) +
-                0.5  * (
+                0.5_rt  * (
                     fine(ii - 1, jj, kk) + fine(ii - 1, jj + 1, kk) +
                     fine(ii + 1, jj, kk) + fine(ii + 1, jj + 1, kk)
                 )
             );
         } else {
-            coarse(i, j, k) = 0.25 * (
+            coarse(i, j, k) = 0.25_rt * (
                 fine(ii, jj, kk) +
-                0.5  * (
+                0.5_rt  * (
                     fine(ii - 1, jj    , kk) + fine(ii + 1, jj    , kk) +
                     fine(ii    , jj - 1, kk) + fine(ii    , jj + 1, kk)
                 ) +
-                0.25 * (
+                0.25_rt * (
                     fine(ii - 1, jj - 1, kk) + fine(ii + 1, jj - 1, kk) +
                     fine(ii - 1, jj + 1, kk) + fine(ii + 1, jj + 1, kk)
                 )
@@ -101,64 +103,64 @@ public:
         }
 #elif AMREX_SPACEDIM == 3
         if (IDim == 0) {
-            coarse(i,j,k) = 0.125 * (
+            coarse(i,j,k) = 0.125_rt * (
                 fine(ii  , jj, kk) +
-                0.5  * (
+                0.5_rt  * (
                     fine(ii  , jj-1, kk  ) + fine(ii  , jj+1, kk  ) +
                     fine(ii  , jj  , kk-1) + fine(ii  , jj  , kk+1)
                 ) +
-                0.25 * (
+                0.25_rt * (
                     fine(ii  , jj-1, kk-1) + fine(ii  , jj+1, kk-1) +
                     fine(ii  , jj-1, kk+1) + fine(ii  , jj+1, kk+1)
                 ) +
                 fine(ii+1, jj, kk) +
-                0.5  * (
+                0.5_rt  * (
                     fine(ii+1, jj-1, kk  ) + fine(ii+1, jj+1, kk  ) +
                     fine(ii+1, jj  , kk-1) + fine(ii+1, jj  , kk+1)
                 ) +
-                0.25 * (
+                0.25_rt * (
                     fine(ii+1, jj-1, kk-1) + fine(ii+1, jj+1, kk-1) +
                     fine(ii+1, jj-1, kk+1) + fine(ii+1, jj+1, kk+1)
                 )
             );
         } else if (IDim == 1) {
-            coarse(i, j, k) = 0.125 * (
+            coarse(i, j, k) = 0.125_rt * (
                 fine(ii, jj , kk) +
-                0.5  * (
+                0.5_rt  * (
                     fine(ii-1, jj  , kk  ) + fine(ii+1, jj  , kk  ) +
                     fine(ii  , jj  , kk-1) + fine(ii  , jj  , kk+1)
-                 ) +
-                 0.25 * (
+                ) +
+                0.25_rt * (
                     fine(ii-1, jj  , kk-1) + fine(ii+1, jj  , kk-1) +
                     fine(ii-1, jj  , kk+1) + fine(ii+1, jj  , kk+1)
-                 ) +
-                 fine(ii, jj+1, kk) +
-                 0.5  * (
+                ) +
+                fine(ii, jj+1, kk) +
+                0.5_rt  * (
                     fine(ii-1, jj+1, kk  ) + fine(ii+1, jj+1, kk  ) +
                     fine(ii  , jj+1, kk-1) + fine(ii  , jj+1, kk+1)
-                 ) +
-                 0.25 * (
+                ) +
+                0.25_rt * (
                     fine(ii-1, jj+1, kk-1) + fine(ii+1, jj+1, kk-1) +
                     fine(ii-1, jj+1, kk+1) + fine(ii+1, jj+1, kk+1)
-                 )
+                )
              );
         } else {
-            coarse(i, j, k) = 0.125 * (
+            coarse(i, j, k) = 0.125_rt * (
                 fine(ii, jj, kk  ) +
-                0.5  * (
+                0.5_rt  * (
                     fine(ii-1, jj  , kk  ) + fine(ii+1, jj  , kk  ) +
                     fine(ii  , jj-1, kk  ) + fine(ii  , jj+1, kk  )
                 ) +
-                0.25 * (
+                0.25_rt * (
                     fine(ii-1, jj-1, kk  ) + fine(ii+1, jj-1, kk  ) +
                     fine(ii-1, jj+1, kk  ) + fine(ii+1, jj+1, kk  )
                 ) +
                 fine(ii, jj, kk+1) +
-                0.5  * (
+                0.5_rt  * (
                     fine(ii-1, jj  , kk+1) + fine(ii+1, jj  , kk+1) +
                     fine(ii  , jj-1, kk+1) + fine(ii  , jj+1, kk+1)
                 ) +
-                0.25 * (
+                0.25_rt * (
                     fine(ii-1, jj-1, kk+1) + fine(ii+1, jj-1, kk+1) +
                     fine(ii-1, jj+1, kk+1) + fine(ii+1, jj+1, kk+1)
                 )
diff --git a/Source/Parallelization/InterpolateDensityFineToCoarse.H b/Source/Parallelization/InterpolateDensityFineToCoarse.H
new file mode 100644
index 000000000..947db2aac
--- /dev/null
+++ b/Source/Parallelization/InterpolateDensityFineToCoarse.H
@@ -0,0 +1,120 @@
+/* Copyright 2019 Axel Huebl, Weiqun Zhang
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef INTERPOLATE_DENSITY_FINE_TO_COARSE_H
+#define INTERPOLATE_DENSITY_FINE_TO_COARSE_H
+
+#include <AMReX_Array4.H>
+#include <AMReX_REAL.H>
+#include <AMReX_BLassert.H>
+#include <AMReX_Extension.H>
+#include <AMReX_GpuQualifiers.H>
+
+#include <utility> // std::move
+
+
+/** Fill a charge density (rho) coarse patch with averaged values from a fine patch
+ *
+ * Fills the values of the charge density on the coarse patch
+ * by averaging the values of the charge density of the fine patch.
+ */
+class InterpolateDensityFineToCoarse
+{
+public:
+    /** Construct with fine and coarse patch and their refinement ratio
+     *
+     * @param[in]     fine read-only fine patch
+     * @param[in,out] coarse overwritten coarse patch
+     * @param[in]     refinement_ratio ratio between coarse and fine patch granularity
+     *                (currently, only a value of is implemented)
+     * @param[in]     number_of_components the number of components
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    InterpolateDensityFineToCoarse(
+        amrex::Array4<amrex::Real const> const fine,
+        amrex::Array4<amrex::Real      > const coarse,
+        int const refinement_ratio,
+        int const number_of_components
+    ) : m_fine(std::move(fine)),
+        m_coarse(std::move(coarse)),
+        m_refinement_ratio(std::move(refinement_ratio)),
+        m_number_of_components(std::move(number_of_components))
+    {
+        //! @note constants and stencils in operator() implementation assume 2x refinement
+        BL_ASSERT(refinement_ratio == 2);
+    }
+
+    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+    void
+    operator()(
+        int const i,
+        int const j,
+        int const k
+    ) const noexcept
+    {
+        using namespace amrex;
+
+        auto const & fine_unsafe = m_fine;   // out-of-bounds access not secured with zero-values yet
+        auto const & coarse      = m_coarse; // out-of-bounds access not secured but will also not occur
+
+        // return zero for out-of-bounds accesses during interpolation
+        // this is efficiently used as a method to add neutral elements beyond guards in the average below
+        auto const fine = [fine_unsafe] AMREX_GPU_DEVICE (int const j, int const k, int const l, int const m) noexcept
+        {
+            return fine_unsafe.contains(j, k, l) ? fine_unsafe(j, k, l, m) :  amrex::Real{0.};
+        };
+
+        int const ii = i * m_refinement_ratio;
+        int const jj = j * m_refinement_ratio;
+        int const kk = k * m_refinement_ratio;
+        for( int m = 0; m < m_number_of_components; ++m ) {
+#if AMREX_SPACEDIM == 2
+            coarse(i,j,kk,m) = 0.25_rt * (
+                fine(ii,jj,kk,m) +
+                0.5_rt  * (
+                    fine(ii-1,jj  ,kk,m) + fine(ii+1,jj  ,kk,m) +
+                    fine(ii  ,jj-1,kk,m) + fine(ii  ,jj+1,kk,m)
+                ) +
+                0.25_rt * (
+                    fine(ii-1,jj-1,kk,m) + fine(ii+1,jj-1,kk,m) +
+                    fine(ii-1,jj+1,kk,m) + fine(ii+1,jj+1,kk,m)
+                )
+            );
+#elif AMREX_SPACEDIM == 3
+            coarse(i,j,k,m) = 0.125_rt * (
+                fine(ii,jj,kk,m) +
+                0.5_rt   * (
+                    fine(ii-1,jj  ,kk  ,m) + fine(ii+1,jj  ,kk  ,m) +
+                    fine(ii  ,jj-1,kk  ,m) + fine(ii  ,jj+1,kk  ,m) +
+                    fine(ii  ,jj  ,kk-1,m) + fine(ii  ,jj  ,kk+1,m)
+                ) +
+                0.25_rt  * (
+                    fine(ii-1,jj-1,kk  ,m) + fine(ii+1,jj-1,kk  ,m) +
+                    fine(ii-1,jj+1,kk  ,m) + fine(ii+1,jj+1,kk  ,m) +
+                    fine(ii-1,jj  ,kk-1,m) + fine(ii+1,jj  ,kk-1,m) +
+                    fine(ii-1,jj  ,kk+1,m) + fine(ii+1,jj  ,kk+1,m) +
+                    fine(ii  ,jj-1,kk-1,m) + fine(ii  ,jj+1,kk-1,m) +
+                    fine(ii  ,jj-1,kk+1,m) + fine(ii  ,jj+1,kk+1,m)
+                ) +
+                0.125_rt * (
+                    fine(ii-1,jj-1,kk-1,m) + fine(ii-1,jj-1,kk+1,m) +
+                    fine(ii-1,jj+1,kk-1,m) + fine(ii-1,jj+1,kk+1,m) +
+                    fine(ii+1,jj-1,kk-1,m) + fine(ii+1,jj-1,kk+1,m) +
+                    fine(ii+1,jj+1,kk-1,m) + fine(ii+1,jj+1,kk+1,m)
+                )
+            );
+#endif
+        }
+    }
+private:
+    amrex::Array4< amrex::Real const > m_fine;
+    amrex::Array4< amrex::Real       > m_coarse;
+    int m_refinement_ratio;
+    int m_number_of_components;
+};
+
+#endif // INTERPOLATE_DENSITY_FINE_TO_COARSE_H
diff --git a/Source/Parallelization/Make.package b/Source/Parallelization/Make.package
index c74583522..7c3c38627 100644
--- a/Source/Parallelization/Make.package
+++ b/Source/Parallelization/Make.package
@@ -2,6 +2,7 @@ CEXE_sources += WarpXComm.cpp
 CEXE_sources += WarpXRegrid.cpp
 CEXE_headers += WarpXSumGuardCells.H
 CEXE_headers += WarpXComm_K.H
+CEXE_headers += WarpXComm.H
 
 INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Parallelization
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Parallelization
diff --git a/Source/Parallelization/WarpXComm.H b/Source/Parallelization/WarpXComm.H
new file mode 100644
index 000000000..81f00088e
--- /dev/null
+++ b/Source/Parallelization/WarpXComm.H
@@ -0,0 +1,33 @@
+#ifndef WARPX_PARALLELIZATION_COMM_H_
+#define WARPX_PARALLELIZATION_COMM_H_
+
+#include <AMReX_MultiFab.H>
+
+/** \brief Fills the values of the current on the coarse patch by
+ *  averaging the values of the current of the fine patch (on the same level).
+ *  Also fills the guards of the coarse patch.
+ *
+ * \param[in]  fine             fine patches to interpolate from
+ * \param[out] coarse           coarse patches to interpolate to
+ * \param[in]  refinement_ratio integer ratio between the two
+ */
+void
+interpolateCurrentFineToCoarse (
+      std::array< amrex::MultiFab const *, 3 > const & fine,
+      std::array< amrex::MultiFab       *, 3 > const & coarse,
+      int const refinement_ratio);
+
+/** \brief Fills the values of the charge density on the coarse patch by
+ *  averaging the values of the charge density of the fine patch (on the same level).
+ *
+ * \param[in]  fine             fine patches to interpolate from
+ * \param[out] coarse           coarse patches to interpolate to
+ * \param[in]  refinement_ratio integer ratio between the two
+ */
+void
+interpolateDensityFineToCoarse (
+     const amrex::MultiFab& fine,
+     amrex::MultiFab& coarse,
+     int const refinement_ratio);
+
+#endif // WARPX_PARALLELIZATION_COMM_H_
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index 4f870e79c..b61ae4fc7 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -1,8 +1,10 @@
+#include <WarpXComm.H>
 #include <WarpXComm_K.H>
 #include <WarpX.H>
 #include <WarpX_f.H>
 #include <WarpXSumGuardCells.H>
-#include <Parallelization/InterpolateCurrentFineToCoarse.H>
+#include <InterpolateCurrentFineToCoarse.H>
+#include <InterpolateDensityFineToCoarse.H>
 
 #include <algorithm>
 #include <cstdlib>
@@ -52,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();
@@ -350,11 +503,11 @@ WarpX::SyncCurrent ()
 }
 
 void
-WarpX::interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 > const & fine,
-                                        std::array< amrex::MultiFab       *, 3 > const & coarse,
-                                        int const refinement_ratio)
+interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 > const & fine,
+                                 std::array< amrex::MultiFab       *, 3 > const & coarse,
+                                 int const refinement_ratio)
 {
-    BL_PROFILE("InterpolateCurrentFineToCoarse()");
+    BL_PROFILE("interpolateCurrentFineToCoarse()");
     BL_ASSERT(refinement_ratio == 2);
     const IntVect& ng = (fine[0]->nGrowVect() + 1) / refinement_ratio; // add equivalent no. of guards to coarse patch
 
@@ -386,6 +539,8 @@ WarpX::interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 >
 void
 WarpX::SyncRho ()
 {
+    BL_PROFILE("SyncRho()");
+
     if (!rho_fp[0]) return;
     const int ncomp = rho_fp[0]->nComp();
 
@@ -395,7 +550,7 @@ WarpX::SyncRho ()
     {
         rho_cp[lev]->setVal(0.0);
         const IntVect& refinement_ratio = refRatio(lev-1);
-        SyncRho(*rho_fp[lev], *rho_cp[lev], refinement_ratio[0]);
+        interpolateDensityFineToCoarse(*rho_fp[lev], *rho_cp[lev], refinement_ratio[0]);
     }
 
     // For each level
@@ -408,29 +563,26 @@ WarpX::SyncRho ()
 }
 
 void
-WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int refinement_ratio)
+interpolateDensityFineToCoarse (const MultiFab& fine, MultiFab& coarse, int const refinement_ratio)
 {
+    BL_PROFILE("interpolateDensityFineToCoarse()");
     BL_ASSERT(refinement_ratio == 2);
-    const IntVect& ng = (fine.nGrowVect()+1)/refinement_ratio;
+    const IntVect& ng = (fine.nGrowVect() + 1) / refinement_ratio;  // add equivalent no. of guards to coarse patch
     const int nc = fine.nComp();
 
 #ifdef _OPEMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
     {
-        FArrayBox ffab;
-        for (MFIter mfi(crse,true); mfi.isValid(); ++mfi)
+        // OMP in-box decomposition of coarse into tilebox
+        for (MFIter mfi(coarse, TilingIfNotGPU()); mfi.isValid(); ++mfi)
         {
-            const Box& bx = mfi.growntilebox(ng);
-            Box fbx = amrex::grow(amrex::refine(bx,refinement_ratio),1);
-            ffab.resize(fbx, nc);
-            fbx &= fine[mfi].box();
-            ffab.setVal(0.0);
-            ffab.copy(fine[mfi], fbx, 0, fbx, 0, nc);
-            WRPX_SYNC_RHO(bx.loVect(), bx.hiVect(),
-                          BL_TO_FORTRAN_ANYD(crse[mfi]),
-                          BL_TO_FORTRAN_ANYD(ffab),
-                          &nc);
+            const Box& bx = mfi.growntilebox(ng); // only grow to outer directions of tileboxes for filling guards
+
+            amrex::ParallelFor(
+                bx,
+                InterpolateDensityFineToCoarse(fine.const_array(mfi), coarse.array(mfi), refinement_ratio, nc)
+            );
         }
     }
 }
@@ -562,7 +714,7 @@ WarpX::RestrictRhoFromFineToCoarsePatch (int lev)
     if (rho_fp[lev]) {
         rho_cp[lev]->setVal(0.0);
         const IntVect& refinement_ratio = refRatio(lev-1);
-        SyncRho(*rho_fp[lev], *rho_cp[lev], refinement_ratio[0]);
+        interpolateDensityFineToCoarse(*rho_fp[lev], *rho_cp[lev], refinement_ratio[0]);
     }
 }
 
diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H
index 093323ec3..169cd0ee1 100644
--- a/Source/Parallelization/WarpXComm_K.H
+++ b/Source/Parallelization/WarpXComm_K.H
@@ -5,38 +5,38 @@
 
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
 void warpx_interp_bfield_x (int j, int k, int l,
-                            amrex::Array4<amrex::Real> const& Bxa,
+                            amrex::Array4<amrex::Real      > const& Bxa,
                             amrex::Array4<amrex::Real const> const& Bxf,
                             amrex::Array4<amrex::Real const> const& Bxc)
 {
     using namespace amrex;
 
-    int lg = amrex::coarsen(l,2);
-    int kg = amrex::coarsen(k,2);
-    int jg = amrex::coarsen(j,2);
+    int const lg = amrex::coarsen(l,2);
+    int const kg = amrex::coarsen(k,2);
+    int const jg = amrex::coarsen(j,2);
 
-    Real wx = (j == jg*2) ? 0.0 : 0.5;
-    Real owx = 1.0-wx;
+    Real const wx = (j == jg*2) ? 0.0_rt : 0.5_rt;
+    Real const owx = 1.0_rt - wx;
     Bxa(j,k,l) = owx * Bxc(jg,kg,lg) + wx * Bxc(jg+1,kg,lg) + Bxf(j,k,l);
 }
 
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
 void warpx_interp_bfield_y (int j, int k, int l,
-                            amrex::Array4<amrex::Real> const& Bya,
+                            amrex::Array4<amrex::Real      > const& Bya,
                             amrex::Array4<amrex::Real const> const& Byf,
                             amrex::Array4<amrex::Real const> const& Byc)
 {
     using namespace amrex;
 
-    int lg = amrex::coarsen(l,2);
-    int kg = amrex::coarsen(k,2);
-    int jg = amrex::coarsen(j,2);
+    int const lg = amrex::coarsen(l,2);
+    int const kg = amrex::coarsen(k,2);
+    int const jg = amrex::coarsen(j,2);
 
     // Note that for 2d, l=0, because the amrex convention is used here.
 
 #if (AMREX_SPACEDIM == 3)
-    Real wy = (k == kg*2) ? 0.0 : 0.5;
-    Real owy = 1.0-wy;
+    Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt;
+    Real const owy = 1.0_rt - wy;
     Bya(j,k,l) = owy * Byc(jg,kg,lg) + wy * Byc(jg,kg+1,lg) + Byf(j,k,l);
 #else
     Bya(j,k,l) = Byc(jg,kg,lg) + Byf(j,k,l);
@@ -45,47 +45,47 @@ void warpx_interp_bfield_y (int j, int k, int l,
 
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
 void warpx_interp_bfield_z (int j, int k, int l,
-                            amrex::Array4<amrex::Real> const& Bza,
+                            amrex::Array4<amrex::Real      > const& Bza,
                             amrex::Array4<amrex::Real const> const& Bzf,
                             amrex::Array4<amrex::Real const> const& Bzc)
 {
     using namespace amrex;
 
-    int lg = amrex::coarsen(l,2);
-    int kg = amrex::coarsen(k,2);
-    int jg = amrex::coarsen(j,2);
+    int const lg = amrex::coarsen(l,2);
+    int const kg = amrex::coarsen(k,2);
+    int const jg = amrex::coarsen(j,2);
 
     // Note that for 2d, l=0, because the amrex convention is used here.
 
 #if (AMREX_SPACEDIM == 3)
-    Real wz = (l == lg*2) ? 0.0 : 0.5;
-    Real owz = 1.0-wz;
+    Real const wz = (l == lg*2) ? 0.0_rt : 0.5_rt;
+    Real const owz = 1.0_rt - wz;
     Bza(j,k,l) = owz * Bzc(jg,kg,lg) + owz * Bzc(jg,kg,lg+1) + Bzf(j,k,l);
 #else
-    Real wy = (k == kg*2) ? 0.0 : 0.5;
-    Real owy = 1.0-wy;
+    Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt;
+    Real const owy = 1.0_rt - wy;
     Bza(j,k,l) = owy * Bzc(jg,kg,lg) + owy * Bzc(jg,kg+1,lg) + Bzf(j,k,l);
 #endif
 }
 
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
 void warpx_interp_efield_x (int j, int k, int l,
-                            amrex::Array4<amrex::Real> const& Exa,
+                            amrex::Array4<amrex::Real      > const& Exa,
                             amrex::Array4<amrex::Real const> const& Exf,
                             amrex::Array4<amrex::Real const> const& Exc)
 {
     using namespace amrex;
 
-    int lg = amrex::coarsen(l,2);
-    int kg = amrex::coarsen(k,2);
-    int jg = amrex::coarsen(j,2);
+    int const lg = amrex::coarsen(l,2);
+    int const kg = amrex::coarsen(k,2);
+    int const jg = amrex::coarsen(j,2);
 
-    Real wy = (k == kg*2) ? 0.0 : 0.5;
-    Real owy = 1.0-wy;
+    Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt;
+    Real const owy = 1.0_rt - wy;
 
 #if (AMREX_SPACEDIM == 3)
-    Real wz = (l == lg*2) ? 0.0 : 0.5;
-    Real owz = 1.0-wz;
+    Real const wz = (l == lg*2) ? 0.0_rt : 0.5_rt;
+    Real const owz = 1.0_rt - wz;
     Exa(j,k,l) = owy * owz * Exc(jg  ,kg  ,lg  )
         +         wy * owz * Exc(jg  ,kg+1,lg  )
         +        owy *  wz * Exc(jg  ,kg  ,lg+1)
@@ -98,30 +98,30 @@ void warpx_interp_efield_x (int j, int k, int l,
 
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
 void warpx_interp_efield_y (int j, int k, int l,
-                            amrex::Array4<amrex::Real> const& Eya,
+                            amrex::Array4<amrex::Real      > const& Eya,
                             amrex::Array4<amrex::Real const> const& Eyf,
                             amrex::Array4<amrex::Real const> const& Eyc)
 {
     using namespace amrex;
 
-    int lg = amrex::coarsen(l,2);
-    int kg = amrex::coarsen(k,2);
-    int jg = amrex::coarsen(j,2);
+    int const lg = amrex::coarsen(l,2);
+    int const kg = amrex::coarsen(k,2);
+    int const jg = amrex::coarsen(j,2);
 
-    Real wx = (j == jg*2) ? 0.0 : 0.5;
-    Real owx = 1.0-wx;
+    Real const wx = (j == jg*2) ? 0.0_rt : 0.5_rt;
+    Real const owx = 1.0_rt - wx;
 
 #if (AMREX_SPACEDIM == 3)
-    Real wz = (l == lg*2) ? 0.0 : 0.5;
-    Real owz = 1.0-wz;
+    Real const wz = (l == lg*2) ? 0.0_rt : 0.5_rt;
+    Real const owz = 1.0_rt - wz;
     Eya(j,k,l) = owx * owz * Eyc(jg  ,kg  ,lg  )
         +         wx * owz * Eyc(jg+1,kg  ,lg  )
         +        owx *  wz * Eyc(jg  ,kg  ,lg+1)
         +         wx *  wz * Eyc(jg+1,kg  ,lg+1)
         + Eyf(j,k,l);
 #else
-    Real wy = (k == kg*2) ? 0.0 : 0.5;
-    Real owy = 1.0-wy;
+    Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt;
+    Real const owy = 1.0_rt - wy;
     Eya(j,k,l) = owx * owy * Eyc(jg  ,kg  ,lg)
         +         wx * owy * Eyc(jg+1,kg  ,lg)
         +        owx *  wy * Eyc(jg  ,kg+1,lg)
@@ -132,22 +132,22 @@ void warpx_interp_efield_y (int j, int k, int l,
 
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
 void warpx_interp_efield_z (int j, int k, int l,
-                            amrex::Array4<amrex::Real> const& Eza,
+                            amrex::Array4<amrex::Real      > const& Eza,
                             amrex::Array4<amrex::Real const> const& Ezf,
                             amrex::Array4<amrex::Real const> const& Ezc)
 {
     using namespace amrex;
 
-    int lg = amrex::coarsen(l,2);
-    int kg = amrex::coarsen(k,2);
-    int jg = amrex::coarsen(j,2);
+    int const lg = amrex::coarsen(l,2);
+    int const kg = amrex::coarsen(k,2);
+    int const jg = amrex::coarsen(j,2);
 
-    Real wx = (j == jg*2) ? 0.0 : 0.5;
-    Real owx = 1.0-wx;
+    Real const wx = (j == jg*2) ? 0.0_rt : 0.5_rt;
+    Real const owx = 1.0_rt - wx;
 
 #if (AMREX_SPACEDIM == 3)
-    Real wy = (k == kg*2) ? 0.0 : 0.5;
-    Real owy = 1.0-wy;
+    Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt;
+    Real owy = 1.0_rt - wy;
     Eza(j,k,l) = owx * owy * Ezc(jg  ,kg  ,lg  )
         +         wx * owy * Ezc(jg+1,kg  ,lg  )
         +        owx *  wy * Ezc(jg  ,kg+1,lg  )
@@ -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()));