Merge pull request #225 from ECP-WarpX/EBpush_to_cpp

EB push to cpp

3 files changed, 515 insertions, 281 deletions

diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index bea008598..4fce4717b 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -6,6 +6,7 @@
 #include <WarpXConst.H>
 #include <WarpX_f.H>
 #include <WarpX_K.H>
+#include <WarpX_FDTD.H>
 #ifdef WARPX_USE_PY
 #include <WarpX_py.H>
 #endif
@@ -88,7 +89,8 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt)
 {
     const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
     const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
-    Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2];
+    const Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2];
+    const Real dxinv = 1./dx[0];
 
     MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz;
     if (patch_type == PatchType::fine)
@@ -129,13 +131,13 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt)
         const Box& tby  = mfi.tilebox(By_nodal_flag);
         const Box& tbz  = mfi.tilebox(Bz_nodal_flag);
 
+        auto const& Bxfab = Bx->array(mfi);
+        auto const& Byfab = By->array(mfi);
+        auto const& Bzfab = Bz->array(mfi);
+        auto const& Exfab = Ex->array(mfi);
+        auto const& Eyfab = Ey->array(mfi);
+        auto const& Ezfab = Ez->array(mfi);
         if (do_nodal) {
-            auto const& Bxfab = Bx->array(mfi);
-            auto const& Byfab = By->array(mfi);
-            auto const& Bzfab = Bz->array(mfi);
-            auto const& Exfab = Ex->array(mfi);
-            auto const& Eyfab = Ey->array(mfi);
-            auto const& Ezfab = Ez->array(mfi);
             amrex::ParallelFor(tbx,
             [=] AMREX_GPU_DEVICE (int j, int k, int l)
             {
@@ -151,21 +153,54 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt)
             {
                 warpx_push_bz_nodal(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy);
             });
-        } else {
-            // 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, &dtsdy, &dtsdz,
-                      &xmin, &dx[0],
-		      &WarpX::maxwell_fdtd_solver_id);
+        } else if (WarpX::maxwell_fdtd_solver_id == 0) {
+            amrex::ParallelFor(tbx,
+            [=] AMREX_GPU_DEVICE (int j, int k, int l)
+            {
+                warpx_push_bx_yee(j,k,l,Bxfab,Eyfab,Ezfab,dtsdy,dtsdz);
+            });
+            amrex::ParallelFor(tby,
+            [=] AMREX_GPU_DEVICE (int j, int k, int l)
+            {
+                warpx_push_by_yee(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz);
+            });
+            amrex::ParallelFor(tbz,
+            [=] AMREX_GPU_DEVICE (int j, int k, int l)
+            {
+                warpx_push_bz_yee(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy,dxinv,xmin);
+            });
+        } else if (WarpX::maxwell_fdtd_solver_id == 1) {
+            Real betaxy, betaxz, betayx, betayz, betazx, betazy;
+            Real gammax, gammay, gammaz;
+            Real alphax, alphay, alphaz;
+            warpx_calculate_ckc_coefficients(dtsdx, dtsdy, dtsdz,
+                                             betaxy, betaxz, betayx, betayz, betazx, betazy,
+                                             gammax, gammay, gammaz,
+                                             alphax, alphay, alphaz);
+            amrex::ParallelFor(tbx,
+            [=] AMREX_GPU_DEVICE (int j, int k, int l)
+            {
+                warpx_push_bx_ckc(j,k,l,Bxfab,Eyfab,Ezfab,
+                                  betaxy, betaxz, betayx, betayz, betazx, betazy,
+                                  gammax, gammay, gammaz,
+                                  alphax, alphay, alphaz);
+            });
+            amrex::ParallelFor(tby,
+            [=] AMREX_GPU_DEVICE (int j, int k, int l)
+            {
+                warpx_push_by_ckc(j,k,l,Byfab,Exfab,Ezfab,
+                                  betaxy, betaxz, betayx, betayz, betazx, betazy,
+                                  gammax, gammay, gammaz,
+                                  alphax, alphay, alphaz);
+            });
+            amrex::ParallelFor(tbz,
+            [=] AMREX_GPU_DEVICE (int j, int k, int l)
+            {
+                warpx_push_bz_ckc(j,k,l,Bzfab,Exfab,Eyfab,
+                                  betaxy, betaxz, betayx, betayz, betazx, betazy,
+                                  gammax, gammay, gammaz,
+                                  alphax, alphay, alphaz);
+            });
         }
 
         if (cost) {
@@ -237,9 +272,10 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt)
     const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * a_dt;
     const Real c2dt = (PhysConst::c*PhysConst::c) * a_dt;
 
-    int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
+    const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
     const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
-    Real dtsdx_c2 = c2dt/dx[0], dtsdy_c2 = c2dt/dx[1], dtsdz_c2 = c2dt/dx[2];
+    const Real dtsdx_c2 = c2dt/dx[0], dtsdy_c2 = c2dt/dx[1], dtsdz_c2 = c2dt/dx[2];
+    const Real dxinv = 1./dx[0];
 
     MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F;
     if (patch_type == PatchType::fine)
@@ -288,16 +324,17 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt)
         const Box& tey  = mfi.tilebox(Ey_nodal_flag);
         const Box& tez  = mfi.tilebox(Ez_nodal_flag);
 
