diff options
Diffstat (limited to 'Source')
24 files changed, 703 insertions, 657 deletions
diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index b34cbe88b..9e04322f5 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -47,10 +47,15 @@ struct SigmaBox using SigmaVect = std::array<Sigma,AMREX_SPACEDIM>; - SigmaVect sigma; // sigma/epsilon - SigmaVect sigma_star; // sigma_star/mu + SigmaVect sigma; + SigmaVect sigma_cumsum; + SigmaVect sigma_star; + SigmaVect sigma_star_cumsum; SigmaVect sigma_fac; + SigmaVect sigma_cumsum_fac; SigmaVect sigma_star_fac; + SigmaVect sigma_star_cumsum_fac; + }; namespace amrex { @@ -102,14 +107,18 @@ public: amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, #endif int do_dive_cleaning, int do_moving_window, - const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi); + int pml_has_particles, int do_pml_in_domain, + const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), + const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); void ComputePMLFactors (amrex::Real dt); std::array<amrex::MultiFab*,3> GetE_fp (); std::array<amrex::MultiFab*,3> GetB_fp (); + std::array<amrex::MultiFab*,3> Getj_fp (); std::array<amrex::MultiFab*,3> GetE_cp (); std::array<amrex::MultiFab*,3> GetB_cp (); + std::array<amrex::MultiFab*,3> Getj_cp (); amrex::MultiFab* GetF_fp (); amrex::MultiFab* GetF_cp (); @@ -125,16 +134,21 @@ public: #endif void 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); void 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); + void CopyJtoPMLs (const std::array<amrex::MultiFab*,3>& j_fp, + const std::array<amrex::MultiFab*,3>& j_cp); + void ExchangeB (PatchType patch_type, - const std::array<amrex::MultiFab*,3>& Bp); + const std::array<amrex::MultiFab*,3>& Bp, int do_pml_in_domain); void ExchangeE (PatchType patch_type, - const std::array<amrex::MultiFab*,3>& Ep); + const std::array<amrex::MultiFab*,3>& Ep, int do_pml_in_domain); + void CopyJtoPMLs (PatchType patch_type, + const std::array<amrex::MultiFab*,3>& jp); - void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp); - void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp); + void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain); + void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain); void FillBoundary (); void FillBoundaryE (); @@ -157,9 +171,11 @@ private: std::array<std::unique_ptr<amrex::MultiFab>,3> pml_E_fp; std::array<std::unique_ptr<amrex::MultiFab>,3> pml_B_fp; + std::array<std::unique_ptr<amrex::MultiFab>,3> pml_j_fp; std::array<std::unique_ptr<amrex::MultiFab>,3> pml_E_cp; std::array<std::unique_ptr<amrex::MultiFab>,3> pml_B_cp; + std::array<std::unique_ptr<amrex::MultiFab>,3> pml_j_cp; std::unique_ptr<amrex::MultiFab> pml_F_fp; std::unique_ptr<amrex::MultiFab> pml_F_cp; @@ -173,11 +189,13 @@ private: #endif static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom, - const amrex::BoxArray& grid_ba, int ncell, - const amrex::IntVect do_pml_Lo, - const amrex::IntVect do_pml_Hi); + const amrex::BoxArray& grid_ba, + int ncell, int do_pml_in_domain, + const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), + const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); - static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); + static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom, int do_pml_in_domain); + static void CopyToPML (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); }; #ifdef WARPX_USE_PSATD 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 () { diff --git a/Source/BoundaryConditions/PML_current.H b/Source/BoundaryConditions/PML_current.H new file mode 100644 index 000000000..910186a96 --- /dev/null +++ b/Source/BoundaryConditions/PML_current.H @@ -0,0 +1,129 @@ +#ifndef PML_CURRENT_H_ +#define PML_CURRENT_H_ + +#include <AMReX_FArrayBox.H> + +using namespace amrex; + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void push_ex_pml_current (int j, int k, int l, Array4<Real> const& Ex, + Array4<Real const> const& jx, + Real const * const sigjy, + Real const * const sigjz, + int ylo, int zlo, + Real mu_c2_dt) +{ +#if (AMREX_SPACEDIM == 3) + Real alpha_xy, alpha_xz; + if (sigjy[k-ylo]+sigjz[l-zlo] == 0){ + alpha_xy = 0.5; + alpha_xz = 0.5; + } + else { + alpha_xy = sigjy[k-ylo]/(sigjy[k-ylo]+sigjz[l-zlo]); + alpha_xz = sigjz[l-zlo]/(sigjy[k-ylo]+sigjz[l-zlo]); + } + Ex(j,k,l,0) = Ex(j,k,l,0) - mu_c2_dt * alpha_xy * jx(j,k,l); + Ex(j,k,l,1) = Ex(j,k,l,1) - mu_c2_dt * alpha_xz * jx(j,k,l); +#else + Ex(j,k,l,1) = Ex(j,k,l,1) - mu_c2_dt * jx(j,k,l); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void push_ey_pml_current (int j, int k, int l, Array4<Real> const& Ey, + Array4<Real const> const& jy, + Real const * const sigjx, + Real const * const sigjz, + int xlo, int zlo, + Real mu_c2_dt) +{ +#if (AMREX_SPACEDIM == 3) + Real alpha_yx, alpha_yz; + if (sigjx[j-xlo]+sigjz[l-zlo] == 0){ + alpha_yx = 0.5; + alpha_yz = 0.5; + } + else { + alpha_yx = sigjx[j-xlo]/(sigjx[j-xlo]+sigjz[l-zlo]); + alpha_yz = sigjz[l-zlo]/(sigjx[j-xlo]+sigjz[l-zlo]); + } + Ey(j,k,l,0) = Ey(j,k,l,0) - mu_c2_dt * alpha_yx * jy(j,k,l); + Ey(j,k,l,1) = Ey(j,k,l,1) - mu_c2_dt * alpha_yz * jy(j,k,l); +#else + Ey(j,k,l,0) = Ey(j,k,l,0) - 0.5 * mu_c2_dt * jy(j,k,l); + Ey(j,k,l,1) = Ey(j,k,l,1) - 0.5 * mu_c2_dt * jy(j,k,l); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void push_ez_pml_current (int j, int k, int l, Array4<Real> const& Ez, + Array4<Real const> const& jz, + Real const * const sigjx, + Real const * const sigjy, + int xlo, int ylo, + Real mu_c2_dt) +{ +#if (AMREX_SPACEDIM == 3) + Real alpha_zx, alpha_zy; + if (sigjx[j-xlo]+sigjy[k-ylo]==0){ + alpha_zx = 0.5; + alpha_zy = 0.5; + } + else { + alpha_zx = sigjx[j-xlo]/(sigjx[j-xlo]+sigjy[k-ylo]); + alpha_zy = sigjy[k-ylo]/(sigjx[j-xlo]+sigjy[k-ylo]); + } + Ez(j,k,l,0) = Ez(j,k,l,0) - mu_c2_dt * alpha_zx * jz(j,k,l); + Ez(j,k,l,1) = Ez(j,k,l,1) - mu_c2_dt * alpha_zy * jz(j,k,l); +#else + Ez(j,k,l,0) = Ez(j,k,l,0) - mu_c2_dt * jz(j,k,l); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void damp_jx_pml (int j, int k, int l, + Array4<Real> const& jx, + Real const* const sigsjx, + Real const* const sigjy, + Real const* const sigjz, + int xlo, int ylo, int zlo) +{ +#if (AMREX_SPACEDIM == 3) + jx(j,k,l) = jx(j,k,l) * sigsjx[j-xlo] * sigjy[k-ylo] * sigjz[l-zlo]; +#else + jx(j,k,l) = jx(j,k,l) * sigsjx[j-xlo] * sigjz[k-zlo]; +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void damp_jy_pml (int j, int k, int l, + Array4<Real> const& jy, + Real const * const sigjx, + Real const * const sigsjy, + Real const * const sigjz, + int xlo, int ylo, int zlo) +{ +#if (AMREX_SPACEDIM == 3) + jy(j,k,l) = jy(j,k,l) * sigjx[j-xlo] * sigsjy[k-ylo] * sigjz[l-zlo]; +#else + jy(j,k,l) = jy(j,k,l) * sigjx[j-xlo] * sigjz[k-zlo]; +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void damp_jz_pml (int j, int k, int l, + Array4<Real> const& jz, + Real const * const sigjx, + Real const * const sigjy, + Real const * const sigsjz, + int xlo, int ylo, int zlo) +{ +#if (AMREX_SPACEDIM == 3) + jz(j,k,l) = jz(j,k,l) * sigjx[j-xlo] * sigjy[k-ylo] * sigsjz[l-zlo]; +#else + jz(j,k,l) = jz(j,k,l) * sigjx[j-xlo] * sigsjz[k-zlo]; +#endif +} + +#endif diff --git a/Source/BoundaryConditions/PML_routines.F90 b/Source/BoundaryConditions/PML_routines.F90 index 380e52934..e2540350e 100644 --- a/Source/BoundaryConditions/PML_routines.F90 +++ b/Source/BoundaryConditions/PML_routines.F90 @@ -1006,4 +1006,5 @@ contains end do end subroutine warpx_damp_pml_f_3d + end module warpx_pml_module 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() }); + } + } +} diff --git a/Source/Diagnostics/SliceDiagnostic.cpp b/Source/Diagnostics/SliceDiagnostic.cpp index 9f365b39d..0ec528e32 100644 --- a/Source/Diagnostics/SliceDiagnostic.cpp +++ b/Source/Diagnostics/SliceDiagnostic.cpp @@ -140,10 +140,10 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, MFIter mfi_dst(mfDst); for (MFIter mfi(mfSrc); mfi.isValid(); ++mfi) { - FArrayBox& Src_fabox = mfSrc[mfi]; + Array4<Real const> const& Src_fabox = mfSrc.const_array(mfi); const Box& Dst_bx = mfi_dst.validbox(); - FArrayBox& Dst_fabox = mfDst[mfi_dst]; + Array4<Real> const& Dst_fabox = mfDst.array(mfi_dst); int scomp = 0; int dcomp = 0; diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index b9aa3b5e8..8f800665d 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -13,6 +13,7 @@ #include <AMReX_AmrMeshInSituBridge.H> #endif + using namespace amrex; void @@ -135,7 +136,7 @@ WarpX::EvolveEM (int numsteps) bool to_make_plot = (plot_int > 0) && ((step+1) % plot_int == 0); // slice generation // - bool to_make_slice_plot = (slice_plot_int > 0) && ( (step+1)% slice_plot_int == 0); + bool to_make_slice_plot = (slice_plot_int > 0) && ( (step+1)% slice_plot_int == 0); bool do_insitu = ((step+1) >= insitu_start) && (insitu_int > 0) && ((step+1) % insitu_int == 0); @@ -153,7 +154,7 @@ WarpX::EvolveEM (int numsteps) // We might need to move j because we are going to make a plotfile. int num_moved = MoveWindow(move_j); - + if (max_level == 0) { int num_redistribute_ghost = num_moved + 1; mypc->RedistributeLocal(num_redistribute_ghost); @@ -228,7 +229,7 @@ WarpX::EvolveEM (int numsteps) // End loop on time steps } - bool write_plot_file = plot_int > 0 && istep[0] > last_plot_file_step + bool write_plot_file = plot_int > 0 && istep[0] > last_plot_file_step && (max_time_reached || istep[0] >= max_step); bool do_insitu = (insitu_start >= istep[0]) && (insitu_int > 0) && @@ -255,7 +256,7 @@ WarpX::EvolveEM (int numsteps) UpdateInSitu(); } - if (check_int > 0 && istep[0] > last_check_file_step && + if (check_int > 0 && istep[0] > last_check_file_step && (max_time_reached || istep[0] >= max_step)) { WriteCheckPointFile(); } @@ -298,6 +299,10 @@ WarpX::OneStep_nosub (Real cur_time) SyncRho(); + // For extended PML: copy J from regular grid to PML, and damp J in PML + if (do_pml && pml_has_particles) CopyJPML(); + if (do_pml && do_pml_j_damping) DampJPML(); + // Push E and B from {n} to {n+1} // (And update guard cells immediately afterwards) #ifdef WARPX_USE_PSATD @@ -310,6 +315,7 @@ WarpX::OneStep_nosub (Real cur_time) FillBoundaryF(); EvolveB(0.5*dt[0]); // We now have B^{n+1/2} FillBoundaryB(); + EvolveE(dt[0]); // We now have E^{n+1} FillBoundaryE(); EvolveF(0.5*dt[0], DtType::SecondHalf); @@ -319,6 +325,7 @@ WarpX::OneStep_nosub (Real cur_time) FillBoundaryE(); } FillBoundaryB(); + #endif } @@ -559,13 +566,13 @@ WarpX::ComputeDt () /* \brief computes max_step for wakefield simulation in boosted frame. * \param geom: Geometry object that contains simulation domain. - * - * max_step is set so that the simulation stop when the lower corner of the + * + * max_step is set so that the simulation stop when the lower corner of the * simulation box passes input parameter zmax_plasma_to_compute_max_step. */ void WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ - // Sanity checks: can use zmax_plasma_to_compute_max_step only if + // Sanity checks: can use zmax_plasma_to_compute_max_step only if // the moving window and the boost are all in z direction. AMREX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::moving_window_dir == AMREX_SPACEDIM-1, @@ -579,7 +586,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ "Can use zmax_plasma_to_compute_max_step only if " + "warpx.boost_direction = z. TODO: all directions."); - // Lower end of the simulation domain. All quantities are given in boosted + // Lower end of the simulation domain. All quantities are given in boosted // frame except zmax_plasma_to_compute_max_step. const Real zmin_domain_boost = a_geom.ProbLo(AMREX_SPACEDIM-1); // End of the plasma: Transform input argument @@ -587,7 +594,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ const Real len_plasma_boost = zmax_plasma_to_compute_max_step/gamma_boost; // Plasma velocity const Real v_plasma_boost = -beta_boost * PhysConst::c; - // Get time at which the lower end of the simulation domain passes the + // Get time at which the lower end of the simulation domain passes the // upper end of the plasma (in the z direction). const Real interaction_time_boost = (len_plasma_boost-zmin_domain_boost)/ (moving_window_v-v_plasma_boost); @@ -605,7 +612,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ /* \brief Apply perfect mirror condition inside the box (not at a boundary). * In practice, set all fields to 0 on a section of the simulation domain - * (as for a perfect conductor with a given thickness). + * (as for a perfect conductor with a given thickness). * The mirror normal direction has to be parallel to the z axis. */ void @@ -622,10 +629,10 @@ WarpX::applyMirrors(Real time){ } // Loop over levels for(int lev=0; lev<=finest_level; lev++){ - // Make sure that the mirror contains at least + // Make sure that the mirror contains at least // mirror_z_npoints[i_mirror] cells Real dz = WarpX::CellSize(lev)[2]; - Real z_max = std::max(z_max_tmp, + Real z_max = std::max(z_max_tmp, z_min+mirror_z_npoints[i_mirror]*dz); // Get fine patch field MultiFabs MultiFab& Ex = *Efield_fp[lev][0].get(); diff --git a/Source/FieldSolver/Make.package b/Source/FieldSolver/Make.package index 25e21c146..076be41bb 100644 --- a/Source/FieldSolver/Make.package +++ b/Source/FieldSolver/Make.package @@ -3,7 +3,7 @@ CEXE_headers += WarpX_FDTD.H CEXE_sources += WarpXPushFieldsEM.cpp ifeq ($(USE_PSATD),TRUE) include $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/Make.package - ifeq ($(USE_CUDA),FALSE) + ifeq ($(USE_PSATD_PICSAR),TRUE) include $(WARPX_HOME)/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package endif endif diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index b876f2751..b3d8e9d76 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -11,6 +11,8 @@ #include <WarpX_py.H> #endif +#include <PML_current.H> + #ifdef BL_USE_SENSEI_INSITU #include <AMReX_AmrMeshInSituBridge.H> #endif @@ -59,7 +61,7 @@ WarpX::PushPSATD (amrex::Real a_dt) for (int lev = 0; lev <= finest_level; ++lev) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dt[lev] == a_dt, "dt must be consistent"); if (fft_hybrid_mpi_decomposition){ -#ifndef AMREX_USE_CUDA // Only available on CPU +#ifdef WARPX_USE_PSATD_HYBRID PushPSATD_hybridFFT(lev, a_dt); #endif } else { @@ -136,7 +138,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) MultiFab* cost = costs[lev].get(); const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); - // xmin is only used by the picsar kernel with cylindrical geometry, + // xmin is only used by the kernel for cylindrical geometry, // in which case it is actually rmin. const Real xmin = Geom(0).ProbLo(0); @@ -324,7 +326,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) MultiFab* cost = costs[lev].get(); const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); - // xmin is only used by the picsar kernel with cylindrical geometry, + // xmin is only used by the kernel for cylindrical geometry, // in which case it is actually rmin. const Real xmin = Geom(0).ProbLo(0); @@ -447,11 +449,14 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) if (do_pml && pml[lev]->ok()) { - if (F) pml[lev]->ExchangeF(patch_type, F); + if (F) pml[lev]->ExchangeF(patch_type, F, do_pml_in_domain); const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_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 @@ -461,6 +466,10 @@ 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& pml_Exfab = pml_E[0]->array(mfi); + auto const& pml_Eyfab = pml_E[1]->array(mfi); + auto const& pml_Ezfab = pml_E[2]->array(mfi); + WRPX_PUSH_PML_EVEC( tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), @@ -471,7 +480,45 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), - &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); + &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); + + if (pml_has_particles) { + // Update the E field in the PML, using the current + // deposited by the particles in the PML + 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* sigmaj_x = sigba[mfi].sigma[0].data(); + const Real* sigmaj_y = sigba[mfi].sigma[1].data(); + const Real* sigmaj_z = sigba[mfi].sigma[2].data(); + + auto const& AMREX_RESTRICT x_lo = sigba[mfi].sigma[0].lo(); +#if (AMREX_SPACEDIM == 3) + auto const& AMREX_RESTRICT y_lo = sigba[mfi].sigma[1].lo(); + auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma[2].lo(); +#else + int y_lo = 0; + auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma[1].lo(); +#endif + + amrex::ParallelFor( tex, tey, tez, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + push_ex_pml_current(i,j,k, + pml_Exfab, pml_jxfab, sigmaj_y, sigmaj_z, + y_lo, z_lo, mu_c2_dt); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + push_ey_pml_current(i,j,k, + pml_Eyfab, pml_jyfab, sigmaj_x, sigmaj_z, + x_lo, z_lo, mu_c2_dt); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + push_ez_pml_current(i,j,k, + pml_Ezfab, pml_jzfab, sigmaj_x, sigmaj_y, + x_lo, y_lo, mu_c2_dt); + } + ); + } if (pml_F) { @@ -483,7 +530,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), BL_TO_FORTRAN_3D((*pml_F )[mfi]), - &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, + &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, &WarpX::maxwell_fdtd_solver_id); } } diff --git a/Source/FortranInterface/Make.package b/Source/FortranInterface/Make.package index fe0c2ba3d..d7b17fa56 100644 --- a/Source/FortranInterface/Make.package +++ b/Source/FortranInterface/Make.package @@ -1,6 +1,5 @@ CEXE_headers += WarpX_f.H F90EXE_sources += WarpX_f.F90 -F90EXE_sources += WarpX_picsar.F90 INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FortranInterface VPATH_LOCATIONS += $(WARPX_HOME)/Source/FortranInterface diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 6dff3469d..26da42606 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -1,4 +1,3 @@ - #include <AMReX_BLFort.H> #ifdef __cplusplus @@ -17,8 +16,6 @@ #if (AMREX_SPACEDIM == 3) -#define WRPX_COMPUTE_DIVB warpx_compute_divb_3d -#define WRPX_COMPUTE_DIVE warpx_compute_dive_3d #define WRPX_SYNC_CURRENT warpx_sync_current_3d #define WRPX_SYNC_RHO warpx_sync_rho_3d @@ -28,7 +25,6 @@ #define WRPX_PUSH_PML_F warpx_push_pml_f_3d #define WRPX_DAMP_PML warpx_damp_pml_3d #define WRPX_DAMP_PML_F warpx_damp_pml_f_3d - #define WRPX_SUM_FINE_TO_CRSE_NODAL warpx_sum_fine_to_crse_nodal_3d #define WRPX_ZERO_OUT_BNDRY warpx_zero_out_bndry_3d #define WRPX_BUILD_MASK warpx_build_mask_3d @@ -41,7 +37,6 @@ #elif (AMREX_SPACEDIM == 2) -#define WRPX_COMPUTE_DIVB warpx_compute_divb_2d #define WRPX_SYNC_CURRENT warpx_sync_current_2d #define WRPX_SYNC_RHO warpx_sync_rho_2d @@ -62,12 +57,6 @@ #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_2d #define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_2d -#ifdef WARPX_DIM_RZ -#define WRPX_COMPUTE_DIVE warpx_compute_dive_rz -#else -#define WRPX_COMPUTE_DIVE warpx_compute_dive_2d -#endif - #endif #ifdef __cplusplus @@ -75,45 +64,7 @@ extern "C" { #endif - // Current deposition - void warpx_current_deposition( - amrex::Real* jx, const long* jx_ng, const int* jx_ntot, - amrex::Real* jy, const long* jy_ng, const int* jy_ntot, - amrex::Real* jz, const long* jz_ng, const int* jz_ntot, - const long* n_rz_azimuthal_modes, - const long* np, - const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp, - const amrex::Real* uxp, const amrex::Real* uyp,const amrex::Real* uzp, - const amrex::Real* gip, const amrex::Real* w, const amrex::Real* q, - const amrex::Real* xmin, const amrex::Real* ymin, const amrex::Real* zmin, - const amrex::Real* dt, - const amrex::Real* dx, const amrex::Real* dy, const amrex::Real* dz, - const long* nox, const long* noy,const long* noz, - const int* l_nodal, const long* lvect, const long* current_depo_algo); - - // Particle pusher (velocity and position) - - void warpx_particle_pusher(const long* np, - amrex::Real* xp, amrex::Real* yp, amrex::Real* zp, - amrex::Real* uxp, amrex::Real* uyp, amrex::Real* uzp, amrex::Real* gaminv, - const amrex::Real* exp, const amrex::Real* eyp,const amrex::Real* ezp, - const amrex::Real* bxp, const amrex::Real* byp,const amrex::Real* bzp, - const amrex::Real* charge, const amrex::Real* mass, const amrex::Real* dt, - const long* particle_pusher_algo); - - // Particle pusher (velocity) - - void warpx_particle_pusher_momenta(const long* np, - amrex::Real* xp, amrex::Real* yp, amrex::Real* zp, - amrex::Real* uxp, amrex::Real* uyp, amrex::Real* uzp, amrex::Real* gaminv, - const amrex::Real* exp, const amrex::Real* eyp,const amrex::Real* ezp, - const amrex::Real* bxp, const amrex::Real* byp,const amrex::Real* bzp, - const amrex::Real* charge, const amrex::Real* mass, const amrex::Real* dt, - const long* particle_pusher_algo); - - // Laser pusher - void warpx_gaussian_laser( const long* np, amrex::Real* Xp, amrex::Real* Yp, amrex::Real* t, amrex::Real* wavelength, amrex::Real* e_max, amrex::Real* waist, amrex::Real* duration, amrex::Real* t_peak, amrex::Real* f, amrex::Real* amplitude, @@ -139,60 +90,6 @@ extern "C" amrex::Real* mobility, amrex::Real* dt, const amrex::Real* c, const amrex::Real* beta_boost, const amrex::Real* gamma_boost ); - // Maxwell solver - - void warpx_push_evec( - const int* xlo, const int* xhi, - const int* ylo, const int* yhi, - const int* zlo, const int* zhi, - const long* n_rz_azimuthal_modes, - BL_FORT_FAB_ARG_3D(ex), - BL_FORT_FAB_ARG_3D(ey), - BL_FORT_FAB_ARG_3D(ez), - const BL_FORT_FAB_ARG_3D(bx), - const BL_FORT_FAB_ARG_3D(by), - const BL_FORT_FAB_ARG_3D(bz), - const BL_FORT_FAB_ARG_3D(jx), - const BL_FORT_FAB_ARG_3D(jy), - const BL_FORT_FAB_ARG_3D(jz), - const amrex::Real* mudt, - const amrex::Real* dtsdx, - const amrex::Real* dtsdy, - const amrex::Real* dtsdz, - const amrex::Real* xmin, - const amrex::Real* dx); - - void warpx_push_bvec( - const int* xlo, const int* xhi, - const int* ylo, const int* yhi, - const int* zlo, const int* zhi, - const long* n_rz_azimuthal_modes, - const BL_FORT_FAB_ARG_3D(ex), - const BL_FORT_FAB_ARG_3D(ey), - const BL_FORT_FAB_ARG_3D(ez), - BL_FORT_FAB_ARG_3D(bx), - BL_FORT_FAB_ARG_3D(by), - BL_FORT_FAB_ARG_3D(bz), - const amrex::Real* dtsdx, - const amrex::Real* dtsdy, - const amrex::Real* dtsdz, - const amrex::Real* xmin, - const amrex::Real* dx, - const int* maxwell_fdtd_solver_id); - - void warpx_push_evec_f( - const int* xlo, const int* xhi, - const int* ylo, const int* yhi, - const int* zlo, const int* zhi, - BL_FORT_FAB_ARG_3D(ex), - BL_FORT_FAB_ARG_3D(ey), - BL_FORT_FAB_ARG_3D(ez), - const BL_FORT_FAB_ARG_3D(f), - const amrex::Real* dtsdx_c2, - const amrex::Real* dtsdy_c2, - const amrex::Real* dtsdz_c2, - const int* maxwell_fdtd_solver_id); - #ifdef USE_OPENBC_POISSON void warpx_openbc_potential (amrex::Real* rho, amrex::Real* phi, const amrex::Real* dx); void warpx_openbc_decompose (const int*, const int*, int*, int*); @@ -288,24 +185,6 @@ extern "C" // These functions are used to evolve E and B in the PML - void WRPX_COMPUTE_DIVB (const int* lo, const int* hi, - BL_FORT_FAB_ARG_ANYD(divb), - const BL_FORT_FAB_ARG_ANYD(Bx), - const BL_FORT_FAB_ARG_ANYD(By), - const BL_FORT_FAB_ARG_ANYD(Bz), - const amrex::Real* dx); - - void WRPX_COMPUTE_DIVE (const int* lo, const int* hi, - BL_FORT_FAB_ARG_ANYD(dive), - const BL_FORT_FAB_ARG_ANYD(ex), - const BL_FORT_FAB_ARG_ANYD(ey), - const BL_FORT_FAB_ARG_ANYD(ez), - const amrex::Real* dx -#ifdef WARPX_DIM_RZ - ,const amrex::Real* rmin -#endif - ); - void WRPX_PUSH_PML_BVEC(const int* xlo, const int* xhi, const int* ylo, const int* yhi, const int* zlo, const int* zhi, @@ -391,6 +270,7 @@ extern "C" #endif const amrex::Real* sigbz, int sigbz_lo, int sigbz_hi); + void WRPX_SYNC_CURRENT (const int* lo, const int* hi, BL_FORT_FAB_ARG_ANYD(crse), const BL_FORT_FAB_ARG_ANYD(fine), diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 deleted file mode 100644 index 3a9f8f41e..000000000 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ /dev/null @@ -1,187 +0,0 @@ -#if (AMREX_SPACEDIM == 3) - -#define WRPX_PXR_CURRENT_DEPOSITION depose_jxjyjz_generic - -#elif (AMREX_SPACEDIM == 2) - -#ifdef WARPX_DIM_RZ - -#define WRPX_PXR_CURRENT_DEPOSITION depose_jrjtjz_generic_rz - -#else - -#define WRPX_PXR_CURRENT_DEPOSITION depose_jxjyjz_generic_2d - -#endif - -#endif - -#define LVECT_CURRDEPO 8_c_long -#define LVECT_FIELDGATHE 64_c_long - -! _________________________________________________________________ -! -!> @brief -!> Module that contains subroutines to be called with AMReX -!> and that uses subroutines of Picsar -!> -!> @details -!> This avoids the use of interface with bind in the core of Picsar -!> This enables the use of integer in AMReX and Logical in Picsar -!> wihtout compatibility issue -!> -!> @author -!> Weiqun Zhang -!> Ann Almgren -!> Remi Lehe -!> Mathieu Lobet -!> -module warpx_to_pxr_module -! _________________________________________________________________ - - use iso_c_binding - use amrex_fort_module, only : amrex_real - - implicit none - - integer, parameter :: pxr_logical = 8 - -contains - - ! _________________________________________________________________ - !> - !> @brief - !> Main subroutine for the current deposition - !> - !> @details - !> This subroutines enable to controle the interpolation order - !> via the parameters nox,noy,noz and the type of algorithm via - !> the parameter current_depo_algo - ! - !> @param[inout] jx,jy,jz current arrays - !> @param[in] np number of particles - !> @param[in] xp,yp,zp particle position arrays - !> @param[in] uxp,uyp,uzp particle momentum arrays - !> @param[in] gaminv inverve of the particle Lorentz factor (array) - !> @param[in] w particle weight arrays - !> @param[in] q particle species charge - !> @param[in] xmin,ymin,zmin tile grid minimum position - !> @param[in] dx,dy,dz space discretization steps - !> @param[in] nx,ny,nz number of cells - !> @param[in] nxguard,nyguard,nzguard number of guard cells - !> @param[in] nox,noy,noz interpolation order - !> @param[in] lvect vector length - !> @param[in] charge_depo_algo algorithm choice for the charge deposition - !> - subroutine warpx_current_deposition( & - jx,jx_ng,jx_ntot,jy,jy_ng,jy_ntot,jz,jz_ng,jz_ntot,n_rz_azimuthal_modes, & - np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin,dt,dx,dy,dz,nox,noy,noz,& - l_nodal,lvect,current_depo_algo) & - bind(C, name="warpx_current_deposition") - - integer, intent(in) :: jx_ntot(AMREX_SPACEDIM), jy_ntot(AMREX_SPACEDIM), jz_ntot(AMREX_SPACEDIM) - integer(c_long), intent(in) :: jx_ng, jy_ng, jz_ng - integer(c_long), intent(IN) :: n_rz_azimuthal_modes - integer(c_long), intent(IN) :: np - integer(c_long), intent(IN) :: nox,noy,noz - integer(c_int), intent(in) :: l_nodal - -#ifdef WARPX_RZ - real(amrex_real), intent(IN OUT):: jx(jx_ntot(1),jx_ntot(2),2,n_rz_azimuthal_modes) - real(amrex_real), intent(IN OUT):: jy(jy_ntot(1),jy_ntot(2),2,n_rz_azimuthal_modes) - real(amrex_real), intent(IN OUT):: jz(jz_ntot(1),jz_ntot(2),2,n_rz_azimuthal_modes) -#else - real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*) -#endif - real(amrex_real), intent(IN) :: q - real(amrex_real), intent(IN) :: dx,dy,dz - real(amrex_real), intent(IN) :: dt - real(amrex_real), intent(IN) :: xmin,ymin,zmin - real(amrex_real), intent(IN), dimension(np) :: xp,yp,zp,w - real(amrex_real), intent(IN), dimension(np) :: uxp,uyp,uzp - real(amrex_real), intent(IN), dimension(np) :: gaminv - integer(c_long), intent(IN) :: lvect - integer(c_long), intent(IN) :: current_depo_algo - logical(pxr_logical) :: pxr_l_nodal - -#ifdef WARPX_RZ - logical(pxr_logical) :: l_particles_weight = .true. - integer(c_long) :: type_rz_depose = 1 - complex(amrex_real), allocatable, dimension(:,:,:) :: jr_c, jt_c, jz_c - integer :: alloc_status -#endif - - ! Compute the number of valid cells and guard cells - integer(c_long) :: jx_nvalid(AMREX_SPACEDIM), jy_nvalid(AMREX_SPACEDIM), jz_nvalid(AMREX_SPACEDIM), & - jx_nguards(AMREX_SPACEDIM), jy_nguards(AMREX_SPACEDIM), jz_nguards(AMREX_SPACEDIM) - jx_nvalid = jx_ntot - 2*jx_ng - jy_nvalid = jy_ntot - 2*jy_ng - jz_nvalid = jz_ntot - 2*jz_ng - jx_nguards = jx_ng - jy_nguards = jy_ng - jz_nguards = jz_ng - pxr_l_nodal = l_nodal .eq. 1 - -! Dimension 3 -#if (AMREX_SPACEDIM==3) - CALL WRPX_PXR_CURRENT_DEPOSITION( & - jx,jx_nguards,jx_nvalid, & - jy,jy_nguards,jy_nvalid, & - jz,jz_nguards,jz_nvalid, & - np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q, & - xmin,ymin,zmin,dt,dx,dy,dz, & - nox,noy,noz,pxr_l_nodal,current_depo_algo) -! Dimension 2 -#elif (AMREX_SPACEDIM==2) -#ifdef WARPX_RZ - if (n_rz_azimuthal_modes > 1) then - - allocate(jr_c(jx_ntot(1),jx_ntot(2),n_rz_azimuthal_modes), & - jt_c(jy_ntot(1),jy_ntot(2),n_rz_azimuthal_modes), & - jz_c(jz_ntot(1),jz_ntot(2),n_rz_azimuthal_modes), stat=alloc_status) - if (alloc_status /= 0) then - print*,"Error: warpx_current_deposition: complex arrays could not be allocated" - stop - endif - - jr_c = 0._amrex_real - jt_c = 0._amrex_real - jz_c = 0._amrex_real - - CALL pxr_depose_jrjtjz_esirkepov_n_2d_circ( & - jr_c,jx_nguards,jx_nvalid, & - jt_c,jy_nguards,jy_nvalid, & - jz_c,jz_nguards,jz_nvalid, & - n_rz_azimuthal_modes, & - np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q, & - xmin,zmin,dt,dx,dz, & - nox,noz,l_particles_weight,type_rz_depose) - - jx(:,:,1,:) = jx(:,:,1,:) + real(jr_c) - jx(:,:,2,:) = jx(:,:,2,:) + aimag(jr_c) - jy(:,:,1,:) = jy(:,:,1,:) + real(jt_c) - jy(:,:,2,:) = jy(:,:,2,:) + aimag(jt_c) - jz(:,:,1,:) = jz(:,:,1,:) + real(jz_c) - jz(:,:,2,:) = jz(:,:,2,:) + aimag(jz_c) - - deallocate(jr_c) - deallocate(jt_c) - deallocate(jz_c) - - else -#endif - CALL WRPX_PXR_CURRENT_DEPOSITION( & - jx,jx_nguards,jx_nvalid, & - jy,jy_nguards,jy_nvalid, & - jz,jz_nguards,jz_nvalid, & - np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q, & - xmin,zmin,dt,dx,dz,nox,noz,pxr_l_nodal, & - lvect,current_depo_algo) -#ifdef WARPX_RZ - endif -#endif -#endif - - end subroutine warpx_current_deposition - -end module warpx_to_pxr_module diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 590c11b84..4c0dfee17 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -142,6 +142,7 @@ WarpX::InitPML () dt[0], nox_fft, noy_fft, noz_fft, do_nodal, #endif do_dive_cleaning, do_moving_window, + pml_has_particles, do_pml_in_domain, do_pml_Lo_corrected, do_pml_Hi)); for (int lev = 1; lev <= finest_level; ++lev) { @@ -159,6 +160,7 @@ WarpX::InitPML () dt[lev], nox_fft, noy_fft, noz_fft, do_nodal, #endif do_dive_cleaning, do_moving_window, + pml_has_particles, do_pml_in_domain, do_pml_Lo_MR, amrex::IntVect::TheUnitVector())); } } diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index c141ae4b0..3fdd8730b 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -397,8 +397,8 @@ LaserParticleContainer::Evolve (int lev, { BL_PROFILE("Laser::Evolve()"); BL_PROFILE_VAR_NS("Laser::Evolve::Copy", blp_copy); - BL_PROFILE_VAR_NS("PICSAR::LaserParticlePush", blp_pxr_pp); - BL_PROFILE_VAR_NS("PICSAR::LaserCurrentDepo", blp_pxr_cd); + BL_PROFILE_VAR_NS("Laser::ParticlePush", blp_pp); + BL_PROFILE_VAR_NS("Laser::CurrentDepo", blp_cd); BL_PROFILE_VAR_NS("Laser::Evolve::Accumulate", blp_accumulate); Real t_lab = t; @@ -466,7 +466,7 @@ LaserParticleContainer::Evolve (int lev, // // Particle Push // - BL_PROFILE_VAR_START(blp_pxr_pp); + BL_PROFILE_VAR_START(blp_pp); // Find the coordinates of the particles in the emission plane calculate_laser_plane_coordinates( &np, m_xp[thread_num].dataPtr(), @@ -506,7 +506,7 @@ LaserParticleContainer::Evolve (int lev, wp.dataPtr(), amplitude_E.dataPtr(), &p_X[0], &p_X[1], &p_X[2], &nvec[0], &nvec[1], &nvec[2], &mobility, &dt, &PhysConst::c, &WarpX::beta_boost, &WarpX::gamma_boost ); - BL_PROFILE_VAR_STOP(blp_pxr_pp); + BL_PROFILE_VAR_STOP(blp_pp); // // Current Deposition diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 1970c8fc6..ea06f8e06 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -116,7 +116,7 @@ ifeq ($(USE_OPENPMD), TRUE) endif DEFINES += -DWARPX_USE_OPENPMD endif - + ifeq ($(USE_PSATD),TRUE) USERSuffix := $(USERSuffix).PSATD @@ -124,15 +124,17 @@ ifeq ($(USE_PSATD),TRUE) ifeq ($(USE_CUDA),FALSE) # Running on CPU # Use FFTW libraries += -lfftw3_mpi -lfftw3 -lfftw3_threads - FFTW_HOME ?= NOT_SET - ifneq ($(FFTW_HOME),NOT_SET) - VPATH_LOCATIONS += $(FFTW_HOME)/include - INCLUDE_LOCATIONS += $(FFTW_HOME)/include - LIBRARY_LOCATIONS += $(FFTW_HOME)/lib + ifeq ($(USE_PSATD_PICSAR),TRUE) + FFTW_HOME ?= NOT_SET + ifneq ($(FFTW_HOME),NOT_SET) + VPATH_LOCATIONS += $(FFTW_HOME)/include + INCLUDE_LOCATIONS += $(FFTW_HOME)/include + LIBRARY_LOCATIONS += $(FFTW_HOME)/lib + endif + # Compile with PICSAR's Hybrid-MPI PSATD + DEFINES += -DWARPX_USE_PSATD_HYBRID + DEFINES += -DFFTW # PICSAR uses it endif - # Compile with PICSAR's Hybrid-MPI PSATD - DEFINES += -DWARPX_USE_PSATD_HYBRID - DEFINES += -DFFTW # PICSAR uses it else # Use cuFFT libraries += -lcufft @@ -159,7 +161,7 @@ ifeq ($(USE_HDF5),TRUE) endif DEFINES += -DWARPX_USE_HDF5 libraries += -lhdf5 -lz -endif +endif # job_info support CEXE_sources += AMReX_buildInfo.cpp @@ -206,7 +208,7 @@ else ifdef WarpxBinDir -all: $(executable) +all: $(executable) $(SILENT) $(RM) AMReX_buildInfo.cpp @if [ ! -d $(WarpxBinDir) ]; then mkdir -p $(WarpxBinDir); fi $(SILENT) mv -f $(executable) $(WarpxBinDir)/ @@ -215,7 +217,7 @@ all: $(executable) else -all: $(executable) +all: $(executable) $(SILENT) $(RM) AMReX_buildInfo.cpp @echo SUCCESS @@ -223,7 +225,7 @@ endif endif -AMReX_buildInfo.cpp: +AMReX_buildInfo.cpp: $(AMREX_HOME)/Tools/C_scripts/makebuildinfo_C.py \ --amrex_home "$(AMREX_HOME)" \ --COMP "$(COMP)" --COMP_VERSION "$(COMP_VERSION)" \ diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 0ca1e8a5d..88adbc147 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -18,7 +18,8 @@ WarpX::ExchangeWithPmlB (int lev) Bfield_fp[lev][2].get() }, { Bfield_cp[lev][0].get(), Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get() }); + Bfield_cp[lev][2].get() }, + do_pml_in_domain); } } @@ -31,7 +32,8 @@ WarpX::ExchangeWithPmlE (int lev) Efield_fp[lev][2].get() }, { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), - Efield_cp[lev][2].get() }); + Efield_cp[lev][2].get() }, + do_pml_in_domain); } } @@ -40,7 +42,8 @@ WarpX::ExchangeWithPmlF (int lev) { if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeF(F_fp[lev].get(), - F_cp[lev].get()); + F_cp[lev].get(), + do_pml_in_domain); } } @@ -250,9 +253,10 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeE(patch_type, - { Efield_fp[lev][0].get(), + { Efield_fp[lev][0].get(), Efield_fp[lev][1].get(), - Efield_fp[lev][2].get() }); + Efield_fp[lev][2].get() }, + do_pml_in_domain); pml[lev]->FillBoundaryE(patch_type); } @@ -265,9 +269,10 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeE(patch_type, - { Efield_cp[lev][0].get(), + { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), - Efield_cp[lev][2].get() }); + Efield_cp[lev][2].get() }, + do_pml_in_domain); pml[lev]->FillBoundaryE(patch_type); } @@ -292,9 +297,10 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeB(patch_type, - { Bfield_fp[lev][0].get(), + { Bfield_fp[lev][0].get(), Bfield_fp[lev][1].get(), - Bfield_fp[lev][2].get() }); + Bfield_fp[lev][2].get() }, + do_pml_in_domain); pml[lev]->FillBoundaryB(patch_type); } const auto& period = Geom(lev).periodicity(); @@ -306,9 +312,10 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeB(patch_type, - { Bfield_cp[lev][0].get(), - Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get() }); + { Bfield_cp[lev][0].get(), + Bfield_cp[lev][1].get(), + Bfield_cp[lev][2].get() }, + do_pml_in_domain); pml[lev]->FillBoundaryB(patch_type); } const auto& cperiod = Geom(lev-1).periodicity(); @@ -331,7 +338,8 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type) { if (do_pml && pml[lev]->ok()) { - pml[lev]->ExchangeF(patch_type, F_fp[lev].get()); + pml[lev]->ExchangeF(patch_type, F_fp[lev].get(), + do_pml_in_domain); pml[lev]->FillBoundaryF(patch_type); } @@ -342,7 +350,8 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type) { if (do_pml && pml[lev]->ok()) { - pml[lev]->ExchangeF(patch_type, F_cp[lev].get()); + pml[lev]->ExchangeF(patch_type, F_cp[lev].get(), + do_pml_in_domain); pml[lev]->FillBoundaryF(patch_type); } diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index aaed96423..1aff7edfb 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -123,7 +123,6 @@ public: /// /// This deposits the particle charge onto a node-centered MultiFab and returns a unique ptr /// to it. The charge density is accumulated over all the particles in the MultiParticleContainer - /// This version uses PICSAR routines to perform the deposition. /// std::unique_ptr<amrex::MultiFab> GetChargeDensity(int lev, bool local = false); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 73724c7da..2dc25e6fa 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -946,7 +946,7 @@ PhysicalParticleContainer::Evolve (int lev, { BL_PROFILE("PPC::Evolve()"); BL_PROFILE_VAR_NS("PPC::Evolve::Copy", blp_copy); - BL_PROFILE_VAR_NS("PICSAR::FieldGather", blp_pxr_fg); + BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); @@ -1212,7 +1212,7 @@ PhysicalParticleContainer::Evolve (int lev, // // Field Gather of Aux Data (i.e., the full solution) // - BL_PROFILE_VAR_START(blp_pxr_fg); + BL_PROFILE_VAR_START(blp_fg); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrow(), e_is_nodal, @@ -1293,7 +1293,7 @@ PhysicalParticleContainer::Evolve (int lev, thread_num, lev, lev-1); } - BL_PROFILE_VAR_STOP(blp_pxr_fg); + BL_PROFILE_VAR_STOP(blp_fg); // // Particle Push @@ -1306,36 +1306,25 @@ PhysicalParticleContainer::Evolve (int lev, // // Current Deposition // - if (WarpX::use_picsar_deposition) { - // Deposit inside domains - DepositCurrentFortran(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, - 0, np_current, thread_num, - lev, lev, dt); - if (has_buffer){ - // Deposit in buffers - DepositCurrentFortran(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, - np_current, np-np_current, thread_num, - lev, lev-1, dt); - } + + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); } else { - int* AMREX_RESTRICT ion_lev; - if (do_field_ionization){ - ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); - } else { - ion_lev = nullptr; - } - - // Deposit inside domains - DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, - 0, np_current, thread_num, - lev, lev, dt); - if (has_buffer){ - // Deposit in buffers - DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz, - np_current, np-np_current, thread_num, - lev, lev-1, dt); - } + ion_lev = nullptr; } + + // Deposit inside domains + DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, + 0, np_current, thread_num, + lev, lev, dt); + if (has_buffer){ + // Deposit in buffers + DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz, + np_current, np-np_current, thread_num, + lev, lev-1, dt); + } + // // copy particle data back diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index bc46a116b..c39eec9dc 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -208,21 +208,6 @@ public: int depos_lev, amrex::Real dt); - virtual void DepositCurrentFortran(WarpXParIter& pti, - RealVector& wp, - RealVector& uxp, - RealVector& uyp, - RealVector& uzp, - amrex::MultiFab* jx, - amrex::MultiFab* jy, - amrex::MultiFab* jz, - const long offset, - const long np_to_depose, - int thread_num, - int lev, - int depos_lev, - amrex::Real dt); - // If particles start outside of the domain, ContinuousInjection // makes sure that they are initialized when they enter the domain, and // NOT before. Virtual function, overriden by derived classes. diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 8d2499a35..4e80374d8 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -263,133 +263,6 @@ WarpXParticleContainer::AddNParticles (int lev, Redistribute(); } -/* \brief Current Deposition for thread thread_num using PICSAR - * \param pti : Particle iterator - * \param wp : Array of particle weights - * \param uxp uyp uzp : Array of particle - * \param jx jy jz : Full array of current density - * \param offset : Index of first particle for which current is deposited - * \param np_to_depose: Number of particles for which current is deposited. - Particles [offset,offset+np_tp_depose] deposit their - current - * \param thread_num : Thread number (if tiling) - * \param lev : Level of box that contains particles - * \param depos_lev : Level on which particles deposit (if buffers are used) - * \param dt : Time step for particle level - */ -void -WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti, - RealVector& wp, RealVector& uxp, - RealVector& uyp, RealVector& uzp, - MultiFab* jx, MultiFab* jy, MultiFab* jz, - const long offset, const long np_to_depose, - int thread_num, int lev, int depos_lev, - Real dt) -{ - AMREX_ALWAYS_ASSERT_WITH_MESSAGE((depos_lev==(lev-1)) || - (depos_lev==(lev )), - "Deposition buffers only work for lev-1"); - // If no particles, do not do anything - if (np_to_depose == 0) return; - - const long ngJ = jx->nGrow(); - const std::array<Real,3>& dx = WarpX::CellSize(std::max(depos_lev,0)); - const long lvect = 8; - int j_is_nodal = jx->is_nodal() and jy->is_nodal() and jz->is_nodal(); - - BL_PROFILE_VAR_NS("PICSAR::CurrentDeposition", blp_pxr_cd); - BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); - - // Get tile box where current is deposited. - // The tile box is different when depositing in the buffers (depos_lev<lev) - // or when depositing inside the level (depos_lev=lev) - Box tilebox; - if (lev == depos_lev) { - tilebox = pti.tilebox(); - } else { - const IntVect& ref_ratio = WarpX::RefRatio(depos_lev); - tilebox = amrex::coarsen(pti.tilebox(),ref_ratio); - } - - // Staggered tile boxes (different in each direction) - Box tbx = convert(tilebox, WarpX::jx_nodal_flag); - Box tby = convert(tilebox, WarpX::jy_nodal_flag); - Box tbz = convert(tilebox, WarpX::jz_nodal_flag); - - // Lower corner of tile box physical domain - const std::array<Real,3>& xyzmin_tile = WarpX::LowerCorner(tilebox, depos_lev); - const std::array<Real, 3>& xyzmin = xyzmin_tile; - -#ifdef AMREX_USE_GPU - // No tiling on GPU: jx_ptr points to the full - // jx array (same for jy_ptr and jz_ptr). - Real* jx_ptr = (*jx)[pti].dataPtr(); - Real* jy_ptr = (*jy)[pti].dataPtr(); - Real* jz_ptr = (*jz)[pti].dataPtr(); - - auto jxntot = (*jx)[pti].length(); - auto jyntot = (*jy)[pti].length(); - auto jzntot = (*jz)[pti].length(); -#else - // Tiling is on: jx_ptr points to local_jx[thread_num] - // (same for jy_ptr and jz_ptr) - tbx.grow(ngJ); - tby.grow(ngJ); - tbz.grow(ngJ); - - local_jx[thread_num].resize(tbx, jx->nComp()); - local_jy[thread_num].resize(tby, jy->nComp()); - local_jz[thread_num].resize(tbz, jz->nComp()); - - Real* jx_ptr = local_jx[thread_num].dataPtr(); - Real* jy_ptr = local_jy[thread_num].dataPtr(); - Real* jz_ptr = local_jz[thread_num].dataPtr(); - - // local_jx[thread_num] is set to zero - local_jx[thread_num].setVal(0.0); - local_jy[thread_num].setVal(0.0); - local_jz[thread_num].setVal(0.0); - - auto jxntot = local_jx[thread_num].length(); - auto jyntot = local_jy[thread_num].length(); - auto jzntot = local_jz[thread_num].length(); -#endif - // GPU, no tiling: deposit directly in jx - // CPU, tiling: deposit into local_jx - // (same for jx and jz) - BL_PROFILE_VAR_START(blp_pxr_cd); - warpx_current_deposition( - jx_ptr, &ngJ, jxntot.getVect(), - jy_ptr, &ngJ, jyntot.getVect(), - jz_ptr, &ngJ, jzntot.getVect(), - &WarpX::n_rz_azimuthal_modes, - &np_to_depose, - m_xp[thread_num].dataPtr() + offset, - m_yp[thread_num].dataPtr() + offset, - m_zp[thread_num].dataPtr() + offset, - uxp.dataPtr() + offset, - uyp.dataPtr() + offset, - uzp.dataPtr() + offset, - m_giv[thread_num].dataPtr() + offset, - wp.dataPtr() + offset, &this->charge, - &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dt, &dx[0], &dx[1], &dx[2], - &WarpX::nox,&WarpX::noy,&WarpX::noz, &j_is_nodal, - &lvect,&WarpX::current_deposition_algo); - - BL_PROFILE_VAR_STOP(blp_pxr_cd); - -#ifndef AMREX_USE_GPU - BL_PROFILE_VAR_START(blp_accumulate); - // CPU, tiling: atomicAdd local_jx into jx - // (same for jx and jz) - (*jx)[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, local_jx[thread_num].nComp()); - (*jy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, local_jy[thread_num].nComp()); - (*jz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, local_jz[thread_num].nComp()); - BL_PROFILE_VAR_STOP(blp_accumulate); -#endif -} - /* \brief Current Deposition for thread thread_num * \param pti : Particle iterator * \param wp : Array of particle weights @@ -430,6 +303,8 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, const Real stagger_shift = j_is_nodal ? 0.0 : 0.5; BL_PROFILE_VAR_NS("PPC::Evolve::Accumulate", blp_accumulate); + BL_PROFILE_VAR_NS("PPC::CurrentDeposition", blp_deposit); + // Get tile box where current is deposited. // The tile box is different when depositing in the buffers (depos_lev<lev) @@ -492,6 +367,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // Better for memory? worth trying? const Dim3 lo = lbound(tilebox); + BL_PROFILE_VAR_START(blp_deposit); if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov) { if (WarpX::nox == 1){ doEsirkepovDepositionShapeN<1>( @@ -533,6 +409,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, stagger_shift, q); } } + BL_PROFILE_VAR_STOP(blp_deposit); #ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index 6a32513b7..54c721abf 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -5,8 +5,6 @@ #include <string> struct MaxwellSolverAlgo { - // These numbers corresponds to the algorithm code in WarpX's - // `warpx_push_bvec` and `warpx_push_evec_f` enum { Yee = 0, CKC = 1 @@ -14,8 +12,6 @@ struct MaxwellSolverAlgo { }; struct ParticlePusherAlgo { - // These numbers correspond to the algorithm code in WarpX's - // `warpx_particle_pusher` enum { Boris = 0, Vay = 1 @@ -24,12 +20,8 @@ struct ParticlePusherAlgo { struct CurrentDepositionAlgo { enum { - // These numbers corresponds to the algorithm code in PICSAR's - // `depose_jxjyjz_generic` and `depose_jxjyjz_generic_2d` - Direct = 3, - DirectVectorized = 2, - EsirkepovNonOptimized = 1, - Esirkepov = 0 + Esirkepov = 0, + Direct = 1 }; }; @@ -41,11 +33,9 @@ struct ChargeDepositionAlgo { }; struct GatheringAlgo { - // These numbers corresponds to the algorithm code in PICSAR's - // `geteb3d_energy_conserving_generic` function + // Only the Standard algorithm is implemented enum { - Vectorized = 0, - Standard = 1 + Standard = 0 }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index 842085a36..216199103 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -4,7 +4,7 @@ #include <cstring> // Define dictionary with correspondance between user-input strings, -// and corresponding integer for use inside the code (e.g. in PICSAR). +// and corresponding integer for use inside the code const std::map<std::string, int> maxwell_solver_algo_to_int = { {"yee", MaxwellSolverAlgo::Yee }, @@ -23,9 +23,6 @@ const std::map<std::string, int> particle_pusher_algo_to_int = { const std::map<std::string, int> current_deposition_algo_to_int = { {"esirkepov", CurrentDepositionAlgo::Esirkepov }, {"direct", CurrentDepositionAlgo::Direct }, -#if (!defined AMREX_USE_GPU)&&(AMREX_SPACEDIM == 3) // Only available on CPU and 3D - {"direct-vectorized", CurrentDepositionAlgo::DirectVectorized }, -#endif {"default", CurrentDepositionAlgo::Esirkepov } }; @@ -36,12 +33,7 @@ const std::map<std::string, int> charge_deposition_algo_to_int = { const std::map<std::string, int> gathering_algo_to_int = { {"standard", GatheringAlgo::Standard }, -#ifndef AMREX_USE_GPU // Only available on CPU - {"vectorized", GatheringAlgo::Vectorized }, - {"default", GatheringAlgo::Vectorized } -#else {"default", GatheringAlgo::Standard } -#endif }; diff --git a/Source/WarpX.H b/Source/WarpX.H index 199630171..9236b8e20 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -83,7 +83,6 @@ public: static amrex::Vector<amrex::Real> B_external; // Algorithms - static long use_picsar_deposition; static long current_deposition_algo; static long charge_deposition_algo; static long field_gathering_algo; @@ -196,6 +195,12 @@ public: void DampPML (int lev); void DampPML (int lev, PatchType patch_type); + void DampJPML (); + void DampJPML (int lev); + void DampJPML (int lev, PatchType patch_type); + + void CopyJPML (); + void PushParticlesandDepose (int lev, amrex::Real cur_time); void PushParticlesandDepose ( amrex::Real cur_time); @@ -504,6 +509,9 @@ private: int do_pml = 1; int pml_ncell = 10; int pml_delta = 10; + int pml_has_particles = 0; + int do_pml_j_damping = 0; + int do_pml_in_domain = 0; amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(); amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector(); amrex::Vector<std::unique_ptr<PML> > pml; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 8e8c25f4d..95826c075 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -38,7 +38,6 @@ Vector<int> WarpX::boost_direction = {0,0,0}; int WarpX::do_compute_max_step_from_zmax = 0; Real WarpX::zmax_plasma_to_compute_max_step = 0.; -long WarpX::use_picsar_deposition = 0; long WarpX::current_deposition_algo; long WarpX::charge_deposition_algo; long WarpX::field_gathering_algo; @@ -398,6 +397,9 @@ WarpX::ReadParameters () pp.query("do_pml", do_pml); pp.query("pml_ncell", pml_ncell); pp.query("pml_delta", pml_delta); + pp.query("pml_has_particles", pml_has_particles); + pp.query("do_pml_j_damping", do_pml_j_damping); + pp.query("do_pml_in_domain", do_pml_in_domain); Vector<int> parse_do_pml_Lo(AMREX_SPACEDIM,1); pp.queryarr("do_pml_Lo", parse_do_pml_Lo); @@ -414,6 +416,9 @@ WarpX::ReadParameters () do_pml_Hi[2] = parse_do_pml_Hi[2]; #endif + if ( (do_pml_j_damping==1)&&(do_pml_in_domain==0) ){ + amrex::Abort("J-damping can only be done when PML are inside simulation domain (do_pml_in_domain=1)"); + } pp.query("dump_openpmd", dump_openpmd); pp.query("openpmd_backend", openpmd_backend); @@ -533,7 +538,6 @@ WarpX::ReadParameters () { ParmParse pp("algo"); - pp.query("use_picsar_deposition", use_picsar_deposition); current_deposition_algo = GetAlgorithmInteger(pp, "current_deposition"); charge_deposition_algo = GetAlgorithmInteger(pp, "charge_deposition"); field_gathering_algo = GetAlgorithmInteger(pp, "field_gathering"); |