Merge branch 'dev' into push_momentum

1 files changed, 137 insertions, 61 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]};