diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 8f3ca1da1d7..69fa14658d2 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -46,15 +46,18 @@ struct SigmaBox { SigmaBox (const amrex::Box& box, const amrex::BoxArray& grids, const amrex::Real* dx, const amrex::IntVect& ncell, const amrex::IntVect& delta, - const amrex::Box& regdomain, amrex::Real v_sigma); + const amrex::Box& regdomain, amrex::Real v_sigma, + const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength); void define_single (const amrex::Box& regdomain, const amrex::IntVect& ncell, const amrex::Array& fac, - amrex::Real v_sigma); + amrex::Real v_sigma, + const int do_cubic_sigma_pml); void define_multiple (const amrex::Box& box, const amrex::BoxArray& grids, const amrex::IntVect& ncell, const amrex::Array& fac, - amrex::Real v_sigma); + amrex::Real v_sigma, + const int do_cubic_sigma_pml); void ComputePMLFactorsB (const amrex::Real* dx, amrex::Real dt); void ComputePMLFactorsE (const amrex::Real* dx, amrex::Real dt); @@ -81,8 +84,9 @@ class SigmaBoxFactory public: SigmaBoxFactory (const amrex::BoxArray& grid_ba, const amrex::Real* dx, const amrex::IntVect& ncell, const amrex::IntVect& delta, - const amrex::Box& regular_domain, const amrex::Real v_sigma_sb) - : m_grids(grid_ba), m_dx(dx), m_ncell(ncell), m_delta(delta), m_regdomain(regular_domain), m_v_sigma_sb(v_sigma_sb) {} + const amrex::Box& regular_domain, const amrex::Real v_sigma_sb, + const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength) + : m_grids(grid_ba), m_dx(dx), m_ncell(ncell), m_delta(delta), m_regdomain(regular_domain), m_v_sigma_sb(v_sigma_sb), m_do_cubic_sigma_pml(do_cubic_sigma_pml), m_pml_damping_strength(pml_damping_strength) {} ~SigmaBoxFactory () override = default; SigmaBoxFactory (const SigmaBoxFactory&) = default; @@ -95,7 +99,7 @@ public: [[nodiscard]] SigmaBox* create (const amrex::Box& box, int /*ncomps*/, const amrex::FabInfo& /*info*/, int /*box_index*/) const final { - return new SigmaBox(box, m_grids, m_dx, m_ncell, m_delta, m_regdomain, m_v_sigma_sb); + return new SigmaBox(box, m_grids, m_dx, m_ncell, m_delta, m_regdomain, m_v_sigma_sb, m_do_cubic_sigma_pml, m_pml_damping_strength); } void destroy (SigmaBox* fab) const final @@ -116,6 +120,8 @@ private: amrex::IntVect m_delta; amrex::Box m_regdomain; amrex::Real m_v_sigma_sb; + int m_do_cubic_sigma_pml; + amrex::Real m_pml_damping_strength; }; class MultiSigmaBox @@ -125,7 +131,8 @@ public: MultiSigmaBox(const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, const amrex::BoxArray& grid_ba, const amrex::Real* dx, const amrex::IntVect& ncell, const amrex::IntVect& delta, - const amrex::Box& regular_domain, amrex::Real v_sigma_sb); + const amrex::Box& regular_domain, amrex::Real v_sigma_sb, + const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength); void ComputePMLFactorsB (const amrex::Real* dx, amrex::Real dt); void ComputePMLFactorsE (const amrex::Real* dx, amrex::Real dt); private: @@ -148,6 +155,7 @@ public: const amrex::IntVect& fill_guards_fields, const amrex::IntVect& fill_guards_current, int max_guard_EB, amrex::Real v_sigma_sb, + const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength, amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 340005e9211..6471a0ce050 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -64,11 +64,11 @@ namespace void FillLo (Sigma& sigma, Sigma& sigma_cumsum, Sigma& sigma_star, Sigma& sigma_star_cumsum, const int olo, const int ohi, const int glo, Real fac, - const amrex::Real v_sigma) + const amrex::Real v_sigma, const int do_cubic_sigma_pml) { const int slo = sigma.m_lo; const int sslo = sigma_star.m_lo; - + int sigma_exp = (do_cubic_sigma_pml==1) ? 3 : 2; const int N = ohi+1-olo+1; Real* p_sigma = sigma.data(); Real* p_sigma_cumsum = sigma_cumsum.data(); @@ -78,14 +78,22 @@ namespace { i += olo; Real offset = static_cast(glo-i); - p_sigma[i-slo] = fac*(offset*offset); + amrex::Real offset_exp = 1._rt; + for (int iexp = 0; iexp < sigma_exp; ++iexp){ + offset_exp *= offset; + } + p_sigma[i-slo] = fac*(offset_exp); // sigma_cumsum is the analytical integral of sigma function at same points than sigma - p_sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3._rt)/v_sigma; + p_sigma_cumsum[i-slo] = (fac*(offset*offset_exp)/(sigma_exp+1._rt))/v_sigma; if (i <= ohi+1) { offset = static_cast(glo-i) - 0.5_rt; - p_sigma_star[i-sslo] = fac*(offset*offset); + offset_exp = 1._rt; + for (int iexp = 0; iexp < sigma_exp; ++iexp){ + offset_exp *= offset; + } + p_sigma_star[i-sslo] = fac*(offset_exp); // sigma_star_cumsum is the analytical integral of sigma function at same points than sigma_star - p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3._rt)/v_sigma; + p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset_exp)/(sigma_exp + 1._rt))/v_sigma; } }); } @@ -93,10 +101,11 @@ namespace void FillHi (Sigma& sigma, Sigma& sigma_cumsum, Sigma& sigma_star, Sigma& sigma_star_cumsum, const int olo, const int ohi, const int ghi, Real fac, - const amrex::Real v_sigma) + const amrex::Real v_sigma, const int do_cubic_sigma_pml) { const int slo = sigma.m_lo; const int sslo = sigma_star.m_lo; + int sigma_exp = (do_cubic_sigma_pml==1) ? 3 : 2; const int N = ohi+1-olo+1; Real* p_sigma = sigma.data(); @@ -107,12 +116,20 @@ namespace { i += olo; Real offset = static_cast(i-ghi-1); - p_sigma[i-slo] = fac*(offset*offset); - p_sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3._rt)/v_sigma; + amrex::Real offset_exp = 1._rt; + for (int iexp = 0; iexp < sigma_exp; ++iexp){ + offset_exp *= offset; + } + p_sigma[i-slo] = fac*(offset_exp); + p_sigma_cumsum[i-slo] = (fac*(offset*offset_exp)/(sigma_exp+1._rt))/v_sigma; if (i <= ohi+1) { offset = static_cast(i-ghi) - 0.5_rt; - p_sigma_star[i-sslo] = fac*(offset*offset); - p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3._rt)/v_sigma; + offset_exp = 1._rt; + for (int iexp = 0; iexp < sigma_exp; ++iexp){ + offset_exp *= offset; + } + p_sigma_star[i-sslo] = fac*(offset_exp); + p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset_exp)/(sigma_exp + 1._rt))/v_sigma; } }); } @@ -146,7 +163,8 @@ namespace SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, const IntVect& ncell, - const IntVect& delta, const amrex::Box& regdomain, const amrex::Real v_sigma_sb) + const IntVect& delta, const amrex::Box& regdomain, const amrex::Real v_sigma_sb, + const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength) { BL_ASSERT(box.cellCentered()); @@ -182,22 +200,27 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, const sigma_star_cumsum_fac[idim].m_lo = lo[idim]; sigma_star_cumsum_fac[idim].m_hi = hi[idim]+1; } - + int sigma_exp = (do_cubic_sigma_pml==1) ? 3 : 2; Array fac; for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - fac[idim] = 4.0_rt*PhysConst::c/(dx[idim]*static_cast(delta[idim]*delta[idim])); + amrex::Real delta_exp = 1._rt; + for (int iexp = 0; iexp < sigma_exp; ++iexp) { + delta_exp *= static_cast(delta[idim]); + } + fac[idim] = pml_damping_strength*PhysConst::c/(dx[idim]*delta_exp); } if (regdomain.ok()) { // The union of the regular grids is a single box - define_single(regdomain, ncell, fac, v_sigma_sb); + define_single(regdomain, ncell, fac, v_sigma_sb, do_cubic_sigma_pml); } else { - define_multiple(box, grids, ncell, fac, v_sigma_sb); + define_multiple(box, grids, ncell, fac, v_sigma_sb, do_cubic_sigma_pml); } } void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell, const Array& fac, - const amrex::Real v_sigma_sb) + const amrex::Real v_sigma_sb, + const int do_cubic_sigma_pml) { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { const int slo = sigma[idim].lo(); @@ -211,7 +234,8 @@ void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell, if (ohi >= olo) { FillLo(sigma[idim], sigma_cumsum[idim], sigma_star[idim], sigma_star_cumsum[idim], - olo, ohi, dlo, fac[idim], v_sigma_sb); + olo, ohi, dlo, fac[idim], v_sigma_sb, + do_cubic_sigma_pml); } #if (AMREX_SPACEDIM != 1) @@ -231,7 +255,7 @@ void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell, if (ohi >= olo) { FillHi(sigma[idim], sigma_cumsum[idim], sigma_star[idim], sigma_star_cumsum[idim], - olo, ohi, dhi, fac[idim], v_sigma_sb); + olo, ohi, dhi, fac[idim], v_sigma_sb, do_cubic_sigma_pml); } } @@ -239,7 +263,8 @@ void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell, } void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const IntVect& ncell, - const Array& fac, const amrex::Real v_sigma_sb) + const Array& fac, const amrex::Real v_sigma_sb, + const int do_cubic_sigma_pml) { const std::vector >& isects = grids.intersections(box, false, ncell); @@ -309,7 +334,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int FillLo(sigma[idim], sigma_cumsum[idim], sigma_star[idim], sigma_star_cumsum[idim], looverlap.smallEnd(idim), looverlap.bigEnd(idim), - grid_box.smallEnd(idim), fac[idim], v_sigma_sb); + grid_box.smallEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml); } Box hibox = amrex::adjCellHi(grid_box, idim, ncell[idim]); @@ -322,7 +347,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int FillHi(sigma[idim], sigma_cumsum[idim], sigma_star[idim], sigma_star_cumsum[idim], hioverlap.smallEnd(idim), hioverlap.bigEnd(idim), - grid_box.bigEnd(idim), fac[idim], v_sigma_sb); + grid_box.bigEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml); } WARPX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -358,7 +383,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int FillLo(sigma[idim], sigma_cumsum[idim], sigma_star[idim], sigma_star_cumsum[idim], looverlap.smallEnd(idim), looverlap.bigEnd(idim), - grid_box.smallEnd(idim), fac[idim], v_sigma_sb); + grid_box.smallEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml); } Box hibox = amrex::adjCellHi(grid_box, idim, ncell[idim]); @@ -367,7 +392,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int FillHi(sigma[idim], sigma_cumsum[idim], sigma_star[idim], sigma_star_cumsum[idim], hioverlap.smallEnd(idim), hioverlap.bigEnd(idim), - grid_box.bigEnd(idim), fac[idim], v_sigma_sb); + grid_box.bigEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml); } WARPX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -408,7 +433,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int FillLo(sigma[idim], sigma_cumsum[idim], sigma_star[idim], sigma_star_cumsum[idim], looverlap.smallEnd(idim), looverlap.bigEnd(idim), - grid_box.smallEnd(idim), fac[idim], v_sigma_sb); + grid_box.smallEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml); } const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell[idim]); @@ -417,7 +442,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int FillHi(sigma[idim], sigma_cumsum[idim], sigma_star[idim], sigma_star_cumsum[idim], hioverlap.smallEnd(idim), hioverlap.bigEnd(idim), - grid_box.bigEnd(idim), fac[idim], v_sigma_sb); + grid_box.bigEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml); } WARPX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -508,9 +533,9 @@ SigmaBox::ComputePMLFactorsE (const Real* a_dx, Real dt) MultiSigmaBox::MultiSigmaBox (const BoxArray& ba, const DistributionMapping& dm, const BoxArray& grid_ba, const Real* dx, const IntVect& ncell, const IntVect& delta, - const amrex::Box& regular_domain, const amrex::Real v_sigma_sb) + const amrex::Box& regular_domain, const amrex::Real v_sigma_sb, const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength) : FabArray(ba,dm,1,0,MFInfo(), - SigmaBoxFactory(grid_ba,dx,ncell,delta, regular_domain, v_sigma_sb)) + SigmaBoxFactory(grid_ba,dx,ncell,delta, regular_domain, v_sigma_sb,do_cubic_sigma_pml,pml_damping_strength)) {} void @@ -557,6 +582,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const amrex::IntVect& fill_guards_fields, const amrex::IntVect& fill_guards_current, int max_guard_EB, const amrex::Real v_sigma_sb, + const int do_cubic_sigma_pml, amrex::Real pml_damping_strength, const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) : m_dive_cleaning(do_pml_dive_cleaning), m_divb_cleaning(do_pml_divb_cleaning), @@ -744,7 +770,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, Box single_domain_box = is_single_box_domain ? domain0 : Box(); // Empty box (i.e., Box()) means it's not a single box domain. sigba_fp = std::make_unique(ba, dm, grid_ba_reduced, geom->CellSize(), - IntVect(ncell), IntVect(delta), single_domain_box, v_sigma_sb); + IntVect(ncell), IntVect(delta), single_domain_box, v_sigma_sb, do_cubic_sigma_pml, pml_damping_strength); if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { #ifndef WARPX_USE_FFT @@ -859,7 +885,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, single_domain_box = is_single_box_domain ? cdomain : Box(); sigba_cp = std::make_unique(cba, cdm, grid_cba_reduced, cgeom->CellSize(), - cncells, cdelta, single_domain_box, v_sigma_sb); + cncells, cdelta, single_domain_box, v_sigma_sb, do_cubic_sigma_pml, pml_damping_strength); if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { #ifndef WARPX_USE_FFT diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp index ec696689373..90b2274d9c5 100644 --- a/Source/BoundaryConditions/WarpXEvolvePML.cpp +++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp @@ -77,7 +77,7 @@ WarpX::DampPML_Cartesian (const int lev, PatchType patch_type) { const bool dive_cleaning = WarpX::do_pml_dive_cleaning; const bool divb_cleaning = WarpX::do_pml_divb_cleaning; - + const bool abc_in_pml = WarpX::do_abc_in_pml; if (pml[lev]->ok()) { const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); @@ -154,19 +154,19 @@ WarpX::DampPML_Cartesian (const int lev, PatchType patch_type) warpx_damp_pml_ex(i, j, k, pml_Exfab, Ex_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo, - dive_cleaning); + dive_cleaning, abc_in_pml); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { warpx_damp_pml_ey(i, j, k, pml_Eyfab, Ey_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo, - dive_cleaning); + dive_cleaning, abc_in_pml); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { warpx_damp_pml_ez(i, j, k, pml_Ezfab, Ez_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo, - dive_cleaning); + dive_cleaning, abc_in_pml); }); amrex::ParallelFor(tbx, tby, tbz, @@ -174,19 +174,19 @@ WarpX::DampPML_Cartesian (const int lev, PatchType patch_type) warpx_damp_pml_bx(i, j, k, pml_Bxfab, Bx_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo, - divb_cleaning); + divb_cleaning, abc_in_pml); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { warpx_damp_pml_by(i, j, k, pml_Byfab, By_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo, - divb_cleaning); + divb_cleaning, abc_in_pml); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { warpx_damp_pml_bz(i, j, k, pml_Bzfab, Bz_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo, - divb_cleaning); + divb_cleaning, abc_in_pml); }); // For warpx_damp_pml_F(), mfi.nodaltilebox is used in the ParallelFor loop and here we diff --git a/Source/BoundaryConditions/WarpX_PML_kernels.H b/Source/BoundaryConditions/WarpX_PML_kernels.H index 7aaf5a57369..97003d8d2b6 100644 --- a/Source/BoundaryConditions/WarpX_PML_kernels.H +++ b/Source/BoundaryConditions/WarpX_PML_kernels.H @@ -24,11 +24,11 @@ void warpx_damp_pml_ex (int i, int j, int k, amrex::Array4 const& E const amrex::Real* const sigma_star_fac_y, const amrex::Real* const sigma_star_fac_z, int xlo, int ylo, int zlo, - const bool dive_cleaning) + const bool dive_cleaning, const bool abc_in_pml) { #if defined(WARPX_DIM_1D_Z) amrex::ignore_unused(i, j, k, Ex, Ex_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, - sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning); + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning, abc_in_pml); amrex::Abort("PML not implemented in Cartesian 1D geometry"); #endif @@ -53,8 +53,14 @@ void warpx_damp_pml_ex (int i, int j, int k, amrex::Array4 const& E // Exz if (sz == 0) { Ex(i,j,k,PMLComp::xz) *= sigma_star_fac_z[j-zlo]; + if (abc_in_pml) { + Ex(i,j,k,PMLComp::xy) *= sigma_star_fac_z[j-zlo]; + } } else { Ex(i,j,k,PMLComp::xz) *= sigma_fac_z[j-zlo]; + if (abc_in_pml) { + Ex(i,j,k,PMLComp::xy) *= sigma_fac_z[j-zlo]; + } } #elif defined(WARPX_DIM_3D) @@ -77,15 +83,27 @@ void warpx_damp_pml_ex (int i, int j, int k, amrex::Array4 const& E // Exy if (sy == 0) { Ex(i,j,k,PMLComp::xy) *= sigma_star_fac_y[j-ylo]; + if (abc_in_pml) { + Ex(i,j,k,PMLComp::xz) *= sigma_star_fac_y[j-ylo]; + } } else { Ex(i,j,k,PMLComp::xy) *= sigma_fac_y[j-ylo]; + if (abc_in_pml) { + Ex(i,j,k,PMLComp::xz) *= sigma_fac_y[j-ylo]; + } } // Exz if (sz == 0) { Ex(i,j,k,PMLComp::xz) *= sigma_star_fac_z[k-zlo]; + if (abc_in_pml) { + Ex(i,j,k,PMLComp::xy) *= sigma_star_fac_z[k-zlo]; + } } else { Ex(i,j,k,PMLComp::xz) *= sigma_fac_z[k-zlo]; + if (abc_in_pml) { + Ex(i,j,k,PMLComp::xy) *= sigma_fac_z[k-zlo]; + } } #endif @@ -101,17 +119,17 @@ void warpx_damp_pml_ey (int i, int j, int k, amrex::Array4 const& E const amrex::Real* const sigma_star_fac_y, const amrex::Real* const sigma_star_fac_z, int xlo, int ylo, int zlo, - const bool dive_cleaning) + const bool dive_cleaning, const bool abc_in_pml) { #if defined(WARPX_DIM_1D_Z) amrex::ignore_unused(i, j, k, Ey, Ey_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, - sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning); + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning, abc_in_pml); amrex::Abort("PML not implemented in Cartesian 1D geometry"); #endif #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo, dive_cleaning); + amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo, dive_cleaning, abc_in_pml); // sx = 0 means that Ey is staggered in x, while sx = 1 means that Ey is nodal in x (same for z) const int sx = Ey_stag[0]; @@ -120,15 +138,27 @@ void warpx_damp_pml_ey (int i, int j, int k, amrex::Array4 const& E // Eyx if (sx == 0) { Ey(i,j,k,PMLComp::yx) *= sigma_star_fac_x[i-xlo]; + if (abc_in_pml) { + Ey(i,j,k,PMLComp::yz) *= sigma_star_fac_x[i-xlo]; + } } else { Ey(i,j,k,PMLComp::yx) *= sigma_fac_x[i-xlo]; + if (abc_in_pml) { + Ey(i,j,k,PMLComp::yz) *= sigma_fac_x[i-xlo]; + } } // Eyz if (sz == 0) { Ey(i,j,k,PMLComp::yz) *= sigma_star_fac_z[j-zlo]; + if (abc_in_pml) { + Ey(i,j,k,PMLComp::yx) *= sigma_star_fac_z[j-zlo]; + } } else { Ey(i,j,k,PMLComp::yz) *= sigma_fac_z[j-zlo]; + if (abc_in_pml) { + Ey(i,j,k,PMLComp::yx) *= sigma_fac_z[j-zlo]; + } } #elif defined(WARPX_DIM_3D) @@ -141,8 +171,14 @@ void warpx_damp_pml_ey (int i, int j, int k, amrex::Array4 const& E // Eyx if (sx == 0) { Ey(i,j,k,PMLComp::yx) *= sigma_star_fac_x[i-xlo]; + if (abc_in_pml) { + Ey(i,j,k,PMLComp::yz) *= sigma_star_fac_x[i-xlo]; + } } else { Ey(i,j,k,PMLComp::yx) *= sigma_fac_x[i-xlo]; + if (abc_in_pml) { + Ey(i,j,k,PMLComp::yz) *= sigma_fac_x[i-xlo]; + } } if (dive_cleaning) @@ -158,8 +194,14 @@ void warpx_damp_pml_ey (int i, int j, int k, amrex::Array4 const& E // Eyz if (sz == 0) { Ey(i,j,k,PMLComp::yz) *= sigma_star_fac_z[k-zlo]; + if (abc_in_pml) { + Ey(i,j,k,PMLComp::yx) *= sigma_star_fac_z[k-zlo]; + } } else { Ey(i,j,k,PMLComp::yz) *= sigma_fac_z[k-zlo]; + if (abc_in_pml) { + Ey(i,j,k,PMLComp::yx) *= sigma_fac_z[k-zlo]; + } } #endif @@ -175,11 +217,11 @@ void warpx_damp_pml_ez (int i, int j, int k, amrex::Array4 const& E const amrex::Real* const sigma_star_fac_y, const amrex::Real* const sigma_star_fac_z, int xlo, int ylo, int zlo, - const bool dive_cleaning) + const bool dive_cleaning, const bool abc_in_pml) { #if defined(WARPX_DIM_1D_Z) amrex::ignore_unused(i, j, k, Ez, Ez_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, - sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning); + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning, abc_in_pml); amrex::Abort("PML not implemented in Cartesian 1D geometry"); #endif @@ -194,8 +236,14 @@ void warpx_damp_pml_ez (int i, int j, int k, amrex::Array4 const& E // Ezx if (sx == 0) { Ez(i,j,k,PMLComp::zx) *= sigma_star_fac_x[i-xlo]; + if (abc_in_pml) { + Ez(i,j,k,PMLComp::zy) *= sigma_star_fac_x[i-xlo]; + } } else { Ez(i,j,k,PMLComp::zx) *= sigma_fac_x[i-xlo]; + if (abc_in_pml) { + Ez(i,j,k,PMLComp::zy) *= sigma_fac_x[i-xlo]; + } } if (dive_cleaning) @@ -218,15 +266,27 @@ void warpx_damp_pml_ez (int i, int j, int k, amrex::Array4 const& E // Ezx if (sx == 0) { Ez(i,j,k,PMLComp::zx) *= sigma_star_fac_x[i-xlo]; + if (abc_in_pml) { + Ez(i,j,k,PMLComp::zy) *= sigma_star_fac_x[i-xlo]; + } } else { Ez(i,j,k,PMLComp::zx) *= sigma_fac_x[i-xlo]; + if (abc_in_pml) { + Ez(i,j,k,PMLComp::zy) *= sigma_fac_x[i-xlo]; + } } // Ezy if (sy == 0) { Ez(i,j,k,PMLComp::zy) *= sigma_star_fac_y[j-ylo]; + if (abc_in_pml) { + Ez(i,j,k,PMLComp::zx) *= sigma_star_fac_y[j-ylo]; + } } else { Ez(i,j,k,PMLComp::zy) *= sigma_fac_y[j-ylo]; + if (abc_in_pml) { + Ez(i,j,k,PMLComp::zx) *= sigma_fac_y[j-ylo]; + } } if (dive_cleaning) @@ -252,11 +312,11 @@ void warpx_damp_pml_bx (int i, int j, int k, amrex::Array4 const& B const amrex::Real* const sigma_star_fac_y, const amrex::Real* const sigma_star_fac_z, int xlo, int ylo, int zlo, - const bool divb_cleaning) + const bool divb_cleaning, const bool abc_in_pml) { #if defined(WARPX_DIM_1D_Z) amrex::ignore_unused(i, j, k, Bx, Bx_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, - sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning); + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning, abc_in_pml); amrex::Abort("PML not implemented in Cartesian 1D geometry"); #endif @@ -281,8 +341,14 @@ void warpx_damp_pml_bx (int i, int j, int k, amrex::Array4 const& B // Bxz if (sz == 0) { Bx(i,j,k,PMLComp::xz) *= sigma_star_fac_z[j-zlo]; + if (abc_in_pml) { + Bx(i,j,k,PMLComp::xy) *= sigma_star_fac_z[j-zlo]; + } } else { Bx(i,j,k,PMLComp::xz) *= sigma_fac_z[j-zlo]; + if (abc_in_pml) { + Bx(i,j,k,PMLComp::xy) *= sigma_fac_z[j-zlo]; + } } #elif defined(WARPX_DIM_3D) @@ -305,15 +371,27 @@ void warpx_damp_pml_bx (int i, int j, int k, amrex::Array4 const& B // Bxy if (sy == 0) { Bx(i,j,k,PMLComp::xy) *= sigma_star_fac_y[j-ylo]; + if (abc_in_pml) { + Bx(i,j,k,PMLComp::xz) *= sigma_star_fac_y[j-ylo]; + } } else { Bx(i,j,k,PMLComp::xy) *= sigma_fac_y[j-ylo]; + if (abc_in_pml) { + Bx(i,j,k,PMLComp::xz) *= sigma_fac_y[j-ylo]; + } } // Bxz if (sz == 0) { Bx(i,j,k,PMLComp::xz) *= sigma_star_fac_z[k-zlo]; + if (abc_in_pml) { + Bx(i,j,k,PMLComp::xy) *= sigma_star_fac_z[k-zlo]; + } } else { Bx(i,j,k,PMLComp::xz) *= sigma_fac_z[k-zlo]; + if (abc_in_pml) { + Bx(i,j,k,PMLComp::xy) *= sigma_fac_z[k-zlo]; + } } #endif @@ -329,11 +407,11 @@ void warpx_damp_pml_by (int i, int j, int k, amrex::Array4 const& B const amrex::Real* const sigma_star_fac_y, const amrex::Real* const sigma_star_fac_z, int xlo, int ylo, int zlo, - const bool divb_cleaning) + const bool divb_cleaning, const bool abc_in_pml) { #if defined(WARPX_DIM_1D_Z) amrex::ignore_unused(i, j, k, By, By_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, - sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning); + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning, abc_in_pml); amrex::Abort("PML not implemented in Cartesian 1D geometry"); #endif @@ -348,15 +426,27 @@ void warpx_damp_pml_by (int i, int j, int k, amrex::Array4 const& B // Byx if (sx == 0) { By(i,j,k,PMLComp::yx) *= sigma_star_fac_x[i-xlo]; + if (abc_in_pml) { + By(i,j,k,PMLComp::yz) *= sigma_star_fac_x[i-xlo]; + } } else { By(i,j,k,PMLComp::yx) *= sigma_fac_x[i-xlo]; + if (abc_in_pml) { + By(i,j,k,PMLComp::yz) *= sigma_fac_x[i-xlo]; + } } // Byz if (sz == 0) { By(i,j,k,PMLComp::yz) *= sigma_star_fac_z[j-zlo]; + if (abc_in_pml) { + By(i,j,k,PMLComp::yx) *= sigma_star_fac_z[j-zlo]; + } } else { By(i,j,k,PMLComp::yz) *= sigma_fac_z[j-zlo]; + if (abc_in_pml) { + By(i,j,k,PMLComp::yx) *= sigma_fac_z[j-zlo]; + } } #elif defined(WARPX_DIM_3D) @@ -369,8 +459,14 @@ void warpx_damp_pml_by (int i, int j, int k, amrex::Array4 const& B // Byx if (sx == 0) { By(i,j,k,PMLComp::yx) *= sigma_star_fac_x[i-xlo]; + if (abc_in_pml) { + By(i,j,k,PMLComp::yz) *= sigma_star_fac_x[i-xlo]; + } } else { By(i,j,k,PMLComp::yx) *= sigma_fac_x[i-xlo]; + if (abc_in_pml) { + By(i,j,k,PMLComp::yz) *= sigma_fac_x[i-xlo]; + } } if (divb_cleaning) @@ -386,8 +482,14 @@ void warpx_damp_pml_by (int i, int j, int k, amrex::Array4 const& B // Byz if (sz == 0) { By(i,j,k,PMLComp::yz) *= sigma_star_fac_z[k-zlo]; + if (abc_in_pml) { + By(i,j,k,PMLComp::yx) *= sigma_star_fac_z[k-zlo]; + } } else { By(i,j,k,PMLComp::yz) *= sigma_fac_z[k-zlo]; + if (abc_in_pml) { + By(i,j,k,PMLComp::yx) *= sigma_fac_z[k-zlo]; + } } #endif @@ -403,11 +505,11 @@ void warpx_damp_pml_bz (int i, int j, int k, amrex::Array4 const& B const amrex::Real* const sigma_star_fac_y, const amrex::Real* const sigma_star_fac_z, int xlo, int ylo, int zlo, - const bool divb_cleaning) + const bool divb_cleaning, const bool abc_in_pml) { #if defined(WARPX_DIM_1D_Z) amrex::ignore_unused(i, j, k, Bz, Bz_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, - sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning); + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning, abc_in_pml); amrex::Abort("PML not implemented in Cartesian 1D geometry"); #endif @@ -422,8 +524,14 @@ void warpx_damp_pml_bz (int i, int j, int k, amrex::Array4 const& B // Bzx if (sx == 0) { Bz(i,j,k,PMLComp::zx) *= sigma_star_fac_x[i-xlo]; + if (abc_in_pml) { + Bz(i,j,k,PMLComp::zy) *= sigma_star_fac_x[i-xlo]; + } } else { Bz(i,j,k,PMLComp::zx) *= sigma_fac_x[i-xlo]; + if (abc_in_pml) { + Bz(i,j,k,PMLComp::zy) *= sigma_fac_x[i-xlo]; + } } if (divb_cleaning) @@ -446,15 +554,27 @@ void warpx_damp_pml_bz (int i, int j, int k, amrex::Array4 const& B // Bzx if (sx == 0) { Bz(i,j,k,PMLComp::zx) *= sigma_star_fac_x[i-xlo]; + if (abc_in_pml) { + Bz(i,j,k,PMLComp::zy) *= sigma_star_fac_x[i-xlo]; + } } else { Bz(i,j,k,PMLComp::zx) *= sigma_fac_x[i-xlo]; + if (abc_in_pml) { + Bz(i,j,k,PMLComp::zy) *= sigma_fac_x[i-xlo]; + } } // Bzy if (sy == 0) { Bz(i,j,k,PMLComp::zy) *= sigma_star_fac_y[j-ylo]; + if (abc_in_pml) { + Bz(i,j,k,PMLComp::zx) *= sigma_star_fac_y[j-ylo]; + } } else { Bz(i,j,k,PMLComp::zy) *= sigma_fac_y[j-ylo]; + if (abc_in_pml) { + Bz(i,j,k,PMLComp::zx) *= sigma_fac_y[j-ylo]; + } } if (divb_cleaning) diff --git a/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt b/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt index 2a5cc87c0cb..7fcdb832a36 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt +++ b/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt @@ -14,5 +14,6 @@ foreach(D IN LISTS WarpX_DIMS) BackTransformParticleFunctor.cpp ParticleReductionFunctor.cpp TemperatureFunctor.cpp + ProcessorMapFunctor.cpp ) endforeach() diff --git a/Source/Diagnostics/ComputeDiagFunctors/ProcessorMapFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/ProcessorMapFunctor.H new file mode 100644 index 00000000000..506f2501d43 --- /dev/null +++ b/Source/Diagnostics/ComputeDiagFunctors/ProcessorMapFunctor.H @@ -0,0 +1,45 @@ +#ifndef WARPX_PROCESSORMAPFUNCTOR_H_ +#define WARPX_PROCESSORMAPFUNCTOR_H_ + +#include "ComputeDiagFunctor.H" + +#include + +/** + * \brief Functor to cell-center MF and store result in mf_out. + */ +class ProcessorMapFunctor final : public ComputeDiagFunctor +{ +public: + /** Constructor. + * + * \param[in] mf_src source multifab. + * \param[in] lev level of multifab. Used for averaging in rz. + * \param[in] crse_ratio for interpolating field values from the simulation MultiFab, src_mf, + to the output diagnostic MultiFab, mf_dst. + * \param[in] convertRZmodes2cartesian (in cylindrical) whether to + * sum all modes in mf_src before cell-centering into dst multifab. + * \param[in] ncomp Number of component of mf_src to cell-center in dst multifab. + */ + ProcessorMapFunctor(const amrex::MultiFab * mf_src, int lev, + amrex::IntVect crse_ratio, + bool convertRZmodes2cartesian=true, int ncomp=1); + /** \brief Cell-center m_mf_src and write the result in mf_dst. + * + * In cylindrical geometry, by default this functor average all components + * of a MultiFab and writes into one single component. + * + * \param[out] mf_dst output MultiFab where the result is written + * \param[in] dcomp first component of mf_dst in which cell-centered + * data is stored + */ + void operator()(amrex::MultiFab& mf_dst, int dcomp, int /*i_buffer=0*/) const override; +private: + /** pointer to source multifab (can be multi-component) */ + amrex::MultiFab const * const m_mf_src = nullptr; + int m_lev; /**< level on which mf_src is defined (used in cylindrical) */ + /**< (for cylindrical) whether to average all modes into 1 comp */ + bool m_convertRZmodes2cartesian; +}; + +#endif // WARPX_CELLCENTERFUNCTOR_H_ diff --git a/Source/Diagnostics/ComputeDiagFunctors/ProcessorMapFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/ProcessorMapFunctor.cpp new file mode 100644 index 00000000000..33cb8cb9233 --- /dev/null +++ b/Source/Diagnostics/ComputeDiagFunctors/ProcessorMapFunctor.cpp @@ -0,0 +1,37 @@ +#include "ProcessorMapFunctor.H" + +#include "WarpX.H" + +#include +#include +#include + +ProcessorMapFunctor::ProcessorMapFunctor(amrex::MultiFab const * mf_src, int lev, + amrex::IntVect crse_ratio, + bool convertRZmodes2cartesian, int ncomp) + : ComputeDiagFunctor(ncomp, crse_ratio), m_mf_src(mf_src), m_lev(lev), + m_convertRZmodes2cartesian(convertRZmodes2cartesian) +{} + +void +ProcessorMapFunctor::operator()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer*/) const +{ + std::unique_ptr tmp; + tmp = std::make_unique(m_mf_src->boxArray(), m_mf_src->DistributionMap(), 1, m_mf_src->nGrowVect()); + + auto& warpx = WarpX::GetInstance(); + const amrex::DistributionMapping& dm = warpx.DistributionMap(m_lev); + for (amrex::MFIter mfi(*tmp, false); mfi.isValid(); ++mfi) + { + auto bx = mfi.growntilebox(tmp->nGrowVect()); + auto arr = tmp->array(mfi); + int iproc = dm[mfi.index()]; + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + arr(i,j,k) = iproc; + }); + } + InterpolateMFForDiag(mf_dst, *tmp, dcomp, warpx.DistributionMap(m_lev), + m_convertRZmodes2cartesian); +} diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp index 55af73d6408..8c3da11042a 100644 --- a/Source/Diagnostics/FullDiagnostics.cpp +++ b/Source/Diagnostics/FullDiagnostics.cpp @@ -10,6 +10,7 @@ #include "ComputeDiagFunctors/ParticleReductionFunctor.H" #include "ComputeDiagFunctors/TemperatureFunctor.H" #include "ComputeDiagFunctors/RhoFunctor.H" +#include "ComputeDiagFunctors/ProcessorMapFunctor.H" #include "Diagnostics/Diagnostics.H" #include "Diagnostics/ParticleDiag/ParticleDiag.H" #include "FieldSolver/Fields.H" @@ -697,6 +698,8 @@ FullDiagnostics::InitializeFieldFunctors (int lev) m_all_field_functors[lev][comp] = std::make_unique(warpx.getFieldPointerArray(FieldType::Bfield_aux, lev), lev, m_crse_ratio); } else if ( m_varnames[comp] == "divE" ){ m_all_field_functors[lev][comp] = std::make_unique(warpx.getFieldPointerArray(FieldType::Efield_aux, lev), lev, m_crse_ratio); + } else if ( m_varnames[comp] == "proc_number" ){ + m_all_field_functors[lev][comp] = std::make_unique(warpx.getFieldPointer(FieldType::Efield_aux, lev, 0), lev, m_crse_ratio); } else { diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index f88f64798f8..715c25d2e8b 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -640,6 +640,7 @@ WarpX::InitPML () amrex::IntVect(0), amrex::IntVect(0), guard_cells.ng_FieldSolver.max(), v_particle_pml, + do_cubic_sigma_pml, pml_damping_strength, do_pml_Lo[0], do_pml_Hi[0]); #endif @@ -679,6 +680,7 @@ WarpX::InitPML () amrex::IntVect(0), amrex::IntVect(0), guard_cells.ng_FieldSolver.max(), v_particle_pml, + do_cubic_sigma_pml, pml_damping_strength, do_pml_Lo[lev], do_pml_Hi[lev]); } } @@ -864,12 +866,12 @@ WarpX::InitLevelData (int lev, Real /*time*/) // Externally imposed fields are only initialized until the user-defined maxlevel_extEMfield_init. // The default maxlevel_extEMfield_init value is the total number of levels in the simulation if ((m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::parse_ext_grid_function) - && (lev > 0) && (lev <= maxlevel_extEMfield_init)) { + && (lev <= maxlevel_extEMfield_init)) { InitializeExternalFieldsOnGridUsingParser( - Bfield_aux[lev][0].get(), - Bfield_aux[lev][1].get(), - Bfield_aux[lev][2].get(), + Bfield_fp[lev][0].get(), + Bfield_fp[lev][1].get(), + Bfield_fp[lev][2].get(), m_p_ext_field_params->Bxfield_parser->compile<3>(), m_p_ext_field_params->Byfield_parser->compile<3>(), m_p_ext_field_params->Bzfield_parser->compile<3>(), @@ -878,17 +880,31 @@ WarpX::InitLevelData (int lev, Real /*time*/) 'B', lev, PatchType::fine); - InitializeExternalFieldsOnGridUsingParser( - Bfield_cp[lev][0].get(), - Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get(), - m_p_ext_field_params->Bxfield_parser->compile<3>(), - m_p_ext_field_params->Byfield_parser->compile<3>(), - m_p_ext_field_params->Bzfield_parser->compile<3>(), - m_edge_lengths[lev], - m_face_areas[lev], - 'B', - lev, PatchType::coarse); + if (lev > 0) { + InitializeExternalFieldsOnGridUsingParser( + Bfield_aux[lev][0].get(), + Bfield_aux[lev][1].get(), + Bfield_aux[lev][2].get(), + m_p_ext_field_params->Bxfield_parser->compile<3>(), + m_p_ext_field_params->Byfield_parser->compile<3>(), + m_p_ext_field_params->Bzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'B', + lev, PatchType::fine); + + InitializeExternalFieldsOnGridUsingParser( + Bfield_cp[lev][0].get(), + Bfield_cp[lev][1].get(), + Bfield_cp[lev][2].get(), + m_p_ext_field_params->Bxfield_parser->compile<3>(), + m_p_ext_field_params->Byfield_parser->compile<3>(), + m_p_ext_field_params->Bzfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'B', + lev, PatchType::coarse); + } } ParmParse pp_warpx("warpx"); int IncludeBfieldPerturbation = 0; @@ -937,6 +953,18 @@ WarpX::InitLevelData (int lev, Real /*time*/) if ((m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::parse_ext_grid_function) && (lev <= maxlevel_extEMfield_init)) { + InitializeExternalFieldsOnGridUsingParser( + Efield_fp[lev][0].get(), + Efield_fp[lev][1].get(), + Efield_fp[lev][2].get(), + m_p_ext_field_params->Exfield_parser->compile<3>(), + m_p_ext_field_params->Eyfield_parser->compile<3>(), + m_p_ext_field_params->Ezfield_parser->compile<3>(), + m_edge_lengths[lev], + m_face_areas[lev], + 'E', + lev, PatchType::fine); + #ifdef AMREX_USE_EB // We initialize ECTRhofield consistently with the Efield if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { @@ -1390,22 +1418,23 @@ WarpX::LoadExternalFields (int const lev) "External fields from file are not compatible with the moving window." ); } - // External grid fields - if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::parse_ext_grid_function) { - // Initialize Bfield_fp_external with external function - InitializeExternalFieldsOnGridUsingParser( - Bfield_fp_external[lev][0].get(), - Bfield_fp_external[lev][1].get(), - Bfield_fp_external[lev][2].get(), - m_p_ext_field_params->Bxfield_parser->compile<3>(), - m_p_ext_field_params->Byfield_parser->compile<3>(), - m_p_ext_field_params->Bzfield_parser->compile<3>(), - m_edge_lengths[lev], - m_face_areas[lev], - 'B', - lev, PatchType::fine); - } - else if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) { + //// External grid fields + //if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::parse_ext_grid_function) { + // // Initialize Bfield_fp_external with external function + // InitializeExternalFieldsOnGridUsingParser( + // Bfield_fp_external[lev][0].get(), + // Bfield_fp_external[lev][1].get(), + // Bfield_fp_external[lev][2].get(), + // m_p_ext_field_params->Bxfield_parser->compile<3>(), + // m_p_ext_field_params->Byfield_parser->compile<3>(), + // m_p_ext_field_params->Bzfield_parser->compile<3>(), + // m_edge_lengths[lev], + // m_face_areas[lev], + // 'B', + // lev, PatchType::fine); + //} + //else + if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) { #if defined(WARPX_DIM_RZ) WARPX_ALWAYS_ASSERT_WITH_MESSAGE(n_rz_azimuthal_modes == 1, "External field reading is not implemented for more than one RZ mode (see #3829)"); @@ -1419,21 +1448,22 @@ WarpX::LoadExternalFields (int const lev) #endif } - if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::parse_ext_grid_function) { - // Initialize Efield_fp_external with external function - InitializeExternalFieldsOnGridUsingParser( - Efield_fp_external[lev][0].get(), - Efield_fp_external[lev][1].get(), - Efield_fp_external[lev][2].get(), - m_p_ext_field_params->Exfield_parser->compile<3>(), - m_p_ext_field_params->Eyfield_parser->compile<3>(), - m_p_ext_field_params->Ezfield_parser->compile<3>(), - m_edge_lengths[lev], - m_face_areas[lev], - 'E', - lev, PatchType::fine); - } - else if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) { + //if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::parse_ext_grid_function) { + // // Initialize Efield_fp_external with external function + // InitializeExternalFieldsOnGridUsingParser( + // Efield_fp_external[lev][0].get(), + // Efield_fp_external[lev][1].get(), + // Efield_fp_external[lev][2].get(), + // m_p_ext_field_params->Exfield_parser->compile<3>(), + // m_p_ext_field_params->Eyfield_parser->compile<3>(), + // m_p_ext_field_params->Ezfield_parser->compile<3>(), + // m_edge_lengths[lev], + // m_face_areas[lev], + // 'E', + // lev, PatchType::fine); + //} + //else + if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) { #if defined(WARPX_DIM_RZ) WARPX_ALWAYS_ASSERT_WITH_MESSAGE(n_rz_azimuthal_modes == 1, "External field reading is not implemented for more than one RZ mode (see #3829)"); diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp index 58df62974bd..594875d55b0 100644 --- a/Source/Parallelization/WarpXRegrid.cpp +++ b/Source/Parallelization/WarpXRegrid.cpp @@ -20,6 +20,8 @@ #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXProfilerWrapper.H" +#include + #include #include #include @@ -87,68 +89,175 @@ WarpX::LoadBalance () int loadBalancedAnyLevel = false; const int nLevels = finestLevel(); - for (int lev = 0; lev <= nLevels; ++lev) - { + if (do_SFC_dm_vectorlevel) { + amrex::Vector rcost; + amrex::Vector current_pmap; + for (int lev = 0; lev <= nLevels; ++lev) + { + amrex::Vector rcost_lev(costs[lev]->size()); + amrex::ParallelDescriptor::GatherLayoutDataToVector(*costs[lev],rcost_lev, + amrex::ParallelDescriptor::IOProcessorNumber()); + rcost.insert(rcost.end(), rcost_lev.begin(), rcost_lev.end()); + auto& pmap_lev = costs[lev]->DistributionMap().ProcessorMap(); + current_pmap.insert(current_pmap.end(), pmap_lev.begin(), pmap_lev.end()); + } int doLoadBalance = false; - - // Compute the new distribution mapping - DistributionMapping newdm; - const amrex::Real nboxes = costs[lev]->size(); - const amrex::Real nprocs = ParallelContext::NProcsSub(); - const int nmax = static_cast(std::ceil(nboxes/nprocs*load_balance_knapsack_factor)); - // These store efficiency (meaning, the average 'cost' over all ranks, - // normalized to max cost) for current and proposed distribution mappings - amrex::Real currentEfficiency = 0.0; - amrex::Real proposedEfficiency = 0.0; - amrex::Print() << " load balance with sfc : " << load_balance_with_sfc << "\n"; - newdm = (load_balance_with_sfc) - ? DistributionMapping::makeSFC(*costs[lev], - currentEfficiency, proposedEfficiency, - false, - ParallelDescriptor::IOProcessorNumber()) - : DistributionMapping::makeKnapSack(*costs[lev], - currentEfficiency, proposedEfficiency, - nmax, - false, - ParallelDescriptor::IOProcessorNumber()); - // As specified in the above calls to makeSFC and makeKnapSack, the new - // distribution mapping is NOT communicated to all ranks; the loadbalanced - // dm is up-to-date only on root, and we can decide whether to broadcast - if ((load_balance_efficiency_ratio_threshold > 0.0) - && (ParallelDescriptor::MyProc() == ParallelDescriptor::IOProcessorNumber())) + amrex::Real currentEfficiency = 0.; + amrex::Real proposedEfficiency = 0.; + const amrex::Real nboxes = rcost.size(); +// const amrex::Real nprocs = ParallelContext::NProcsSub(); +// const int nmax = static_cast(std::ceil(nboxes/nprocs*load_balance_knapsack_factor)); + + amrex::BoxArray refined_ba = boxArray(0); + for (int lev = 1; lev <= nLevels; ++lev) { - doLoadBalance = (proposedEfficiency > load_balance_efficiency_ratio_threshold*currentEfficiency); + refined_ba.refine(refRatio(lev-1)); + amrex::BoxList refined_bl = refined_ba.boxList(); + refined_bl.join(boxArray(lev).boxList()); + refined_ba = amrex::BoxArray(refined_bl); } - amrex::Print() << "Do loadbalance is : " << doLoadBalance << " nmax : " << nmax << " current efficiency : " << currentEfficiency << " proposed efficiency : " << proposedEfficiency << "\n"; + amrex::Vector newdm(nLevels+1); + amrex::DistributionMapping r; + if (ParallelDescriptor::MyProc() == ParallelDescriptor::IOProcessorNumber()) + { + std::vector cost(rcost.size()); + + amrex::Real wmax = *std::max_element(rcost.begin(), rcost.end()); + amrex::Real scale = (wmax == 0) ? 1.e9_rt : 1.e9_rt/wmax; + + for (int i = 0; i < rcost.size(); ++i) { + cost[i] = amrex::Long(rcost[i]*scale) + 1L; + } + // `sort` needs to be false here since there's a parallel reduce function + // in the processor map function, but we are executing only on root + int nprocs = ParallelDescriptor::NProcs(); + r.SFCProcessorMap(refined_ba, cost, nprocs, proposedEfficiency, false); + + amrex::DistributionMapping::ComputeDistributionMappingEfficiency(r, + rcost, + ¤tEfficiency); + + if ((load_balance_efficiency_ratio_threshold > 0.0)) + { + doLoadBalance = (proposedEfficiency > load_balance_efficiency_ratio_threshold*currentEfficiency); + if ((getistep(0)+1) == 100) { + doLoadBalance = true; + } + } + amrex::Print() << " doloadbalance : " << doLoadBalance << "\n"; + amrex::Print() << proposedEfficiency << "\n"; + amrex::Print() << currentEfficiency << "\n"; + } ParallelDescriptor::Bcast(&doLoadBalance, 1, ParallelDescriptor::IOProcessorNumber()); if (doLoadBalance) { - Vector pmap; + amrex::Vector pmap(rcost.size()); if (ParallelDescriptor::MyProc() == ParallelDescriptor::IOProcessorNumber()) { - pmap = newdm.ProcessorMap(); + pmap = r.ProcessorMap(); } else { pmap.resize(static_cast(nboxes)); } - ParallelDescriptor::Bcast(pmap.data(), pmap.size(), ParallelDescriptor::IOProcessorNumber()); + // Broadcast vector from which to construct new distribution mapping + amrex::ParallelDescriptor::Bcast(&pmap[0], pmap.size(), ParallelDescriptor::IOProcessorNumber()); - if (ParallelDescriptor::MyProc() != ParallelDescriptor::IOProcessorNumber()) + int lev_start = 0; + for (int lev = 0; lev <= nLevels; ++lev) { - newdm = DistributionMapping(pmap); + amrex::Vector pmap_lev(pmap.begin() + lev_start, + pmap.begin() + lev_start + costs[lev]->size()); + newdm[lev] = amrex::DistributionMapping(pmap_lev); + lev_start += costs[lev]->size(); } - RemakeLevel(lev, t_new[lev], boxArray(lev), newdm); - - // Record the load balance efficiency - setLoadBalanceEfficiency(lev, proposedEfficiency); + for (int lev = 0; lev <= nLevels; ++lev) + { + RemakeLevel(lev, t_new[lev], boxArray(lev), newdm[lev]); + setLoadBalanceEfficiency(lev, proposedEfficiency); + } } - loadBalancedAnyLevel = loadBalancedAnyLevel || doLoadBalance; + } else { + //if (do_similar_dm_refpatch) { + // for (int lev = nLevels; lev > 0; --lev) { + // ablastr::coarsen::sample::Coarsen(*costs[lev-1], *costs[lev],0,0,1,0,WarpX::RefRatio(lev-1)); + // } + //} + int startlev = WarpX::load_balance_startlevel; + for (int lev = startlev; lev <= nLevels; ++lev) + { + int doLoadBalance = false; + + // Compute the new distribution mapping + DistributionMapping newdm; + const amrex::Real nboxes = costs[lev]->size(); + const amrex::Real nprocs = ParallelContext::NProcsSub(); + const int nmax = static_cast(std::ceil(nboxes/nprocs*load_balance_knapsack_factor)); + // These store efficiency (meaning, the average 'cost' over all ranks, + // normalized to max cost) for current and proposed distribution mappings + amrex::Real currentEfficiency = 0.0; + amrex::Real proposedEfficiency = 0.0; + + //if (lev == 0 || !do_similar_dm_refpatch) { + newdm = (load_balance_with_sfc) + ? DistributionMapping::makeSFC(*costs[lev], + currentEfficiency, proposedEfficiency, + false, + ParallelDescriptor::IOProcessorNumber()) + : DistributionMapping::makeKnapSack(*costs[lev], + currentEfficiency, proposedEfficiency, + nmax, + false, + ParallelDescriptor::IOProcessorNumber()); + //} else { + // amrex::BoxArray coarse_ba = boxArray(lev-1); + // amrex::DistributionMapping coarse_dm = DistributionMap(lev-1); + // amrex::BoxArray ba = boxArray(lev); + // ba.coarsen(WarpX::RefRatio(lev-1)); + // newdm = amrex::MakeSimilarDM(ba, coarse_ba, coarse_dm, getngEB()); + //} + // As specified in the above calls to makeSFC and makeKnapSack, the new + // distribution mapping is NOT communicated to all ranks; the loadbalanced + // dm is up-to-date only on root, and we can decide whether to broadcast + if ((load_balance_efficiency_ratio_threshold > 0.0) + && (ParallelDescriptor::MyProc() == ParallelDescriptor::IOProcessorNumber())) + { + doLoadBalance = (proposedEfficiency > load_balance_efficiency_ratio_threshold*currentEfficiency); + } + + ParallelDescriptor::Bcast(&doLoadBalance, 1, + ParallelDescriptor::IOProcessorNumber()); + + if (doLoadBalance) + { + Vector pmap; + if (ParallelDescriptor::MyProc() == ParallelDescriptor::IOProcessorNumber()) + { + pmap = newdm.ProcessorMap(); + } else + { + pmap.resize(static_cast(nboxes)); + } + ParallelDescriptor::Bcast(pmap.data(), pmap.size(), ParallelDescriptor::IOProcessorNumber()); + + if (ParallelDescriptor::MyProc() != ParallelDescriptor::IOProcessorNumber()) + { + newdm = DistributionMapping(pmap); + } + + RemakeLevel(lev, t_new[lev], boxArray(lev), newdm); + + // Record the load balance efficiency + setLoadBalanceEfficiency(lev, proposedEfficiency); + } + + loadBalancedAnyLevel = loadBalancedAnyLevel || doLoadBalance; + } } if (loadBalancedAnyLevel) { diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index f1900257035..b7410920edb 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -989,6 +989,26 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int rrfac = m_gdb->refRatio(0); fine_injection_box.coarsen(rrfac); } + bool refineplasma = false; + amrex::ParticleLocator > refinepatch_locator; + amrex::ParticleLocator > parent_locator; + parent_locator.build(ParticleBoxArray(0),Geom(0)); + if (WarpX::refineAddplasma) + { + refineplasma = true; + rrfac = WarpX::AddplasmaRefRatio; + auto fineba = ParticleBoxArray(1); + auto coarsened_fineba = fineba; + coarsened_fineba.coarsen(m_gdb->refRatio(0)); + if (!refinepatch_locator.isValid(coarsened_fineba)) { + refinepatch_locator.build(coarsened_fineba, Geom(0)); + } + refinepatch_locator.setGeometry(Geom(0)); + } + auto assignpartgrid = (WarpX::refineAddplasma) ? refinepatch_locator.getGridAssignor() + : parent_locator.getGridAssignor(); + // if assign_grid(ijk_vec) > 0, then we are in refinement patch. therefore refine plasma particles + // else, usual num_part InjectorPosition* inj_pos = plasma_injector.getInjectorPosition(); InjectorDensity* inj_rho = plasma_injector.getInjectorDensity(); @@ -1094,11 +1114,15 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int lo.z = applyBallisticCorrection(lo, inj_mom, gamma_boost, beta_boost, t); hi.z = applyBallisticCorrection(hi, inj_mom, gamma_boost, beta_boost, t); - if (inj_pos->overlapsWith(lo, hi)) { auto index = overlap_box.index(iv); - const amrex::Long r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv))? + amrex::IntVect glo_iv = iv + tile_box.smallEnd(); + bool in_refpatch = false; + if ( !(assignpartgrid(glo_iv)<0) && refineplasma) { + in_refpatch = true; + } + const amrex::Long r = ( (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) || (refineplasma && in_refpatch) ) ? (AMREX_D_TERM(lrrfac[0],*lrrfac[1],*lrrfac[2])) : (1); pcounts[index] = num_ppc*r; // update pcount by checking if cell-corners or cell-center @@ -1137,7 +1161,6 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int // Max number of new particles. All of them are created, // and invalid ones are then discarded const amrex::Long max_new_particles = Scan::ExclusiveSum(counts.size(), counts.data(), offset.data()); - // Update NextID to include particles created in this function amrex::Long pid; #ifdef AMREX_USE_OMP @@ -1256,6 +1279,11 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int { const IntVect iv = IntVect(AMREX_D_DECL(i, j, k)); const auto index = overlap_box.index(iv); + amrex::IntVect glo_iv = iv + tile_box.smallEnd(); + bool in_refpatch = false; + if ( !(assignpartgrid(glo_iv)<0) && refineplasma) { + in_refpatch = true; + } #ifdef WARPX_DIM_RZ Real theta_offset = 0._rt; if (rz_random_theta) { theta_offset = amrex::Random(engine) * 2._rt * MathConst::pi; } @@ -1276,13 +1304,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int { long ip = poffset[index] + i_part; pa_idcpu[ip] = amrex::SetParticleIDandCPU(pid+ip, cpuid); - const XDim3 r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ? + const XDim3 r = ( (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) || (refineplasma && in_refpatch)) ? // In the refined injection region: use refinement ratio `lrrfac` inj_pos->getPositionUnitBox(i_part, lrrfac, engine) : // Otherwise: use 1 as the refinement ratio inj_pos->getPositionUnitBox(i_part, amrex::IntVect::TheUnitVector(), engine); auto pos = getCellCoords(overlap_corner, dx, r, iv); - #if defined(WARPX_DIM_3D) if (!tile_realbox.contains(XDim3{pos.x,pos.y,pos.z})) { ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi @@ -2028,12 +2055,10 @@ PhysicalParticleContainer::Evolve (int lev, WARPX_PROFILE_VAR_NS("PhysicalParticleContainer::Evolve::GatherAndPush", blp_fg); BL_ASSERT(OnSameGrids(lev,jx)); - amrex::LayoutData* cost = WarpX::getCosts(lev); const iMultiFab* current_masks = WarpX::CurrentBufferMasks(lev); const iMultiFab* gather_masks = WarpX::GatherBufferMasks(lev); - const bool has_buffer = cEx || cjx; if (m_do_back_transformed_particles) @@ -2062,7 +2087,8 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox filtered_Ex, filtered_Ey, filtered_Ez; FArrayBox filtered_Bx, filtered_By, filtered_Bz; - + FArrayBox bufferEx, bufferEy, bufferEz; + FArrayBox bufferBx, bufferBy, bufferBz; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) @@ -2132,7 +2158,7 @@ PhysicalParticleContainer::Evolve (int lev, DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); - if (has_buffer){ + if ((np-np_current)> 0 ){ DepositCharge(pti, wp, ion_lev, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); } @@ -2177,6 +2203,7 @@ PhysicalParticleContainer::Evolve (int lev, if (WarpX::use_fdtd_nci_corr) { + // should this be bufferEFields* // Filter arrays (*cEx)[pti], store the result in // filtered_Ex and update pointer cexfab so that it // points to filtered_Ex (and do the same for all @@ -2188,7 +2215,6 @@ PhysicalParticleContainer::Evolve (int lev, (*cBx)[pti], (*cBy)[pti], (*cBz)[pti], cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab); } - // Field gather and push for particles in gather buffers e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); if (push_type == PushType::Explicit) { @@ -2205,7 +2231,6 @@ PhysicalParticleContainer::Evolve (int lev, lev, lev-1, dt, ScaleFields(false), a_dt_type); } } - WARPX_PROFILE_VAR_STOP(blp_fg); // Current Deposition @@ -2226,7 +2251,7 @@ PhysicalParticleContainer::Evolve (int lev, 0, np_current, thread_num, lev, lev, dt / WarpX::n_subcycle_current, relative_time, push_type); - if (has_buffer) + if ((np-np_current)>0) { // Deposit in buffers DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz, @@ -2249,7 +2274,7 @@ PhysicalParticleContainer::Evolve (int lev, DepositCharge(pti, wp, ion_lev, rho, 1, 0, np_current, thread_num, lev, lev); - if (has_buffer){ + if ((np-np_current)>0 ){ DepositCharge(pti, wp, ion_lev, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); } diff --git a/Source/Particles/Sorting/Partition.cpp b/Source/Particles/Sorting/Partition.cpp index 9f92739a5d5..c94c5e185e9 100644 --- a/Source/Particles/Sorting/Partition.cpp +++ b/Source/Particles/Sorting/Partition.cpp @@ -86,7 +86,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers( // Second, among particles that are in the larger buffer, partition // particles into the smaller buffer - if (WarpX::n_current_deposition_buffer == WarpX::n_field_gather_buffer) { + if (WarpX::n_current_deposition_buffer == WarpX::n_field_gather_buffer ) { // No need to do anything if the buffers have the same size nfine_current = nfine_gather = iteratorDistance(pid.begin(), sep); } else if (sep == pid.end()) { diff --git a/Source/WarpX.H b/Source/WarpX.H index 3361869f4c8..39dc429698b 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -350,6 +350,8 @@ public: static bool do_dynamic_scheduling; static bool refine_plasma; + static bool refineAddplasma; + static amrex::IntVect AddplasmaRefRatio; static utils::parser::IntervalsParser sort_intervals; static amrex::IntVect sort_bin_size; @@ -370,10 +372,13 @@ public: //! #n_field_gather_buffer cells of the edge of the patch, will gather the fields //! from the lower refinement level instead of the refinement patch itself static int n_field_gather_buffer; + static amrex::IntVect n_field_gather_buffer_each_dir; + //! With mesh refinement, particles located inside a refinement patch, but within //! #n_current_deposition_buffer cells of the edge of the patch, will deposit their charge //! and current onto the lower refinement level instead of the refinement patch itself static int n_current_deposition_buffer; + static amrex::IntVect n_current_deposition_buffer_each_dir; //! Integer that corresponds to the type of grid used in the simulation //! (collocated, staggered, hybrid) @@ -699,7 +704,7 @@ public: { return load_balance_intervals; } - + static int load_balance_startlevel; /** * \brief Private function for spectral solver * Applies a damping factor in the guards cells that extend @@ -909,6 +914,10 @@ public: static amrex::XDim3 InvCellSize (int lev); static amrex::RealBox getRealBox(const amrex::Box& bx, int lev); +// [[nodiscard]] const amrex::MultiFab* getInterpWeightGBuffer (int lev) const +// { +// return interp_weight_gbuffer[lev].get(); +// } /** * \brief Return the lower corner of the box in real units. * \param bx The box @@ -1117,7 +1126,8 @@ public: // for cuda void BuildBufferMasksInBox ( amrex::Box tbx, amrex::IArrayBox &buffer_mask, - const amrex::IArrayBox &guard_mask, int ng ); + const amrex::IArrayBox &guard_mask, const int ng ); + #ifdef AMREX_USE_EB amrex::EBFArrayBoxFactory const& fieldEBFactory (int lev) const noexcept { return static_cast(*m_field_factory[lev]); @@ -1422,6 +1432,7 @@ private: return gather_buffer_masks[lev].get(); } + /** * \brief Re-orders the Fornberg coefficients so that they can be used more conveniently for * finite-order centering operations. For example, for finite-order centering of order 6, @@ -1591,6 +1602,11 @@ private: amrex::Vector, 3 > > Bfield_cax; amrex::Vector > current_buffer_masks; amrex::Vector > gather_buffer_masks; +// /** Multifab that stores weights for interpolating fine and coarse solutions +// * in the buffer gather region, such that, the weight for fine solution is +// * 0 near the coarse-fine interface, and 1 as the buffer gather region terminates in the fine patch +// */ +// amrex::Vector > interp_weight_gbuffer; // If charge/current deposition buffers are used amrex::Vector, 3 > > current_buf; @@ -1618,9 +1634,12 @@ private: int pml_has_particles = 0; int do_pml_j_damping = 0; int do_pml_in_domain = 0; - bool do_similar_dm_pml = true; + bool do_similar_dm_pml = false; + int do_cubic_sigma_pml = 0; + amrex::Real pml_damping_strength = 4.; bool do_pml_dive_cleaning; // default set in WarpX.cpp bool do_pml_divb_cleaning; // default set in WarpX.cpp + static bool do_abc_in_pml; amrex::Vector do_pml_Lo; amrex::Vector do_pml_Hi; amrex::Vector > pml; @@ -1681,6 +1700,10 @@ private: * uniform plasma on a domain of size 128 by 128 by 128, from which the approximate * time per iteration per particle is computed. */ amrex::Real costs_heuristic_particles_wt = amrex::Real(0); + /** Whether to use MakeSimilarDM when load balancing lev > 1 */ + int do_similar_dm_refpatch = 0; + /** Whether to use SFC on vectorized BoxArrays from all levels */ + int do_SFC_dm_vectorlevel = 0; // Determines timesteps for override sync utils::parser::IntervalsParser override_sync_intervals; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index a9aacdf1fc6..cf1aa908202 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -128,6 +128,8 @@ int WarpX::macroscopic_solver_algo; bool WarpX::do_single_precision_comms = false; bool WarpX::do_subcycle_current = false; int WarpX::n_subcycle_current = 1; +bool WarpX::do_abc_in_pml = false; +int WarpX::load_balance_startlevel = 0; bool WarpX::do_shared_mem_charge_deposition = false; bool WarpX::do_shared_mem_current_deposition = false; @@ -175,6 +177,8 @@ bool WarpX::use_filter_compensation = false; bool WarpX::serialize_initial_conditions = false; bool WarpX::refine_plasma = false; +bool WarpX::refineAddplasma = false; +amrex::IntVect WarpX::AddplasmaRefRatio(AMREX_D_DECL(1,1,1)); int WarpX::num_mirrors = 0; @@ -208,6 +212,8 @@ std::map WarpX::imultifab_map; IntVect WarpX::filter_npass_each_dir(1); +amrex::IntVect WarpX::n_field_gather_buffer_each_dir(-1); +amrex::IntVect WarpX::n_current_deposition_buffer_each_dir(-1); int WarpX::n_field_gather_buffer = -1; int WarpX::n_current_deposition_buffer = -1; @@ -396,6 +402,7 @@ WarpX::WarpX () gather_buffer_masks.resize(nlevs_max); current_buf.resize(nlevs_max); charge_buf.resize(nlevs_max); +// interp_weight_gbuffer.resize(nlevs_max); pml.resize(nlevs_max); #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_FFT) @@ -898,12 +905,42 @@ WarpX::ReadParameters () pp_warpx.query("serialize_initial_conditions", serialize_initial_conditions); pp_warpx.query("refine_plasma", refine_plasma); + pp_warpx.query("refineAddplasma", refineAddplasma); + amrex::Vector addplasma_ref_ratio(AMREX_SPACEDIM,1); + const bool addplasma_ref_ratio_specified = + utils::parser::queryArrWithParser( + pp_warpx, "refineplasma_refratio", + addplasma_ref_ratio, 0, AMREX_SPACEDIM); + if (addplasma_ref_ratio_specified){ + amrex::Print() << " refratio for addplasma is : "; + for (int i=0; i 0) AddplasmaRefRatio = RefRatio(0); + } + pp_warpx.query("do_dive_cleaning", do_dive_cleaning); pp_warpx.query("do_divb_cleaning", do_divb_cleaning); utils::parser::queryWithParser( pp_warpx, "n_field_gather_buffer", n_field_gather_buffer); + amrex::Vector nfieldgatherbuffer_eachdir(AMREX_SPACEDIM,n_field_gather_buffer); + utils::parser::queryArrWithParser( + pp_warpx, "n_field_gather_buffer_each_dir", nfieldgatherbuffer_eachdir, 0, AMREX_SPACEDIM); + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + n_field_gather_buffer_each_dir[i] = nfieldgatherbuffer_eachdir[i]; + } utils::parser::queryWithParser( pp_warpx, "n_current_deposition_buffer", n_current_deposition_buffer); + amrex::Vector ncurrentdepositionbuffer_eachdir(AMREX_SPACEDIM,n_current_deposition_buffer); + utils::parser::queryArrWithParser( + pp_warpx, "n_current_deposition_buffer_each_dir", + ncurrentdepositionbuffer_eachdir, 0, AMREX_SPACEDIM); + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + n_current_deposition_buffer_each_dir[i] = ncurrentdepositionbuffer_eachdir[i]; + } amrex::Real quantum_xi_tmp; const auto quantum_xi_is_specified = @@ -944,6 +981,9 @@ WarpX::ReadParameters () pp_warpx.query("do_pml_j_damping", do_pml_j_damping); pp_warpx.query("do_pml_in_domain", do_pml_in_domain); pp_warpx.query("do_similar_dm_pml", do_similar_dm_pml); + pp_warpx.query("do_cubic_sigma_pml",do_cubic_sigma_pml); + pp_warpx.query("pml_damping_strength",pml_damping_strength); + pp_warpx.query("do_abc_in_pml",do_abc_in_pml); // Read `v_particle_pml` in units of the speed of light v_particle_pml = 1._rt; utils::parser::queryWithParser(pp_warpx, "v_particle_pml", v_particle_pml); @@ -1362,6 +1402,9 @@ WarpX::ReadParameters () load_balance_intervals = utils::parser::IntervalsParser( load_balance_intervals_string_vec); pp_algo.query("load_balance_with_sfc", load_balance_with_sfc); + pp_algo.query("do_similar_dm_refpatch", do_similar_dm_refpatch); + pp_algo.query("do_SFC_dm_vectorlevel",do_SFC_dm_vectorlevel); + pp_algo.query("load_balance_startlevel",load_balance_startlevel); // Knapsack factor only used with non-SFC strategy if (!load_balance_with_sfc) { pp_algo.query("load_balance_knapsack_factor", load_balance_knapsack_factor); @@ -2117,6 +2160,7 @@ WarpX::ClearLevel (int lev) current_buffer_masks[lev].reset(); gather_buffer_masks[lev].reset(); +// interp_weight_gbuffer[lev].reset(); F_fp [lev].reset(); G_fp [lev].reset(); @@ -2202,7 +2246,16 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d // Field gather buffer should be larger than current deposition buffers n_field_gather_buffer = n_current_deposition_buffer + 1; } - + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + if (n_current_deposition_buffer_each_dir[i] < 0) { + n_current_deposition_buffer_each_dir[i] = guard_cells.ng_alloc_J.max(); + } + } + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + if (n_field_gather_buffer_each_dir[i] < 0) { + n_field_gather_buffer_each_dir[i] = n_current_deposition_buffer_each_dir[i] + 1; + } + } AllocLevelMFs(lev, ba, dm, guard_cells.ng_alloc_EB, guard_cells.ng_alloc_J, guard_cells.ng_alloc_Rho, guard_cells.ng_alloc_F, guard_cells.ng_alloc_G, aux_is_nodal); @@ -2746,7 +2799,8 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm // // Copy of the coarse aux // - if (lev > 0 && (n_field_gather_buffer > 0 || n_current_deposition_buffer > 0 || + if (lev > 0 && (n_field_gather_buffer > 0 || + n_current_deposition_buffer > 0 || mypc->nSpeciesGatherFromMainGrid() > 0)) { BoxArray cba = ba; @@ -2772,20 +2826,23 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm AllocInitMultiFab(Efield_cax[lev][1], amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngEB,lev, "Efield_cax[y]", 0.0_rt); AllocInitMultiFab(Efield_cax[lev][2], amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngEB,lev, "Efield_cax[z]", 0.0_rt); } - - AllocInitMultiFab(gather_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "gather_buffer_masks"); + amrex::Print() << " ba for allocating gahter buffer mask " << "\n"; // Gather buffer masks have 1 ghost cell, because of the fact // that particles may move by more than one cell when using subcycling. + AllocInitMultiFab(gather_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "gather_buffer_masks"); +// AllocInitMultiFab(interp_weight_gbuffer[lev], amrex::convert(ba, IntVect::TheNodeVector()), dm, ncomps, ngEB, lev, "interp_weight_gbuffer", 0.0_rt); } if (n_current_deposition_buffer > 0) { + amrex::Print() << " current dep buffer : " << n_current_deposition_buffer << "\n"; AllocInitMultiFab(current_buf[lev][0], amrex::convert(cba,jx_nodal_flag),dm,ncomps,ngJ,lev, "current_buf[x]"); AllocInitMultiFab(current_buf[lev][1], amrex::convert(cba,jy_nodal_flag),dm,ncomps,ngJ,lev, "current_buf[y]"); AllocInitMultiFab(current_buf[lev][2], amrex::convert(cba,jz_nodal_flag),dm,ncomps,ngJ,lev, "current_buf[z]"); if (rho_cp[lev]) { AllocInitMultiFab(charge_buf[lev], amrex::convert(cba,rho_nodal_flag),dm,2*ncomps,ngRho,lev, "charge_buf"); } - AllocInitMultiFab(current_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "current_buffer_masks"); + amrex::Print() << " allocate current buffer mask \n"; + AllocInitMultiFab(current_buffer_masks[lev], ba, dm, ncomps, ngJ, lev, "current_buffer_masks"); // Current buffer masks have 1 ghost cell, because of the fact // that particles may move by more than one cell when using subcycling. }