+        auto const& Exfab = Ex->array(mfi);
+        auto const& Eyfab = Ey->array(mfi);
+        auto const& Ezfab = Ez->array(mfi);
+        auto const& Bxfab = Bx->array(mfi);
+        auto const& Byfab = By->array(mfi);
+        auto const& Bzfab = Bz->array(mfi);
+        auto const& jxfab = jx->array(mfi);
+        auto const& jyfab = jy->array(mfi);
+        auto const& jzfab = jz->array(mfi);
+
         if (do_nodal) {
-            auto const& Exfab = Ex->array(mfi);
-            auto const& Eyfab = Ey->array(mfi);
-            auto const& Ezfab = Ez->array(mfi);
-            auto const& Bxfab = Bx->array(mfi);
-            auto const& Byfab = By->array(mfi);
-            auto const& Bzfab = Bz->array(mfi);
-            auto const& jxfab = jx->array(mfi);
-            auto const& jyfab = jy->array(mfi);
-            auto const& jzfab = jz->array(mfi);
             amrex::ParallelFor(tex,
             [=] AMREX_GPU_DEVICE (int j, int k, int l)
             {
@@ -314,37 +351,76 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt)
                 warpx_push_ez_nodal(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2);
             });
         } else {
-            // 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, &dtsdy_c2, &dtsdz_c2,
-		      &xmin, &dx[0]);
+            amrex::ParallelFor(tex,
+            [=] AMREX_GPU_DEVICE (int j, int k, int l)
+            {
+                warpx_push_ex_yee(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdy_c2,dtsdz_c2);
+            });
+            amrex::ParallelFor(tey,
+            [=] AMREX_GPU_DEVICE (int j, int k, int l)
+            {
+                warpx_push_ey_yee(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,mu_c2_dt,dtsdx_c2,dtsdz_c2,xmin);
+            });
+            amrex::ParallelFor(tez,
+            [=] AMREX_GPU_DEVICE (int j, int k, int l)
+            {
+                warpx_push_ez_yee(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2,dxinv,xmin);
+            });
         }
 
         if (F)
         {
-            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, &dtsdy_c2, &dtsdz_c2,
-			  &WarpX::maxwell_fdtd_solver_id);
+            auto const& Ffab = F->array(mfi);
+            if (WarpX::maxwell_fdtd_solver_id == 0) {
+                amrex::ParallelFor(tex,
+                [=] AMREX_GPU_DEVICE (int j, int k, int l)
+                {
+                    warpx_push_ex_f_yee(j,k,l,Exfab,Ffab,dtsdx_c2);
+                });
+                amrex::ParallelFor(tey,
+                [=] AMREX_GPU_DEVICE (int j, int k, int l)
+                {
+                    warpx_push_ey_f_yee(j,k,l,Eyfab,Ffab,dtsdy_c2);
+                });
+                amrex::ParallelFor(tez,
+                [=] AMREX_GPU_DEVICE (int j, int k, int l)
+                {
+                    warpx_push_ez_f_yee(j,k,l,Ezfab,Ffab,dtsdz_c2);
+                });
+            }
+            else if (WarpX::maxwell_fdtd_solver_id == 1) {
+                Real betaxy, betaxz, betayx, betayz, betazx, betazy;
+                Real gammax, gammay, gammaz;
+                Real alphax, alphay, alphaz;
+                warpx_calculate_ckc_coefficients(dtsdx_c2, dtsdy_c2, dtsdz_c2,
+                                                 betaxy, betaxz, betayx, betayz, betazx, betazy,
+                                                 gammax, gammay, gammaz,
+                                                 alphax, alphay, alphaz);
+                amrex::ParallelFor(tex,
+                [=] AMREX_GPU_DEVICE (int j, int k, int l)
+                {
+                    warpx_push_ex_f_ckc(j,k,l,Exfab,Ffab,
+                                        betaxy, betaxz, betayx, betayz, betazx, betazy,
+                                        gammax, gammay, gammaz,
+                                        alphax, alphay, alphaz);
+                });
+                amrex::ParallelFor(tey,
+                [=] AMREX_GPU_DEVICE (int j, int k, int l)
+                {
+                    warpx_push_ey_f_ckc(j,k,l,Eyfab,Ffab,
+                                        betaxy, betaxz, betayx, betayz, betazx, betazy,
+                                        gammax, gammay, gammaz,
+                                        alphax, alphay, alphaz);
+                });
+                amrex::ParallelFor(tez,
+                [=] AMREX_GPU_DEVICE (int j, int k, int l)
+                {
+                    warpx_push_ez_f_ckc(j,k,l,Ezfab,Ffab,
+                                        betaxy, betaxz, betayx, betayz, betazx, betazy,
+                                        gammax, gammay, gammaz,
+                                        alphax, alphay, alphaz);
+                });
+            }
         }
 
         if (cost) {
@@ -434,7 +510,7 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type)
 
     static constexpr Real mu_c2 = PhysConst::mu0*PhysConst::c*PhysConst::c;
 
-    int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
+    const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
     const auto& dx = WarpX::CellSize(patch_level);
     const std::array<Real,3> dtsdx {a_dt/dx[0], a_dt/dx[1], a_dt/dx[2]};
 
diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H
index 865277446..e2d40f4c2 100644
--- a/Source/FieldSolver/WarpX_FDTD.H
+++ b/Source/FieldSolver/WarpX_FDTD.H
@@ -5,50 +5,411 @@
 
 using namespace amrex;
 
