diff options
Diffstat (limited to 'Source/BoundaryConditions/PML.cpp')
| -rw-r--r-- | Source/BoundaryConditions/PML.cpp | 346 |
1 files changed, 264 insertions, 82 deletions
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 21d348482..8f8a2608e 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -1,4 +1,3 @@ - #include <PML.H> #include <WarpX.H> #include <WarpXConst.H> @@ -16,54 +15,70 @@ using namespace amrex; namespace { - static void FillLo (int idim, Sigma& sigma, Sigma& sigma_star, + static void FillLo (int idim, Sigma& sigma, Sigma& sigma_cumsum, + Sigma& sigma_star, Sigma& sigma_star_cumsum, const Box& overlap, const Box& grid, Real fac) { int glo = grid.smallEnd(idim); int olo = overlap.smallEnd(idim); int ohi = overlap.bigEnd(idim); int slo = sigma.m_lo; + int shi = sigma.m_hi; int sslo = sigma_star.m_lo; + for (int i = olo; i <= ohi+1; ++i) { Real offset = static_cast<Real>(glo-i); sigma[i-slo] = fac*(offset*offset); + // sigma_cumsum is the analytical integral of sigma function at same points than sigma + sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3.)/PhysConst::c; } + for (int i = olo; i <= ohi; ++i) { Real offset = static_cast<Real>(glo-i) - 0.5; sigma_star[i-sslo] = fac*(offset*offset); + // sigma_star_cumsum is the analytical integral of sigma function at same points than sigma_star + sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3.)/PhysConst::c; } } - static void FillHi (int idim, Sigma& sigma, Sigma& sigma_star, + static void FillHi (int idim, Sigma& sigma, Sigma& sigma_cumsum, + Sigma& sigma_star, Sigma& sigma_star_cumsum, const Box& overlap, const Box& grid, Real fac) { int ghi = grid.bigEnd(idim); int olo = overlap.smallEnd(idim); int ohi = overlap.bigEnd(idim); int slo = sigma.m_lo; + int shi = sigma.m_hi; int sslo = sigma_star.m_lo; for (int i = olo; i <= ohi+1; ++i) { Real offset = static_cast<Real>(i-ghi-1); sigma[i-slo] = fac*(offset*offset); + sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3.)/PhysConst::c; } for (int i = olo; i <= ohi; ++i) { Real offset = static_cast<Real>(i-ghi) - 0.5; sigma_star[i-sslo] = fac*(offset*offset); + sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3.)/PhysConst::c; } } - static void FillZero (int idim, Sigma& sigma, Sigma& sigma_star, const Box& overlap) + static void FillZero (int idim, Sigma& sigma, Sigma& sigma_cumsum, + Sigma& sigma_star, Sigma& sigma_star_cumsum, + const Box& overlap) { int olo = overlap.smallEnd(idim); int ohi = overlap.bigEnd(idim); int slo = sigma.m_lo; int sslo = sigma_star.m_lo; std::fill(sigma.begin()+(olo-slo), sigma.begin()+(ohi+2-slo), 0.0); + std::fill(sigma_cumsum.begin()+(olo-slo), sigma_cumsum.begin()+(ohi+2-slo), 0.0); std::fill(sigma_star.begin()+(olo-sslo), sigma_star.begin()+(ohi+1-sslo), 0.0); + std::fill(sigma_star_cumsum.begin()+(olo-sslo), sigma_star_cumsum.begin()+(ohi+1-sslo), 0.0); } } @@ -77,19 +92,31 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - sigma [idim].resize(sz[idim]+1); - sigma_star [idim].resize(sz[idim] ); - sigma_fac [idim].resize(sz[idim]+1); - sigma_star_fac[idim].resize(sz[idim] ); - - sigma [idim].m_lo = lo[idim]; - sigma [idim].m_hi = hi[idim]+1; - sigma_star [idim].m_lo = lo[idim]; - sigma_star [idim].m_hi = hi[idim]; - sigma_fac [idim].m_lo = lo[idim]; - sigma_fac [idim].m_hi = hi[idim]+1; - sigma_star_fac[idim].m_lo = lo[idim]; - sigma_star_fac[idim].m_hi = hi[idim]; + sigma [idim].resize(sz[idim]+1); + sigma_cumsum [idim].resize(sz[idim]+1); + sigma_star [idim].resize(sz[idim]); + sigma_star_cumsum [idim].resize(sz[idim]); + sigma_fac [idim].resize(sz[idim]+1); + sigma_cumsum_fac [idim].resize(sz[idim]+1); + sigma_star_fac [idim].resize(sz[idim]); + sigma_star_cumsum_fac[idim].resize(sz[idim]); + + sigma [idim].m_lo = lo[idim]; + sigma [idim].m_hi = hi[idim]+1; + sigma_cumsum [idim].m_lo = lo[idim]; + sigma_cumsum [idim].m_hi = hi[idim]+1; + sigma_star [idim].m_lo = lo[idim]; + sigma_star [idim].m_hi = hi[idim]; + sigma_star_cumsum [idim].m_lo = lo[idim]; + sigma_star_cumsum [idim].m_hi = hi[idim]; + sigma_fac [idim].m_lo = lo[idim]; + sigma_fac [idim].m_hi = hi[idim]+1; + sigma_cumsum_fac [idim].m_lo = lo[idim]; + sigma_cumsum_fac [idim].m_hi = hi[idim]+1; + sigma_star_fac [idim].m_lo = lo[idim]; + sigma_star_fac [idim].m_hi = hi[idim]; + sigma_star_cumsum_fac[idim].m_lo = lo[idim]; + sigma_star_cumsum_fac[idim].m_hi = hi[idim]; } Array<Real,AMREX_SPACEDIM> fac; @@ -157,7 +184,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n #endif Box looverlap = lobox & box; if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); + FillLo(idim, sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + looverlap, grid_box, fac[idim]); } Box hibox = amrex::adjCellHi(grid_box, idim, ncell); @@ -167,7 +196,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n #endif Box hioverlap = hibox & box; if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); + FillHi(idim, sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + hioverlap, grid_box, fac[idim]); } if (!looverlap.ok() && !hioverlap.ok()) { @@ -181,8 +212,10 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n const Box& grid_box = grids[gid]; const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box; if (overlap.ok()) { - FillZero(idim, sigma[idim], sigma_star[idim], overlap); - } else { + FillZero(idim, sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], overlap); + } + else { amrex::Abort("SigmaBox::SigmaBox(): side_side_edges, how did this happen?\n"); } } @@ -194,13 +227,17 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n Box lobox = amrex::adjCellLo(grid_box, idim, ncell); Box looverlap = lobox.grow(jdim,ncell).grow(kdim,ncell) & box; if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); + FillLo(idim, sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + looverlap, grid_box, fac[idim]); } Box hibox = amrex::adjCellHi(grid_box, idim, ncell); Box hioverlap = hibox.grow(jdim,ncell).grow(kdim,ncell) & box; if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); + FillHi(idim, sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + hioverlap, grid_box, fac[idim]); } if (!looverlap.ok() && !hioverlap.ok()) { @@ -218,7 +255,8 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box; #endif if (overlap.ok()) { - FillZero(idim, sigma[idim], sigma_star[idim], overlap); + FillZero(idim, sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], overlap); } else { amrex::Abort("SigmaBox::SigmaBox(): side_faces, how did this happen?\n"); } @@ -231,13 +269,17 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n const Box& lobox = amrex::adjCellLo(grid_box, idim, ncell); Box looverlap = lobox & box; if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); + FillLo(idim, sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + looverlap, grid_box, fac[idim]); } const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell); Box hioverlap = hibox & box; if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); + FillHi(idim, sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + hioverlap, grid_box, fac[idim]); } if (!looverlap.ok() && !hioverlap.ok()) { @@ -251,6 +293,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n } } + void SigmaBox::ComputePMLFactorsB (const Real* dx, Real dt) { @@ -259,6 +302,7 @@ SigmaBox::ComputePMLFactorsB (const Real* dx, Real dt) for (int i = 0, N = sigma_star[idim].size(); i < N; ++i) { sigma_star_fac[idim][i] = std::exp(-sigma_star[idim][i]*dt); + sigma_star_cumsum_fac[idim][i] = std::exp(-sigma_star_cumsum[idim][i]*dx[idim]); } } } @@ -271,6 +315,7 @@ SigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) for (int i = 0, N = sigma[idim].size(); i < N; ++i) { sigma_fac[idim][i] = std::exp(-sigma[idim][i]*dt); + sigma_cumsum_fac[idim][i] = std::exp(-sigma_cumsum[idim][i]*dx[idim]); } } } @@ -320,11 +365,35 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, #endif int do_dive_cleaning, int do_moving_window, + int pml_has_particles, int do_pml_in_domain, const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) : m_geom(geom), m_cgeom(cgeom) { - const BoxArray& ba = MakeBoxArray(*geom, grid_ba, ncell, do_pml_Lo, do_pml_Hi); + + // When `do_pml_in_domain` is true, the PML overlap with the last `ncell` of the physical domain + // (instead of extending `ncell` outside of the physical domain) + // In order to implement this, a reduced domain is created here (decreased by ncells in all direction) + // and passed to `MakeBoxArray`, which surrounds it by PML boxes + // (thus creating the PML boxes at the right position, where they overlap with the original domain) + Box domain0 = geom->Domain(); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if ( ! geom->isPeriodic(idim)) { + if (do_pml_Lo[idim]){ + domain0.growLo(idim, -ncell); + } + if (do_pml_Hi[idim]){ + domain0.growHi(idim, -ncell); + } + + } + } + const BoxArray grid_ba_reduced = BoxArray(grid_ba.boxList().intersect(domain0)); + + const BoxArray& ba = (do_pml_in_domain)? + MakeBoxArray(*geom, grid_ba_reduced, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi) : + MakeBoxArray(*geom, grid_ba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi); + if (ba.size() == 0) { m_ok = false; return; @@ -366,6 +435,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, pml_B_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::By_nodal_flag), dm, 2, ngb)); pml_B_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::Bz_nodal_flag), dm, 2, ngb)); + pml_E_fp[0]->setVal(0.0); pml_E_fp[1]->setVal(0.0); pml_E_fp[2]->setVal(0.0); @@ -373,15 +443,30 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, pml_B_fp[1]->setVal(0.0); pml_B_fp[2]->setVal(0.0); + pml_j_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::jx_nodal_flag), dm, 1, ngb)); + pml_j_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::jy_nodal_flag), dm, 1, ngb)); + pml_j_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::jz_nodal_flag), dm, 1, ngb)); + pml_j_fp[0]->setVal(0.0); + pml_j_fp[1]->setVal(0.0); + pml_j_fp[2]->setVal(0.0); + if (do_dive_cleaning) { pml_F_fp.reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()), dm, 3, ngf)); pml_F_fp->setVal(0.0); } - sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta)); + if (do_pml_in_domain){ + sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba_reduced, geom->CellSize(), ncell, delta)); + } + else { + sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta)); + } + #ifdef WARPX_USE_PSATD + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( do_pml_in_domain==false, + "PSATD solver cannot be used with `do_pml_in_domain`."); const bool in_pml = true; // Tells spectral solver to use split-PML equations const RealVect dx{AMREX_D_DECL(geom->CellSize(0), geom->CellSize(1), geom->CellSize(2))}; // Get the cell-centered box, with guard cells @@ -400,7 +485,11 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, BoxArray grid_cba = grid_ba; grid_cba.coarsen(ref_ratio); - const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba, ncell, do_pml_Lo, do_pml_Hi); + const BoxArray grid_cba_reduced = BoxArray(grid_cba.boxList().intersect(domain0)); + + const BoxArray& cba = (do_pml_in_domain) ? + MakeBoxArray(*cgeom, grid_cba_reduced, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi) : + MakeBoxArray(*cgeom, grid_cba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi); DistributionMapping cdm{cba}; @@ -422,9 +511,20 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, { pml_F_cp.reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()), cdm, 3, ngf)); pml_F_cp->setVal(0.0); - } - sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); + } + pml_j_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::jx_nodal_flag), cdm, 1, ngb)); + pml_j_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::jy_nodal_flag), cdm, 1, ngb)); + pml_j_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::jz_nodal_flag), cdm, 1, ngb)); + pml_j_cp[0]->setVal(0.0); + pml_j_cp[1]->setVal(0.0); + pml_j_cp[2]->setVal(0.0); + + if (do_pml_in_domain){ + sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba_reduced, cgeom->CellSize(), ncell, delta)); + } else { + sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); + } #ifdef WARPX_USE_PSATD const bool in_pml = true; // Tells spectral solver to use split-PML equations @@ -439,7 +539,8 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, } BoxArray -PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell, +PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, + int ncell, int do_pml_in_domain, const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) { Box domain = geom.Domain(); @@ -453,14 +554,18 @@ PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, } } } - BoxList bl; for (int i = 0, N = grid_ba.size(); i < N; ++i) { const Box& grid_bx = grid_ba[i]; const IntVect& grid_bx_sz = grid_bx.size(); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell, - "Consider using larger amr.blocking_factor"); + + if (do_pml_in_domain == 0) { + // Make sure that, in the case of several distinct refinement patches, + // the PML cells surrounding these patches cannot overlap + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell, + "Consider using larger amr.blocking_factor"); + } Box bx = grid_bx; bx.grow(ncell); @@ -530,6 +635,12 @@ PML::GetB_fp () } std::array<MultiFab*,3> +PML::Getj_fp () +{ + return {pml_j_fp[0].get(), pml_j_fp[1].get(), pml_j_fp[2].get()}; +} + +std::array<MultiFab*,3> PML::GetE_cp () { return {pml_E_cp[0].get(), pml_E_cp[1].get(), pml_E_cp[2].get()}; @@ -541,6 +652,12 @@ PML::GetB_cp () return {pml_B_cp[0].get(), pml_B_cp[1].get(), pml_B_cp[2].get()}; } +std::array<MultiFab*,3> +PML::Getj_cp () +{ + return {pml_j_cp[0].get(), pml_j_cp[1].get(), pml_j_cp[2].get()}; +} + MultiFab* PML::GetF_fp () { @@ -555,116 +672,181 @@ PML::GetF_cp () void PML::ExchangeB (const std::array<amrex::MultiFab*,3>& B_fp, - const std::array<amrex::MultiFab*,3>& B_cp) + const std::array<amrex::MultiFab*,3>& B_cp, + int do_pml_in_domain) { - ExchangeB(PatchType::fine, B_fp); - ExchangeB(PatchType::coarse, B_cp); + ExchangeB(PatchType::fine, B_fp, do_pml_in_domain); + ExchangeB(PatchType::coarse, B_cp, do_pml_in_domain); } void PML::ExchangeB (PatchType patch_type, - const std::array<amrex::MultiFab*,3>& Bp) + const std::array<amrex::MultiFab*,3>& Bp, + int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0]) { - Exchange(*pml_B_fp[0], *Bp[0], *m_geom); - Exchange(*pml_B_fp[1], *Bp[1], *m_geom); - Exchange(*pml_B_fp[2], *Bp[2], *m_geom); + Exchange(*pml_B_fp[0], *Bp[0], *m_geom, do_pml_in_domain); + Exchange(*pml_B_fp[1], *Bp[1], *m_geom, do_pml_in_domain); + Exchange(*pml_B_fp[2], *Bp[2], *m_geom, do_pml_in_domain); } else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0]) { - Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom); - Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom); - Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom); + Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom, do_pml_in_domain); + Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom, do_pml_in_domain); + Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom, do_pml_in_domain); } } void PML::ExchangeE (const std::array<amrex::MultiFab*,3>& E_fp, - const std::array<amrex::MultiFab*,3>& E_cp) + const std::array<amrex::MultiFab*,3>& E_cp, + int do_pml_in_domain) { - ExchangeE(PatchType::fine, E_fp); - ExchangeE(PatchType::coarse, E_cp); + ExchangeE(PatchType::fine, E_fp, do_pml_in_domain); + ExchangeE(PatchType::coarse, E_cp, do_pml_in_domain); } void PML::ExchangeE (PatchType patch_type, - const std::array<amrex::MultiFab*,3>& Ep) + const std::array<amrex::MultiFab*,3>& Ep, + int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0]) { - Exchange(*pml_E_fp[0], *Ep[0], *m_geom); - Exchange(*pml_E_fp[1], *Ep[1], *m_geom); - Exchange(*pml_E_fp[2], *Ep[2], *m_geom); + Exchange(*pml_E_fp[0], *Ep[0], *m_geom, do_pml_in_domain); + Exchange(*pml_E_fp[1], *Ep[1], *m_geom, do_pml_in_domain); + Exchange(*pml_E_fp[2], *Ep[2], *m_geom, do_pml_in_domain); } else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0]) { - Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom); - Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom); - Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom); + Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom, do_pml_in_domain); + Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom, do_pml_in_domain); + Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom, do_pml_in_domain); } } void -PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp) +PML::CopyJtoPMLs (PatchType patch_type, + const std::array<amrex::MultiFab*,3>& jp) +{ + if (patch_type == PatchType::fine && pml_j_fp[0] && jp[0]) + { + CopyToPML(*pml_j_fp[0], *jp[0], *m_geom); + CopyToPML(*pml_j_fp[1], *jp[1], *m_geom); + CopyToPML(*pml_j_fp[2], *jp[2], *m_geom); + } + else if (patch_type == PatchType::coarse && pml_j_cp[0] && jp[0]) + { + CopyToPML(*pml_j_cp[0], *jp[0], *m_cgeom); + CopyToPML(*pml_j_cp[1], *jp[1], *m_cgeom); + CopyToPML(*pml_j_cp[2], *jp[2], *m_cgeom); + } +} + +void +PML::CopyJtoPMLs (const std::array<amrex::MultiFab*,3>& j_fp, + const std::array<amrex::MultiFab*,3>& j_cp) +{ + CopyJtoPMLs(PatchType::fine, j_fp); + CopyJtoPMLs(PatchType::coarse, j_cp); +} + + +void +PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain) { - ExchangeF(PatchType::fine, F_fp); - ExchangeF(PatchType::coarse, F_cp); + ExchangeF(PatchType::fine, F_fp, do_pml_in_domain); + ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain); } void -PML::ExchangeF (PatchType patch_type, MultiFab* Fp) +PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_F_fp && Fp) { - Exchange(*pml_F_fp, *Fp, *m_geom); + Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain); } else if (patch_type == PatchType::coarse && pml_F_cp && Fp) { - Exchange(*pml_F_cp, *Fp, *m_cgeom); + Exchange(*pml_F_cp, *Fp, *m_cgeom, do_pml_in_domain); } } + void -PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom) +PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom, + int do_pml_in_domain) { + BL_PROFILE("PML::Exchange"); + const IntVect& ngr = reg.nGrowVect(); const IntVect& ngp = pml.nGrowVect(); const int ncp = pml.nComp(); const auto& period = geom.periodicity(); + // Create temporary MultiFab to copy to and from the PML MultiFab tmpregmf(reg.boxArray(), reg.DistributionMap(), ncp, ngr); - if (ngp.max() > 0) // Copy from pml to the ghost cells of regular data - { - MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, 0); - MultiFab::LinComb(totpmlmf, 1.0, pml, 0, 1.0, pml, 1, 0, 1, 0); - if (ncp == 3) { - MultiFab::Add(totpmlmf,pml,2,0,1,0); - } - - MultiFab::Copy(tmpregmf, reg, 0, 0, 1, ngr); - tmpregmf.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), ngr, period); + // Create the sum of the split fields, in the PML + MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, 0); // Allocate + MultiFab::LinComb(totpmlmf, 1.0, pml, 0, 1.0, pml, 1, 0, 1, 0); // Sum + if (ncp == 3) { + MultiFab::Add(totpmlmf,pml,2,0,1,0); // Sum the third split component + } + // Copy from the sum of PML split field to valid cells of regular grid + if (do_pml_in_domain){ + // Valid cells of the PML and of the regular grid overlap + // Copy from valid cells of the PML to valid cells of the regular grid + reg.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), IntVect(0), period); + } else { + // Valid cells of the PML only overlap with guard cells of regular grid + // (and outermost valid cell of the regular grid, for nodal direction) + // Copy from valid cells of PML to ghost cells of regular grid + // but avoid updating the outermost valid cell + if (ngr.max() > 0) { + MultiFab::Copy(tmpregmf, reg, 0, 0, 1, ngr); + tmpregmf.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), ngr, period); #ifdef _OPENMP #pragma omp parallel #endif - for (MFIter mfi(reg); mfi.isValid(); ++mfi) - { - const FArrayBox& src = tmpregmf[mfi]; - FArrayBox& dst = reg[mfi]; - const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox()); - for (const Box& bx : bl) + for (MFIter mfi(reg); mfi.isValid(); ++mfi) { - dst.copy(src, bx, 0, bx, 0, 1); + const FArrayBox& src = tmpregmf[mfi]; + FArrayBox& dst = reg[mfi]; + const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox()); + // boxDiff avoids the outermost valid cell + for (const Box& bx : bl) { + dst.copy(src, bx, 0, bx, 0, 1); + } } } } - // Copy from regular data to PML's first component + // Copy from valid cells of the regular grid to guard cells of the PML + // (and outermost valid cell in the nodal direction) + // More specifically, copy from regular data to PML's first component // Zero out the second (and third) component - MultiFab::Copy(tmpregmf,reg,0,0,1,0); - tmpregmf.setVal(0.0, 1, ncp-1, 0); + MultiFab::Copy(tmpregmf,reg,0,0,1,0); // Fill first component of tmpregmf + tmpregmf.setVal(0.0, 1, ncp-1, 0); // Zero out the second (and third) component + if (do_pml_in_domain){ + // Where valid cells of tmpregmf overlap with PML valid cells, + // copy the PML (this is order to avoid overwriting PML valid cells, + // in the next `ParallelCopy`) + tmpregmf.ParallelCopy(pml,0, 0, ncp, IntVect(0), IntVect(0), period); + } pml.ParallelCopy(tmpregmf, 0, 0, ncp, IntVect(0), ngp, period); } + +void +PML::CopyToPML (MultiFab& pml, MultiFab& reg, const Geometry& geom) +{ + const IntVect& ngr = reg.nGrowVect(); + const IntVect& ngp = pml.nGrowVect(); + const auto& period = geom.periodicity(); + + pml.ParallelCopy(reg, 0, 0, 1, ngr, ngp, period); +} + void PML::FillBoundary () { |