diff options
Diffstat (limited to 'Source/BoundaryConditions/WarpXEvolvePML.cpp')
| -rw-r--r-- | Source/BoundaryConditions/WarpXEvolvePML.cpp | 115 |
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() }); + } + } +} |