-AMREX_GPU_HOST_DEVICE AMREX_INLINE
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_bx_yee(int i, int j, int k, Array4<Real> const& Bx,
+                       Array4<Real const> const& Ey, Array4<Real const> const& Ez,
+                       Real dtsdy, Real dtsdz)
+{
+#if defined WARPX_DIM_3D
+    Bx(i,j,k) += - dtsdy * (Ez(i,j+1,k  ) - Ez(i,j,k))
+                 + dtsdz * (Ey(i,j  ,k+1) - Ey(i,j,k));
+#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ)
+    // Note that the 2D Cartesian and RZ versions are the same
+    Bx(i,j,0) += + dtsdz * (Ey(i,j+1,0) - Ey(i,j,0));
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_by_yee(int i, int j, int k, Array4<Real> const& By,
+                       Array4<Real const> const& Ex, Array4<Real const> const& Ez,
+                       Real dtsdx, Real dtsdz)
+{
+#if defined WARPX_DIM_3D
+    By(i,j,k) += + dtsdx * (Ez(i+1,j,k  ) - Ez(i,j,k))
+                 - dtsdz * (Ex(i  ,j,k+1) - Ex(i,j,k));
+#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ)
+    // Note that the 2D Cartesian and RZ versions are the same
+    By(i,j,0) += + dtsdx * (Ez(i+1,j  ,0) - Ez(i,j,0))
+                 - dtsdz * (Ex(i  ,j+1,0) - Ex(i,j,0));
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_bz_yee(int i, int j, int k, Array4<Real> const& Bz,
+                       Array4<Real const> const& Ex, Array4<Real const> const& Ey,
+                       Real dtsdx, Real dtsdy, Real dxinv, Real rmin)
+{
+#if defined WARPX_DIM_3D
+    Bz(i,j,k) += - dtsdx * (Ey(i+1,j  ,k) - Ey(i,j,k))
+                 + dtsdy * (Ex(i  ,j+1,k) - Ex(i,j,k));
+#elif defined WARPX_DIM_2D
+    Bz(i,j,0) += - dtsdx * (Ey(i+1,j,0) - Ey(i,j,0));
+#elif defined WARPX_DIM_RZ
+    const Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5);
+    const Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5);
+    Bz(i,j,0) += - dtsdx * (ru*Ey(i+1,j,0) - rd*Ey(i,j,0));
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_ex_yee(int i, int j, int k, Array4<Real> const& Ex,
+                       Array4<Real const> const& By, Array4<Real const> const& Bz, Array4<Real const> const& Jx,
+                       Real mu_c2_dt, Real dtsdy_c2, Real dtsdz_c2)
+{
+#if defined WARPX_DIM_3D
+    Ex(i,j,k) += + dtsdy_c2 * (Bz(i,j,k) - Bz(i,j-1,k  ))
+                 - dtsdz_c2 * (By(i,j,k) - By(i,j  ,k-1))
+                 - mu_c2_dt  * Jx(i,j,k);
+#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ)
+    // Note that the 2D Cartesian and RZ versions are the same
+    Ex(i,j,0) += - dtsdz_c2 * (By(i,j,0) - By(i,j-1,0))
+                 - mu_c2_dt  * Jx(i,j,0);
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_ey_yee(int i, int j, int k, Array4<Real> const& Ey,
+                       Array4<Real const> const& Bx, Array4<Real const> const& Bz, Array4<Real const> const& Jy,
+                       Real mu_c2_dt, Real dtsdx_c2, Real dtsdz_c2, Real rmin)
+{
+#if defined WARPX_DIM_3D
+    Ey(i,j,k) += - dtsdx_c2 * (Bz(i,j,k) - Bz(i-1,j,k))
+                 + dtsdz_c2 * (Bx(i,j,k) - Bx(i,j,k-1))
+                 - mu_c2_dt  * Jy(i,j,k);
+#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ)
+    // 2D Cartesian and RZ are the same, except that the axis is skipped with RZ
+#ifdef WARPX_DIM_RZ
+    if (i != 0 || rmin != 0.)
+#endif
+    {
+    Ey(i,j,0) += - dtsdx_c2 * (Bz(i,j,0) - Bz(i-1,j,0))
+                 + dtsdz_c2 * (Bx(i,j,0) - Bx(i,j-1,0))
+                 - mu_c2_dt  * Jy(i,j,0);
+    }
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_ez_yee(int i, int j, int k, Array4<Real> const& Ez,
+                       Array4<Real const> const& Bx, Array4<Real const> const& By, Array4<Real const> const& Jz,
+                       Real mu_c2_dt, Real dtsdx_c2, Real dtsdy_c2, Real dxinv, Real rmin)
+{
+#if defined WARPX_DIM_3D
+    Ez(i,j,k) += + dtsdx_c2 * (By(i,j,k) - By(i-1,j  ,k))
+                 - dtsdy_c2 * (Bx(i,j,k) - Bx(i  ,j-1,k))
+                 - mu_c2_dt  * Jz(i,j,k);
+#elif defined WARPX_DIM_2D
+    Ez(i,j,0) += + dtsdx_c2 * (By(i,j,0) - By(i-1,j,0))
+                 - mu_c2_dt  * Jz(i,j,0);
+#elif defined WARPX_DIM_RZ
+    if (i != 0 || rmin != 0.) {
+        const Real ru = 1. + 0.5/(rmin*dxinv + i);
+        const Real rd = 1. - 0.5/(rmin*dxinv + i);
+        Ez(i,j,0) += + dtsdx_c2 * (ru*By(i,j,0) - rd*By(i-1,j,0))
+                     - mu_c2_dt  * Jz(i,j,0);
+    } else {
+        Ez(i,j,0) += + 4.*dtsdx_c2 * By(i,j,0)
+                       - mu_c2_dt  * Jz(i,j,0);
+    }
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_ex_f_yee(int j, int k, int l, Array4<Real> const& Ex,
+                         Array4<Real const> const& F, Real dtsdx_c2)
+{
+#if defined WARPX_DIM_3D
+    Ex(j,k,l) += + dtsdx_c2 * (F(j+1,k,l) - F(j,k,l));
+#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ)
+    Ex(j,k,0) += + dtsdx_c2 * (F(j+1,k,0) - F(j  ,k,0));
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_ey_f_yee(int j, int k, int l, Array4<Real> const& Ey,
+                         Array4<Real const> const& F, Real dtsdy_c2)
+{
+#if defined WARPX_DIM_3D
+    Ey(j,k,l) += + dtsdy_c2 * (F(j,k+1,l) - F(j,k,l));
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_ez_f_yee(int j, int k, int l, Array4<Real> const& Ez,
+                         Array4<Real const> const& F, Real dtsdz_c2)
+{
+#if defined WARPX_DIM_3D
+    Ez(j,k,l) += + dtsdz_c2 * (F(j,k,l+1) - F(j,k,l));
+#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ)
+    Ez(j,k,0) += + dtsdz_c2 * (F(j,k+1,0) - F(j,k,0));
+#endif
+}
+
+static void warpx_calculate_ckc_coefficients(Real dtsdx, Real dtsdy, Real dtsdz,
+                                             Real &betaxy, Real &betaxz, Real &betayx, Real &betayz, Real &betazx, Real &betazy,
+                                             Real &gammax, Real &gammay, Real &gammaz,
+                                             Real &alphax, Real &alphay, Real &alphaz)
+{
+    // Cole-Karkkainen-Cowan push
+    // computes coefficients according to Cowan - PRST-AB 16, 041303 (2013)
+#if defined WARPX_DIM_3D
+    const Real delta = std::max( { dtsdx,dtsdy,dtsdz } );
+    const Real rx = (dtsdx/delta)*(dtsdx/delta);
+    const Real ry = (dtsdy/delta)*(dtsdy/delta);
+    const Real rz = (dtsdz/delta)*(dtsdz/delta);
+    const Real beta = 0.125*(1. - rx*ry*rz/(ry*rz + rz*rx + rx*ry));
+    betaxy = ry*beta;
+    betaxz = rz*beta;
+    betayx = rx*beta;
+    betayz = rz*beta;
+    betazx = rx*beta;
+    betazy = ry*beta;
+    gammax = ry*rz*(0.0625 - 0.125*ry*rz/(ry*rz + rz*rx + rx*ry));
+    gammay = rx*rz*(0.0625 - 0.125*rx*rz/(ry*rz + rz*rx + rx*ry));
+    gammaz = rx*ry*(0.0625 - 0.125*rx*ry/(ry*rz + rz*rx + rx*ry));
+    alphax = 1. - 2.*betaxy - 2.*betaxz - 4.*gammax;
+    alphay = 1. - 2.*betayx - 2.*betayz - 4.*gammay;
+    alphaz = 1. - 2.*betazx - 2.*betazy - 4.*gammaz;
+
+    betaxy *= dtsdx;
+    betaxz *= dtsdx;
+    betayx *= dtsdy;
+    betayz *= dtsdy;
+    betazx *= dtsdz;
+    betazy *= dtsdz;
+    alphax *= dtsdx;
+    alphay *= dtsdy;
+    alphaz *= dtsdz;
+    gammax *= dtsdx;
+    gammay *= dtsdy;
+    gammaz *= dtsdz;
+#elif defined WARPX_DIM_2D
+    const Real delta = std::max(dtsdx,dtsdz);
+    const Real rx = (dtsdx/delta)*(dtsdx/delta);
+    const Real rz = (dtsdz/delta)*(dtsdz/delta);
+    betaxz = 0.125*rz;
+    betazx = 0.125*rx;
+    alphax = 1. - 2.*betaxz;
+    alphaz = 1. - 2.*betazx;
+
+    betaxz *= dtsdx;
+    betazx *= dtsdz;
+    alphax *= dtsdx;
+    alphaz *= dtsdz;
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_bx_ckc(int j, int k, int l, Array4<Real> const& Bx,
+                       Array4<Real const> const& Ey, Array4<Real const> const& Ez,
+                       Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy,
+                       Real gammax, Real gammay, Real gammaz,
+                       Real alphax, Real alphay, Real alphaz)
+{
+#if defined WARPX_DIM_3D
+    Bx(j,k,l) += - alphay * (Ez(j  ,k+1,l  ) - Ez(j,  k  ,l  ))
+                 - betayx * (Ez(j+1,k+1,l  ) - Ez(j+1,k  ,l  )
+                          +  Ez(j-1,k+1,l  ) - Ez(j-1,k  ,l  ))
+                 - betayz * (Ez(j  ,k+1,l+1) - Ez(j  ,k  ,l+1)
+                          +  Ez(j  ,k+1,l-1) - Ez(j  ,k  ,l-1))
+                 - gammay * (Ez(j+1,k+1,l+1) - Ez(j+1,k  ,l+1)
+                          +  Ez(j-1,k+1,l+1) - Ez(j-1,k  ,l+1)
+                          +  Ez(j+1,k+1,l-1) - Ez(j+1,k  ,l-1)
+                          +  Ez(j-1,k+1,l-1) - Ez(j-1,k  ,l-1))
+                 + alphaz * (Ey(j  ,k  ,l+1) - Ey(j,  k,  l  ))
+                 + betazx * (Ey(j+1,k  ,l+1) - Ey(j+1,k  ,l  )
+                          +  Ey(j-1,k  ,l+1) - Ey(j-1,k  ,l  ))
+                 + betazy * (Ey(j  ,k+1,l+1) - Ey(j  ,k+1,l  )
+                          +  Ey(j  ,k-1,l+1) - Ey(j  ,k-1,l  ))
+                 + gammaz * (Ey(j+1,k+1,l+1) - Ey(j+1,k+1,l  )
+                          +  Ey(j-1,k+1,l+1) - Ey(j-1,k+1,l  )
+                          +  Ey(j+1,k-1,l+1) - Ey(j+1,k-1,l  )
+                          +  Ey(j-1,k-1,l+1) - Ey(j-1,k-1,l  ));
+#elif defined WARPX_DIM_2D
+    Bx(j,k,0) += + alphaz * (Ey(j  ,k+1,0) - Ey(j,  k,0))
+                 + betazx * (Ey(j+1,k+1,0) - Ey(j+1,k,0)
+                          +  Ey(j-1,k+1,0) - Ey(j-1,k,0));
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_by_ckc(int j, int k, int l, Array4<Real> const& By,
+                       Array4<Real const> const& Ex, Array4<Real const> const& Ez,
+                       Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy,
+                       Real gammax, Real gammay, Real gammaz,
+                       Real alphax, Real alphay, Real alphaz)
+{
+#if defined WARPX_DIM_3D
+    By(j,k,l) += + alphax * (Ez(j+1,k  ,l  ) - Ez(j,  k,  l  ))
+                 + betaxy * (Ez(j+1,k+1,l  ) - Ez(j  ,k+1,l  )
+                          +  Ez(j+1,k-1,l  ) - Ez(j  ,k-1,l  ))
+                 + betaxz * (Ez(j+1,k  ,l+1) - Ez(j  ,k  ,l+1)
+                          +  Ez(j+1,k  ,l-1) - Ez(j  ,k  ,l-1))
+                 + gammax * (Ez(j+1,k+1,l+1) - Ez(j  ,k+1,l+1)
+                          +  Ez(j+1,k-1,l+1) - Ez(j  ,k-1,l+1)
+                          +  Ez(j+1,k+1,l-1) - Ez(j  ,k+1,l-1)
+                          +  Ez(j+1,k-1,l-1) - Ez(j  ,k-1,l-1))
+                 - alphaz * (Ex(j  ,k  ,l+1) - Ex(j  ,k  ,l  ))
+                 - betazx * (Ex(j+1,k  ,l+1) - Ex(j+1,k  ,l  )
+                          +  Ex(j-1,k  ,l+1) - Ex(j-1,k  ,l  ))
+                 - betazy * (Ex(j  ,k+1,l+1) - Ex(j  ,k+1,l  )
+                          +  Ex(j  ,k-1,l+1) - Ex(j  ,k-1,l  ))
+                 - gammaz * (Ex(j+1,k+1,l+1) - Ex(j+1,k+1,l  )
+                          +  Ex(j-1,k+1,l+1) - Ex(j-1,k+1,l  )
+                          +  Ex(j+1,k-1,l+1) - Ex(j+1,k-1,l  )
+                          +  Ex(j-1,k-1,l+1) - Ex(j-1,k-1,l  ));
+#elif defined WARPX_DIM_2D
+    By(j,k,0) += + alphax * (Ez(j+1,k  ,0) - Ez(j,k  ,0))
+                 + betaxz * (Ez(j+1,k+1,0) - Ez(j,k+1,0)
+                          +  Ez(j+1,k-1,0) - Ez(j,k-1,0))
+                 - alphaz * (Ex(j  ,k+1,0) - Ex(j,k  ,0))
+                 - betazx * (Ex(j+1,k+1,0) - Ex(j+1,k,0)
+                          +  Ex(j-1,k+1,0) - Ex(j-1,k,0));
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_bz_ckc(int j, int k, int l, Array4<Real> const& Bz,
+                       Array4<Real const> const& Ex, Array4<Real const> const& Ey,
+                       Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy,
+                       Real gammax, Real gammay, Real gammaz,
+                       Real alphax, Real alphay, Real alphaz)
+{
+#if defined WARPX_DIM_3D
+    Bz(j,k,l) += - alphax * (Ey(j+1,k  ,l  ) - Ey(j  ,k  ,l  ))
+                 - betaxy * (Ey(j+1,k+1,l  ) - Ey(j  ,k+1,l  )
+                          +  Ey(j+1,k-1,l  ) - Ey(j  ,k-1,l  ))
+                 - betaxz * (Ey(j+1,k  ,l+1) - Ey(j  ,k  ,l+1)
+                          +  Ey(j+1,k  ,l-1) - Ey(j  ,k  ,l-1))
+                 - gammax * (Ey(j+1,k+1,l+1) - Ey(j  ,k+1,l+1)
+                          +  Ey(j+1,k-1,l+1) - Ey(j  ,k-1,l+1)
+                          +  Ey(j+1,k+1,l-1) - Ey(j  ,k+1,l-1)
+                          +  Ey(j+1,k-1,l-1) - Ey(j  ,k-1,l-1))
+                 + alphay * (Ex(j  ,k+1,l  ) - Ex(j  ,k  ,l  ))
+                 + betayx * (Ex(j+1,k+1,l  ) - Ex(j+1,k  ,l  )
+                          +  Ex(j-1,k+1,l  ) - Ex(j-1,k  ,l  ))
+                 + betayz * (Ex(j  ,k+1,l+1) - Ex(j  ,k  ,l+1)
+                          +  Ex(j  ,k+1,l-1) - Ex(j  ,k  ,l-1))
+                 + gammay * (Ex(j+1,k+1,l+1) - Ex(j+1,k  ,l+1)
+                          +  Ex(j-1,k+1,l+1) - Ex(j-1,k  ,l+1)
+                          +  Ex(j+1,k+1,l-1) - Ex(j+1,k  ,l-1)
+                          +  Ex(j-1,k+1,l-1) - Ex(j-1,k  ,l-1));
+#elif defined WARPX_DIM_2D
+    Bz(j,k,0) += - alphax * (Ey(j+1,k  ,0) - Ey(j,k  ,0))
+                 - betaxz * (Ey(j+1,k+1,0) - Ey(j,k+1,0)
+                          +  Ey(j+1,k-1,0) - Ey(j,k-1,0));
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_ex_f_ckc(int j, int k, int l, Array4<Real> const& Ex,
+                         Array4<Real const> const& F,
+                         Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy,
+                         Real gammax, Real gammay, Real gammaz,
+                         Real alphax, Real alphay, Real alphaz)
+{
+#if defined WARPX_DIM_3D
+    Ex(j,k,l) += + alphax * (F(j+1,k  ,l  ) - F(j,  k,  l  ))
+                 + betaxy * (F(j+1,k+1,l  ) - F(j,k+1,l  )
+                          +  F(j+1,k-1,l  ) - F(j,k-1,l  ))
+                 + betaxz * (F(j+1,k  ,l+1) - F(j,k  ,l+1)
+                          +  F(j+1,k  ,l-1) - F(j,k  ,l-1))
+                 + gammax * (F(j+1,k+1,l+1) - F(j,k+1,l+1)
+                          +  F(j+1,k-1,l+1) - F(j,k-1,l+1)
+                          +  F(j+1,k+1,l-1) - F(j,k+1,l-1)
+                          +  F(j+1,k-1,l-1) - F(j,k-1,l-1));
+#elif defined WARPX_DIM_2D
+    Ex(j,k,0) += + alphax * (F(j+1,k  ,0) - F(j,k  ,0))
+                 + betaxz * (F(j+1,k+1,0) - F(j,k+1,0)
+                          +  F(j+1,k-1,0) - F(j,k-1,0));
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_ey_f_ckc(int j, int k, int l, Array4<Real> const& Ey,
+                         Array4<Real const> const& F,
+                         Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy,
+                         Real gammax, Real gammay, Real gammaz,
+                         Real alphax, Real alphay, Real alphaz)
+{
+#if defined WARPX_DIM_3D
+    Ey(j,k,l) += + alphay * (F(j  ,k+1,l  ) - F(j  ,k,l  ))
+                 + betayx * (F(j+1,k+1,l  ) - F(j+1,k,l  )
+                          +  F(j-1,k+1,l  ) - F(j-1,k,l  ))
+                 + betayz * (F(j  ,k+1,l+1) - F(j  ,k,l+1)
+                          +  F(j  ,k+1,l-1) - F(j  ,k,l-1))
+                 + gammay * (F(j+1,k+1,l+1) - F(j+1,k,l+1)
+                          +  F(j-1,k+1,l+1) - F(j-1,k,l+1)
+                          +  F(j+1,k+1,l-1) - F(j+1,k,l-1)
+                          +  F(j-1,k+1,l-1) - F(j-1,k,l-1));
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void warpx_push_ez_f_ckc(int j, int k, int l, Array4<Real> const& Ez,
+                         Array4<Real const> const& F,
+                         Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy,
+                         Real gammax, Real gammay, Real gammaz,
+                         Real alphax, Real alphay, Real alphaz)
+{
+#if defined WARPX_DIM_3D
+    Ez(j,k,l) += + alphaz * (F(j  ,k  ,l+1) - F(j,  k,  l))
+                 + betazx * (F(j+1,k  ,l+1) - F(j+1,k  ,l)
+                          +  F(j-1,k  ,l+1) - F(j-1,k  ,l))
+                 + betazy * (F(j  ,k+1,l+1) - F(j  ,k+1,l)
+                          +  F(j  ,k-1,l+1) - F(j  ,k-1,l))
+                 + gammaz * (F(j+1,k+1,l+1) - F(j+1,k+1,l)
+                          +  F(j-1,k+1,l+1) - F(j-1,k+1,l)
+                          +  F(j+1,k-1,l+1) - F(j+1,k-1,l)
+                          +  F(j-1,k-1,l+1) - F(j-1,k-1,l));
+#elif defined WARPX_DIM_2D
+    Ez(j,k,0) += + alphaz * (F(j  ,k+1,0) - F(j  ,k,0))
+                 + betazx * (F(j+1,k+1,0) - F(j+1,k,0)
+                          +  F(j-1,k+1,0) - F(j-1,k,0));
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void warpx_computedivb(int i, int j, int k, int dcomp, Array4<Real> const& divB,
                        Array4<Real const> const& Bx, Array4<Real const> const& By, Array4<Real const> const& Bz,
                        Real dxinv, Real dyinv, Real dzinv
-#ifdef WARPX_RZ
+#ifdef WARPX_DIM_RZ
                        ,const Real rmin
 #endif
                        )
 {
-#if (AMREX_SPACEDIM == 3)
+#if defined WARPX_DIM_3D
     divB(i,j,k,dcomp) = (Bx(i+1,j  ,k  ) - Bx(i,j,k))*dxinv
                +        (By(i  ,j+1,k  ) - By(i,j,k))*dyinv
                +        (Bz(i  ,j  ,k+1) - Bz(i,j,k))*dzinv;
-#else
-#ifndef WARPX_RZ
+#elif defined WARPX_DIM_2D
     divB(i,j,0,dcomp) = (Bx(i+1,j  ,0) - Bx(i,j,0))*dxinv
                +        (Bz(i  ,j+1,0) - Bz(i,j,0))*dzinv;
-#else
-    Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5);
-    Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5);
+#elif defined WARPX_DIM_RZ
+    const Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5);
+    const Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5);
     divB(i,j,0,dcomp) = (ru*Bx(i+1,j,0) - rd*Bx(i,j,0))*dxinv
                        + (Bz(i,j+1,0) - Bz(i,j,0))*dzinv;
 #endif
-#endif
 }
 
-AMREX_GPU_HOST_DEVICE AMREX_INLINE
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void warpx_computedive(int i, int j, int k, int dcomp, Array4<Real> const& divE,
                        Array4<Real const> const& Ex, Array4<Real const> const& Ey, Array4<Real const> const& Ez,
                        Real dxinv, Real dyinv, Real dzinv
-#ifdef WARPX_RZ
+#ifdef WARPX_DIM_RZ
                        ,const Real rmin
 #endif
                        )
 {
-#if (AMREX_SPACEDIM == 3)
+#if defined WARPX_DIM_3D
     divE(i,j,k,dcomp) = (Ex(i,j,k) - Ex(i-1,j,k))*dxinv
                +        (Ey(i,j,k) - Ey(i,j-1,k))*dyinv
                +        (Ez(i,j,k) - Ez(i,j,k-1))*dzinv;
-#else
-#ifndef WARPX_RZ
+#elif defined WARPX_DIM_2D
     divE(i,j,0,dcomp) = (Ex(i,j,0) - Ex(i-1,j,0))*dxinv
                  +      (Ez(i,j,0) - Ez(i,j-1,0))*dzinv;
-#else
+#elif defined WARPX_DIM_RZ
     if (i == 0 && rmin == 0.) {
         // the bulk equation diverges on axis
         // (due to the 1/r terms). The following expressions regularize
@@ -56,13 +417,12 @@ void warpx_computedive(int i, int j, int k, int dcomp, Array4<Real> const& divE,
         divE(i,j,0,dcomp) = 4.*Ex(i,j,0)*dxinv
                             + (Ez(i,j,0) - Ez(i,j-1,0))*dzinv;
     } else {
-        Real ru = 1. + 0.5/(rmin*dxinv + i);
-        Real rd = 1. - 0.5/(rmin*dxinv + i);
+        const Real ru = 1. + 0.5/(rmin*dxinv + i);
+        const Real rd = 1. - 0.5/(rmin*dxinv + i);
         divE(i,j,0,dcomp) = (ru*Ex(i,j,0) - rd*Ex(i-1,j,0))*dxinv
                            + (Ez(i,j,0) - Ez(i,j-1,0))*dzinv;
     }
 #endif
-#endif
 }
 
 #endif
diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90
index 33f85c633..dc47245dd 100644
--- a/Source/FortranInterface/WarpX_picsar.F90
+++ b/Source/FortranInterface/WarpX_picsar.F90
@@ -2,33 +2,20 @@
 
 #define WRPX_PXR_GETEB_ENERGY_CONSERVING  geteb3d_energy_conserving_generic
 #define WRPX_PXR_CURRENT_DEPOSITION       depose_jxjyjz_generic
-#define WRPX_PXR_PUSH_BVEC                pxrpush_em3d_bvec
-#define WRPX_PXR_PUSH_BVEC_CKC            pxrpush_em3d_bvec_ckc
-#define WRPX_PXR_PUSH_EVEC_F              pxrpush_em3d_evec_f
-#define WRPX_PXR_PUSH_EVEC_F_CKC          pxrpush_em3d_evec_f_ckc
-#define WRPX_PXR_PUSH_EVEC                pxrpush_em3d_evec
 
 #elif (AMREX_SPACEDIM == 2)
 
-#define WRPX_PXR_PUSH_BVEC_CKC            pxrpush_em2d_bvec_ckc
-#define WRPX_PXR_PUSH_EVEC_F              pxrpush_em2d_evec_f
-#define WRPX_PXR_PUSH_EVEC_F_CKC          pxrpush_em2d_evec_f_ckc
-
 #ifdef WARPX_RZ
 
 #define WRPX_PXR_GETEB_ENERGY_CONSERVING  geteb2drz_energy_conserving_generic
 #define WRPX_PXR_CURRENT_DEPOSITION       depose_jrjtjz_generic_rz
 #define WRPX_PXR_RZ_VOLUME_SCALING_RHO    apply_rz_volume_scaling_rho
 #define WRPX_PXR_RZ_VOLUME_SCALING_J      apply_rz_volume_scaling_j
-#define WRPX_PXR_PUSH_BVEC                pxrpush_emrz_bvec
-#define WRPX_PXR_PUSH_EVEC                pxrpush_emrz_evec
 
 #else
 
 #define WRPX_PXR_GETEB_ENERGY_CONSERVING  geteb2dxz_energy_conserving_generic
 #define WRPX_PXR_CURRENT_DEPOSITION       depose_jxjyjz_generic_2d
-#define WRPX_PXR_PUSH_BVEC                pxrpush_em2d_bvec
-#define WRPX_PXR_PUSH_EVEC                pxrpush_em2d_evec
 
 #endif
 
@@ -529,193 +516,4 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n
 
   end subroutine warpx_particle_pusher_momenta
 
-  ! _________________________________________________________________
-  !>
-  !> @brief
-  !> Main subroutine for the Maxwell solver (E field)
-  !>
-  !> @param[in] xlo, xhi, ylo, yhi, zlo, zhi lower and higher bounds
-  !>            where to apply the solver (PICSAR's valid cells)
-  !> @param[inout] ex, ey, ez arrays for electric field to update
-  !> @param[in] bx, by, bz, jx, jy, jz arrays for magnetic field and current
-  !> @param[in] exlo, exhi, eylo, eyhi, ezlo, ezhi lower and higher bound for
-  !>            the electric field arrays
-  !> @param[in] bxlo, bxhi, bylo, byhi, bzlo, bzhi lower and higher bound for
-  !>            the magnetic field arrays
-  !> @param[in] jxlo, jxhi, jylo, jyhi, jzlo, jzhi lower and higher bound for
-  !>            the current arrays
-  !> @param[in] mudt normalized time step mu_0 * c**2 * dt
-  !> @param[in] dtsdx, dtsdy, dtsdz factors c**2 * dt/(dx, dy, dz)
-  subroutine warpx_push_evec( &
-    xlo, xhi, ylo, yhi, zlo, zhi, &
-    ex, exlo, exhi, &
-    ey, eylo, eyhi, &
-    ez, ezlo, ezhi, &
-    bx, bxlo, bxhi, &
-    by, bylo, byhi, &
-    bz, bzlo, bzhi, &
-    jx, jxlo, jxhi, &
-    jy, jylo, jyhi, &
-    jz, jzlo, jzhi, &
-    mudt, dtsdx, dtsdy, dtsdz, xmin, dx) bind(C, name="warpx_push_evec")
-
-    integer(c_int), intent(in) :: xlo(BL_SPACEDIM), xhi(BL_SPACEDIM), &
-      ylo(BL_SPACEDIM), yhi(BL_SPACEDIM), zlo(BL_SPACEDIM), zhi(BL_SPACEDIM), &
-      exlo(BL_SPACEDIM), exhi(BL_SPACEDIM), eylo(BL_SPACEDIM), eyhi(BL_SPACEDIM), &
-      ezlo(BL_SPACEDIM), ezhi(BL_SPACEDIM), bxlo(BL_SPACEDIM), bxhi(BL_SPACEDIM), &
-      bylo(BL_SPACEDIM), byhi(BL_SPACEDIM), bzlo(BL_SPACEDIM), bzhi(BL_SPACEDIM), &
-      jxlo(BL_SPACEDIM), jxhi(BL_SPACEDIM), jylo(BL_SPACEDIM), jyhi(BL_SPACEDIM), &
-      jzlo(BL_SPACEDIM), jzhi(BL_SPACEDIM)
-
-    real(amrex_real), intent(IN OUT):: ex(*), ey(*), ez(*)
-
-    real(amrex_real), intent(IN):: bx(*), by(*), bz(*), jx(*), jy(*), jz(*)
-
-    real(amrex_real), intent(IN) :: mudt, dtsdx, dtsdy, dtsdz
-
-    real(amrex_real), intent(IN) :: xmin, dx
-
-    CALL WRPX_PXR_PUSH_EVEC(&
-        xlo, xhi, ylo, yhi, zlo, zhi, &
-        ex, exlo, exhi,&
-        ey, eylo, eyhi,&
-        ez, ezlo, ezhi,&
-        bx, bxlo, bxhi,&
-        by, bylo, byhi,&
-        bz, bzlo, bzhi,&
-        jx, jxlo, jxhi,&
-        jy, jylo, jyhi,&
-        jz, jzlo, jzhi,&
-        mudt, dtsdx, dtsdy, dtsdz &
-#ifdef WARPX_RZ
-        ,xmin,dx &
-#endif
-        )
-
-  end subroutine warpx_push_evec
-
-  ! _________________________________________________________________
-  !>
-  !> @brief
-  !> Main subroutine for the Maxwell solver (B field)
-  !>
-  !> @param[in] xlo, xhi, ylo, yhi, zlo, zhi lower and higher bounds
-  !>            where to apply the solver (PICSAR's valid cells)
-  !> @param[in] ex, ey, ez arrays for electric field
-  !> @param[inout] bx, by, bz arrays for magnetic field to update
-  !> @param[in] exlo, exhi, eylo, eyhi, ezlo, ezhi lower and higher bound for
-  !>            the electric field arrays
-  !> @param[in] bxlo, bxhi, bylo, byhi, bzlo, bzhi lower and higher bound for
-  !>            the magnetic field arrays
-  !> @param[in] dtsdx, dtsdy, dtsdz factors 0.5 * dt/(dx, dy, dz)
-  subroutine warpx_push_bvec( &
-    xlo, xhi, ylo, yhi, zlo, zhi, &
-    ex, exlo, exhi, &
-    ey, eylo, eyhi, &
-    ez, ezlo, ezhi, &
-    bx, bxlo, bxhi, &
-    by, bylo, byhi, &
-    bz, bzlo, bzhi, &
-    dtsdx, dtsdy, dtsdz, xmin, dx, &
-    maxwell_fdtd_solver_id) bind(C, name="warpx_push_bvec")
-
-    integer(c_int), intent(in) :: xlo(BL_SPACEDIM), xhi(BL_SPACEDIM), &
-      ylo(BL_SPACEDIM), yhi(BL_SPACEDIM), zlo(BL_SPACEDIM), zhi(BL_SPACEDIM), &
-      exlo(BL_SPACEDIM), exhi(BL_SPACEDIM), eylo(BL_SPACEDIM), eyhi(BL_SPACEDIM), &
-      ezlo(BL_SPACEDIM), ezhi(BL_SPACEDIM), bxlo(BL_SPACEDIM), bxhi(BL_SPACEDIM), &
-      bylo(BL_SPACEDIM), byhi(BL_SPACEDIM), bzlo(BL_SPACEDIM), bzhi(BL_SPACEDIM), &
-      maxwell_fdtd_solver_id
-
-    real(amrex_real), intent(IN OUT):: ex(*), ey(*), ez(*)
-
-    real(amrex_real), intent(IN):: bx(*), by(*), bz(*)
-
-    real(amrex_real), intent(IN) :: dtsdx, dtsdy, dtsdz
-
-    real(amrex_real), intent(IN) :: xmin, dx
-
-    IF (maxwell_fdtd_solver_id .eq. 0) THEN
-      ! Yee FDTD solver
-      CALL WRPX_PXR_PUSH_BVEC( &
-        xlo, xhi, ylo, yhi, zlo, zhi, &
-      	ex, exlo, exhi, &
-      	ey, eylo, eyhi, &
-      	ez, ezlo, ezhi, &
-      	bx, bxlo, bxhi, &
-      	by, bylo, byhi, &
-      	bz, bzlo, bzhi, &
-      	dtsdx,dtsdy,dtsdz &
-#ifdef WARPX_RZ
-        ,xmin,dx &
-#endif
-        )
-    ELSE IF (maxwell_fdtd_solver_id .eq. 1) THEN
-      ! Cole-Karkkainen FDTD solver
-      CALL WRPX_PXR_PUSH_BVEC_CKC( &
-        xlo, xhi, ylo, yhi, zlo, zhi, &
-      	ex, exlo, exhi, &
-      	ey, eylo, eyhi, &
-      	ez, ezlo, ezhi, &
-      	bx, bxlo, bxhi, &
-      	by, bylo, byhi, &
-      	bz, bzlo, bzhi, &
-      	dtsdx,dtsdy,dtsdz)
-    ENDIF
-  end subroutine warpx_push_bvec
-
-  ! _________________________________________________________________
-  !>
-  !> @brief
-  !> Subroutine for pushing E with the grad F terms
-  !>
-  !> @param[in] xlo, xhi, ylo, yhi, zlo, zhi lower and higher bounds
-  !>            where to apply the solver (PICSAR's valid cells)
-  !> @param[inout] ex, ey, ez arrays for electric field to update
-  !> @param[in] f array for the charge-correction term
-  !> @param[in] exlo, exhi, eylo, eyhi, ezlo, ezhi lower and higher bound for
-  !>            the electric field arrays
-  !> @param[in] flo, fhi, lower and higher bound for the charge-correction term
-  !> @param[in] dtsdx_c2, dtsdy_c2, dtsdz_c2 factors dt/(dx, dy, dz)*c2
-  subroutine warpx_push_evec_f( &
-    xlo, xhi, ylo, yhi, zlo, zhi, &
-    ex, exlo, exhi, &
-    ey, eylo, eyhi, &
-    ez, ezlo, ezhi, &
-    f, flo, fhi, &
-    dtsdx_c2, dtsdy_c2, dtsdz_c2, &
-    maxwell_fdtd_solver_id) bind(C, name="warpx_push_evec_f")
-
-    integer(c_int), intent(in) :: xlo(BL_SPACEDIM), xhi(BL_SPACEDIM), &
-      ylo(BL_SPACEDIM), yhi(BL_SPACEDIM), zlo(BL_SPACEDIM), zhi(BL_SPACEDIM), &
-      exlo(BL_SPACEDIM), exhi(BL_SPACEDIM), eylo(BL_SPACEDIM), eyhi(BL_SPACEDIM), &
-      ezlo(BL_SPACEDIM), ezhi(BL_SPACEDIM), flo(BL_SPACEDIM), fhi(BL_SPACEDIM), &
-      maxwell_fdtd_solver_id
-
-    real(amrex_real), intent(IN OUT):: ex(*), ey(*), ez(*)
-
-    real(amrex_real), intent(IN):: f
-
-    real(amrex_real), intent(IN) :: dtsdx_c2, dtsdy_c2, dtsdz_c2
-
-    IF (maxwell_fdtd_solver_id .eq. 0) THEN
-      ! Yee FDTD solver
-      CALL WRPX_PXR_PUSH_EVEC_F( &
-        xlo, xhi, ylo, yhi, zlo, zhi, &
-      	ex, exlo, exhi, &
-      	ey, eylo, eyhi, &
-      	ez, ezlo, ezhi, &
-      	f, flo, fhi, &
-      	dtsdx_c2, dtsdy_c2, dtsdz_c2)
-    ELSE IF (maxwell_fdtd_solver_id .eq. 1) THEN
-      ! Cole-Karkkainen FDTD solver
-      CALL WRPX_PXR_PUSH_EVEC_F_CKC( &
-        xlo, xhi, ylo, yhi, zlo, zhi, &
-      	ex, exlo, exhi, &
-      	ey, eylo, eyhi, &
-      	ez, ezlo, ezhi, &
-      	f, flo, fhi, &
-      	dtsdx_c2, dtsdy_c2, dtsdz_c2)
-    ENDIF
-  end subroutine warpx_push_evec_f
-
 end module warpx_to_pxr_module