From 1f04d310001f76ca15d0409de20d2250c9abadf3 Mon Sep 17 00:00:00 2001 From: RevathiJambunathan Date: Sat, 2 Mar 2024 14:07:21 -0800 Subject: [PATCH] add pml parameters for ref patch : abcinpml, cubicsigma, dampingstrength --- Source/BoundaryConditions/PML.H | 22 ++- Source/BoundaryConditions/PML.cpp | 86 +++++++---- Source/BoundaryConditions/WarpXEvolvePML.cpp | 14 +- Source/BoundaryConditions/WarpX_PML_kernels.H | 146 ++++++++++++++++-- Source/Initialization/WarpXInitData.cpp | 2 + Source/WarpX.H | 3 + Source/WarpX.cpp | 4 + 7 files changed, 220 insertions(+), 57 deletions(-) diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 5f41ecbd706..2cba95f779b 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -44,15 +44,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); @@ -79,8 +82,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; @@ -93,7 +97,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 @@ -114,6 +118,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 @@ -123,7 +129,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: @@ -147,6 +154,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 805b4fec181..611f67485dc 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -61,11 +61,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(); @@ -75,14 +75,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; } }); } @@ -90,10 +98,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(); @@ -104,12 +113,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; } }); } @@ -143,7 +160,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()); @@ -179,22 +197,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(); @@ -208,7 +231,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) @@ -228,7 +252,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); } } @@ -236,7 +260,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); @@ -306,7 +331,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]); @@ -319,7 +344,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( @@ -355,7 +380,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]); @@ -364,7 +389,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( @@ -405,7 +430,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]); @@ -414,7 +439,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( @@ -505,9 +530,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 @@ -552,6 +577,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri 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), @@ -737,7 +763,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri 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_PSATD @@ -851,7 +877,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri 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_PSATD diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp index af721d70b6d..5952d59e8df 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/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index e32f5e50798..619a72514a5 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -579,6 +579,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 @@ -617,6 +618,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]); } } diff --git a/Source/WarpX.H b/Source/WarpX.H index 79352274068..63f22742e14 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1587,8 +1587,11 @@ private: int do_pml_j_damping = 0; int do_pml_in_domain = 0; static int do_similar_dm_pml; + 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; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 1136e29b739..00faf46e9e2 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -126,6 +126,7 @@ 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; bool WarpX::do_shared_mem_charge_deposition = false; bool WarpX::do_shared_mem_current_deposition = false; @@ -915,6 +916,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);