Merge branch 'dev' into picsar_cleanup

1 files changed, 114 insertions, 1 deletions

diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp
index b0688b2c1..f5c231ddf 100644
--- a/Source/BoundaryConditions/WarpXEvolvePML.cpp
+++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp
@@ -12,6 +12,8 @@
 #include <AMReX_AmrMeshInSituBridge.H>
 #endif
 
+#include <PML_current.H>
+
 using namespace amrex;
 
 void
@@ -55,7 +57,6 @@ WarpX::DampPML (int lev, PatchType patch_type)
             const Box& tbx  = mfi.tilebox(Bx_nodal_flag);
             const Box& tby  = mfi.tilebox(By_nodal_flag);
             const Box& tbz  = mfi.tilebox(Bz_nodal_flag);
-
             WRPX_DAMP_PML(tex.loVect(), tex.hiVect(),
 			    tey.loVect(), tey.hiVect(),
 			    tez.loVect(), tez.hiVect(),
@@ -79,3 +80,115 @@ WarpX::DampPML (int lev, PatchType patch_type)
         }
     }
 }
+
+void
+WarpX::DampJPML ()
+{
+    for (int lev = 0; lev <= finest_level; ++lev) {
+        DampJPML(lev);
+    }
+}
+
+void
+WarpX::DampJPML (int lev)
+{
+    DampJPML(lev, PatchType::fine);
+    if (lev > 0) DampJPML(lev, PatchType::coarse);
+}
+
+void
+WarpX::DampJPML (int lev, PatchType patch_type)
+{
+    if (!do_pml) return;
+    if (!do_pml_j_damping) return;
+
+    BL_PROFILE("WarpX::DampJPML()");
+
+    if (pml[lev]->ok())
+    {
+
+        const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp();
+        const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp()
+                                                            : pml[lev]->GetMultiSigmaBox_cp();
+
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+        for ( MFIter mfi(*pml_j[0], TilingIfNotGPU()); mfi.isValid(); ++mfi )
+        {
+            auto const& pml_jxfab = pml_j[0]->array(mfi);
+            auto const& pml_jyfab = pml_j[1]->array(mfi);
+            auto const& pml_jzfab = pml_j[2]->array(mfi);
+            const Real* sigma_cumsum_fac_j_x = sigba[mfi].sigma_cumsum_fac[0].data();
+            const Real* sigma_star_cumsum_fac_j_x = sigba[mfi].sigma_star_cumsum_fac[0].data();
+#if (AMREX_SPACEDIM == 3)
+            const Real* sigma_cumsum_fac_j_y = sigba[mfi].sigma_cumsum_fac[1].data();
+            const Real* sigma_star_cumsum_fac_j_y = sigba[mfi].sigma_star_cumsum_fac[1].data();
+            const Real* sigma_cumsum_fac_j_z = sigba[mfi].sigma_cumsum_fac[2].data();
+            const Real* sigma_star_cumsum_fac_j_z = sigba[mfi].sigma_star_cumsum_fac[2].data();
+#else
+            const Real* sigma_cumsum_fac_j_y = nullptr;
+            const Real* sigma_star_cumsum_fac_j_y = nullptr;
+            const Real* sigma_cumsum_fac_j_z = sigba[mfi].sigma_cumsum_fac[1].data();
+            const Real* sigma_star_cumsum_fac_j_z = sigba[mfi].sigma_star_cumsum_fac[1].data();
+#endif
+            const Box& tjx  = mfi.tilebox(jx_nodal_flag);
+            const Box& tjy  = mfi.tilebox(jy_nodal_flag);
+            const Box& tjz  = mfi.tilebox(jz_nodal_flag);
+
+            auto const& AMREX_RESTRICT x_lo = sigba[mfi].sigma_cumsum_fac[0].lo();
+#if (AMREX_SPACEDIM == 3)
+            auto const& AMREX_RESTRICT y_lo = sigba[mfi].sigma_cumsum_fac[1].lo();
+            auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma_cumsum_fac[2].lo();
+#else
+            int y_lo = 0;
+            auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma_cumsum_fac[1].lo();
+#endif
+
+            auto const& AMREX_RESTRICT xs_lo = sigba[mfi].sigma_star_cumsum_fac[0].lo();
+#if (AMREX_SPACEDIM == 3)
+            auto const& AMREX_RESTRICT ys_lo = sigba[mfi].sigma_star_cumsum_fac[1].lo();
+            auto const& AMREX_RESTRICT zs_lo = sigba[mfi].sigma_star_cumsum_fac[2].lo();
+#else
+            int ys_lo = 0;
+            auto const& AMREX_RESTRICT zs_lo = sigba[mfi].sigma_star_cumsum_fac[1].lo();
+#endif
+
+            amrex::ParallelFor( tjx, tjy, tjz,
+                [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+                    damp_jx_pml(i, j, k, pml_jxfab, sigma_star_cumsum_fac_j_x,
+                                sigma_cumsum_fac_j_y, sigma_cumsum_fac_j_z,
+                                xs_lo,y_lo, z_lo);
+                },
+                [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+                    damp_jy_pml(i, j, k, pml_jyfab, sigma_cumsum_fac_j_x,
+                                sigma_star_cumsum_fac_j_y, sigma_cumsum_fac_j_z,
+                                x_lo,ys_lo, z_lo);
+                },
+                [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+                    damp_jz_pml(i, j, k, pml_jzfab, sigma_cumsum_fac_j_x,
+                                sigma_cumsum_fac_j_y, sigma_star_cumsum_fac_j_z,
+                                x_lo,y_lo, zs_lo);
+                }
+            );
+        }
+
+    }
+}
+
+/* \brief Copy the current J from the regular grid to the PML */
+void
+WarpX::CopyJPML ()
+{
+    for (int lev = 0; lev <= finest_level; ++lev)
+    {
+        if (pml[lev]->ok()){
+            pml[lev]->CopyJtoPMLs({ current_fp[lev][0].get(),
+                                  current_fp[lev][1].get(),
+                                  current_fp[lev][2].get() },
+                                { current_cp[lev][0].get(),
+                                  current_cp[lev][1].get(),
+                                  current_cp[lev][2].get() });
+        }
+    }
+}