From c71681262d21edca6807e144e6b90306fff22d63 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Fri, 30 Oct 2026 14:36:11 -0700 Subject: [PATCH 1/3] CPML first cut --- Docs/source/usage/parameters.rst | 39 +++ Source/BoundaryConditions/PML.H | 35 +++ Source/BoundaryConditions/PML.cpp | 253 ++++++++++++++---- Source/BoundaryConditions/WarpXEvolvePML.cpp | 4 +- Source/Evolve/WarpXEvolve.cpp | 8 +- .../FiniteDifferenceSolver/CMakeLists.txt | 1 + .../FiniteDifferenceSolver.H | 54 ++++ .../FiniteDifferenceSolver/Make.package | 1 + Source/FieldSolver/WarpXExternalEMFields.cpp | 8 + Source/FieldSolver/WarpXPushFieldsEM.cpp | 131 +++++++-- Source/Utils/WarpXAlgorithmSelection.H | 13 + Source/WarpX.H | 23 +- Source/WarpX.cpp | 31 +++ 13 files changed, 529 insertions(+), 72 deletions(-) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index e167679e6..7427ad612 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -381,6 +381,45 @@ Additional PML parameters The characteristic depth, in number of cells, over which the absorption coefficients of the PML increases. +* ``warpx.pml_type`` (`string`; default: ``split``) + PML formulation used with the FDTD solver. + + * ``split``: split-field (Berenger) PML (original implementation). + + * ``convolutional`` (or ``cpml``): convolutional PML (Roden & Gedney, 2000), using the + complex-frequency-shifted coordinate stretching + :math:`s_i = \kappa_i + \sigma_i/(\alpha_i + i\omega)`. The total (unsplit) field is + advanced in the PML with the *same* material update as in the interior (local + :math:`\epsilon`, :math:`\mu`, :math:`\sigma`, including spatially varying properties and + the ``laxwendroff``/``backwardeuler`` macroscopic schemes); only the spatial derivatives are + stretched, through auxiliary convolution variables. The PEC mask built from the GDS/numpy + layer (``algo.use_PEC_mask``) and its extension into the PML is applied to the CPML fields + in the same way as for the split-field PML. The :math:`\sigma` profile is the same as for the + split-field PML (set by ``warpx.pml_ncell`` and ``warpx.pml_delta``). + Only implemented for the (staggered) Yee solver, without divergence cleaning in the PML, + J-damping, or moving window. + With the default parameters (``kappa_max = 1``, ``alpha_max = 0``) its absorption is + essentially the same as the split-field PML. Using ``warpx.pml_delta`` around 15-20 together + with a small ``warpx.cpml_alpha_max`` (and optionally ``warpx.cpml_kappa_max = 2``) + typically lowers the reflection further; a too-large ``alpha_max`` degrades absorption + of the low frequencies. + +* ``warpx.cpml_kappa_max`` (`float`; default: 1) + CPML only. Maximum of the real stretching factor :math:`\kappa`, graded as + :math:`\kappa(d) = 1 + (\kappa_{max}-1)\,d^m` with :math:`d \in [0,1]` the normalized + depth into the PML. Values > 1 improve absorption of evanescent waves. + +* ``warpx.cpml_alpha_max`` (`float`, in 1/s; default: 0) + CPML only. Maximum of the complex-frequency-shift parameter, given as the rate + :math:`\alpha_{SI}/\epsilon_0`, graded as :math:`\alpha(d) = \alpha_{max} (1-d)^{m_a}`. + A typical choice is :math:`2\pi f_{min}`, with :math:`f_{min}` the lowest frequency of interest. + +* ``warpx.cpml_grading_order`` (`float`; default: 2) + CPML only. Polynomial order :math:`m` of the :math:`\kappa` grading. + +* ``warpx.cpml_alpha_grading_order`` (`float`; default: 1) + CPML only. Polynomial order :math:`m_a` of the :math:`\alpha` grading. + * ``warpx.do_pml_in_domain`` (`int`; default: 0) Whether to create the PML inside the simulation area or outside. If inside, it allows the user to propagate particles in PML and to use extended PML diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 512244f45..b13dd23bf 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -70,6 +70,20 @@ struct SigmaBox SigmaVect sigma_star_cumsum_fac; amrex::Real v_sigma; + // Convolutional PML (CPML) coordinate-stretching parameters, + // s_i = kappa_i + sigma_i/(alpha_i + j*omega), using the same staggering as + // sigma (nodal) and sigma_star (cell-centered). sigma and alpha are rates (1/s), + // i.e. sigma_SI/eps0 and alpha_SI/eps0. + SigmaVect kappa; + SigmaVect kappa_star; + SigmaVect alpha; + SigmaVect alpha_star; + amrex::IntVect m_ncell; + amrex::Real m_cpml_kappa_max = 1.0; + amrex::Real m_cpml_alpha_max = 0.0; + amrex::Real m_cpml_grading_order = 2.0; + amrex::Real m_cpml_alpha_grading_order = 1.0; + }; class SigmaBoxFactory @@ -163,6 +177,20 @@ public: amrex::MultiFab* GetE_fp (int comp); amrex::MultiFab* GetPEC_fp (int comp); + /** Convolutional PML (warpx.pml_type = convolutional) auxiliary fields. + * psi_E[d] lives at the location of E_d and stores, in its two components, + * the recursive convolutions of the two curl(H) derivatives entering the + * update of E_d (same component ordering as PMLComp, e.g. Ex: xy=0 (d/dy), xz=1 (d/dz)). + * psi_B[d] is the analog for the B (or H, with LLG) update, using curl(E). + * These return nullptr with the split-field PML. */ + std::array GetPsiE_fp (); + std::array GetPsiB_fp (); + std::array GetPsiE_cp (); + std::array GetPsiB_cp (); + + /** True if this PML uses the convolutional (CPML) formulation */ + bool isConvolutional () const { return m_convolutional; } + // Used when WarpX::do_pml_dive_cleaning = true amrex::MultiFab* GetF_fp (); amrex::MultiFab* GetF_cp (); @@ -252,6 +280,13 @@ private: std::array,3> pml_PEC_fp; + // CPML auxiliary (memory) variables + bool m_convolutional = false; + std::array,3> pml_psiE_fp; + std::array,3> pml_psiB_fp; + std::array,3> pml_psiE_cp; + std::array,3> pml_psiB_cp; + std::array,3> pml_edge_lengths; std::array,3> pml_E_cp; diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 9c36950fe..a76ed5d09 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -60,19 +60,47 @@ using namespace amrex; namespace { - static void FillLo (Sigma& sigma, Sigma& sigma_cumsum, - Sigma& sigma_star, Sigma& sigma_star_cumsum, + // Grading of the CPML stretching parameters kappa and alpha, as a function of the + // normalized depth d in [0,1] into the PML (d=0 at the PML/domain interface, + // d=1 at the outer edge of the PML). + // kappa(d) = 1 + (kappa_max - 1) d^m, alpha(d) = alpha_max (1 - d)^ma + // (Roden & Gedney, 2000). These arrays are only used by the convolutional PML; + // the split-field PML only uses sigma and sigma_star. + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + Real cpml_kappa (Real d, Real kappa_max, Real m) noexcept + { + d = amrex::max(Real(0.), amrex::min(Real(1.), d)); + return Real(1.) + (kappa_max - Real(1.))*std::pow(d, m); + } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + Real cpml_alpha (Real d, Real alpha_max, Real ma) noexcept + { + d = amrex::max(Real(0.), amrex::min(Real(1.), d)); + return alpha_max*std::pow(Real(1.) - d, ma); + } + + static void FillLo (SigmaBox& sb, const int idim, const int olo, const int ohi, const int glo, Real fac, const amrex::Real v_sigma) { - const int slo = sigma.m_lo; - const int sslo = sigma_star.m_lo; + const int slo = sb.sigma[idim].m_lo; + const int sslo = sb.sigma_star[idim].m_lo; const int N = ohi+1-olo+1; - Real* p_sigma = sigma.data(); - Real* p_sigma_cumsum = sigma_cumsum.data(); - Real* p_sigma_star = sigma_star.data(); - Real* p_sigma_star_cumsum = sigma_star_cumsum.data(); + Real* p_sigma = sb.sigma[idim].data(); + Real* p_sigma_cumsum = sb.sigma_cumsum[idim].data(); + Real* p_sigma_star = sb.sigma_star[idim].data(); + Real* p_sigma_star_cumsum = sb.sigma_star_cumsum[idim].data(); + Real* p_kappa = sb.kappa[idim].data(); + Real* p_kappa_star = sb.kappa_star[idim].data(); + Real* p_alpha = sb.alpha[idim].data(); + Real* p_alpha_star = sb.alpha_star[idim].data(); + const Real inv_ncell = Real(1.)/static_cast(sb.m_ncell[idim]); + const Real kmax = sb.m_cpml_kappa_max; + const Real amax = sb.m_cpml_alpha_max; + const Real m = sb.m_cpml_grading_order; + const Real ma = sb.m_cpml_alpha_grading_order; amrex::ParallelFor(N, [=] AMREX_GPU_DEVICE (int i) noexcept { i += olo; @@ -80,63 +108,86 @@ namespace p_sigma[i-slo] = fac*(offset*offset); // 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_kappa[i-slo] = cpml_kappa(offset*inv_ncell, kmax, m); + p_alpha[i-slo] = cpml_alpha(offset*inv_ncell, amax, ma); if (i <= ohi+1) { offset = static_cast(glo-i) - 0.5_rt; p_sigma_star[i-sslo] = fac*(offset*offset); // 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_kappa_star[i-sslo] = cpml_kappa(offset*inv_ncell, kmax, m); + p_alpha_star[i-sslo] = cpml_alpha(offset*inv_ncell, amax, ma); } }); } - static void FillHi (Sigma& sigma, Sigma& sigma_cumsum, - Sigma& sigma_star, Sigma& sigma_star_cumsum, + static void FillHi (SigmaBox& sb, const int idim, const int olo, const int ohi, const int ghi, Real fac, const amrex::Real v_sigma) { - const int slo = sigma.m_lo; - const int sslo = sigma_star.m_lo; + const int slo = sb.sigma[idim].m_lo; + const int sslo = sb.sigma_star[idim].m_lo; const int N = ohi+1-olo+1; - Real* p_sigma = sigma.data(); - Real* p_sigma_cumsum = sigma_cumsum.data(); - Real* p_sigma_star = sigma_star.data(); - Real* p_sigma_star_cumsum = sigma_star_cumsum.data(); + Real* p_sigma = sb.sigma[idim].data(); + Real* p_sigma_cumsum = sb.sigma_cumsum[idim].data(); + Real* p_sigma_star = sb.sigma_star[idim].data(); + Real* p_sigma_star_cumsum = sb.sigma_star_cumsum[idim].data(); + Real* p_kappa = sb.kappa[idim].data(); + Real* p_kappa_star = sb.kappa_star[idim].data(); + Real* p_alpha = sb.alpha[idim].data(); + Real* p_alpha_star = sb.alpha_star[idim].data(); + const Real inv_ncell = Real(1.)/static_cast(sb.m_ncell[idim]); + const Real kmax = sb.m_cpml_kappa_max; + const Real amax = sb.m_cpml_alpha_max; + const Real m = sb.m_cpml_grading_order; + const Real ma = sb.m_cpml_alpha_grading_order; amrex::ParallelFor(N, [=] AMREX_GPU_DEVICE (int i) noexcept { 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; + p_kappa[i-slo] = cpml_kappa(offset*inv_ncell, kmax, m); + p_alpha[i-slo] = cpml_alpha(offset*inv_ncell, amax, ma); 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; + p_kappa_star[i-sslo] = cpml_kappa(offset*inv_ncell, kmax, m); + p_alpha_star[i-sslo] = cpml_alpha(offset*inv_ncell, amax, ma); } }); } #if (AMREX_SPACEDIM != 1) - static void FillZero (Sigma& sigma, Sigma& sigma_cumsum, - Sigma& sigma_star, Sigma& sigma_star_cumsum, + static void FillZero (SigmaBox& sb, const int idim, const int olo, const int ohi) { - const int slo = sigma.m_lo; - const int sslo = sigma_star.m_lo; + const int slo = sb.sigma[idim].m_lo; + const int sslo = sb.sigma_star[idim].m_lo; const int N = ohi+1-olo+1; - Real* p_sigma = sigma.data(); - Real* p_sigma_cumsum = sigma_cumsum.data(); - Real* p_sigma_star = sigma_star.data(); - Real* p_sigma_star_cumsum = sigma_star_cumsum.data(); + Real* p_sigma = sb.sigma[idim].data(); + Real* p_sigma_cumsum = sb.sigma_cumsum[idim].data(); + Real* p_sigma_star = sb.sigma_star[idim].data(); + Real* p_sigma_star_cumsum = sb.sigma_star_cumsum[idim].data(); + Real* p_kappa = sb.kappa[idim].data(); + Real* p_kappa_star = sb.kappa_star[idim].data(); + Real* p_alpha = sb.alpha[idim].data(); + Real* p_alpha_star = sb.alpha_star[idim].data(); amrex::ParallelFor(N, [=] AMREX_GPU_DEVICE (int i) noexcept { i += olo; p_sigma[i-slo] = Real(0.0); p_sigma_cumsum[i-slo] = Real(0.0); + p_kappa[i-slo] = Real(1.0); + p_alpha[i-slo] = Real(0.0); if (i <= ohi+1) { p_sigma_star[i-sslo] = Real(0.0); p_sigma_star_cumsum[i-sslo] = Real(0.0); + p_kappa_star[i-sslo] = Real(1.0); + p_alpha_star[i-sslo] = Real(0.0); } }); } @@ -180,8 +231,28 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, const sigma_star_fac [idim].m_hi = hi[idim]+1; sigma_star_cumsum_fac[idim].m_lo = lo[idim]; sigma_star_cumsum_fac[idim].m_hi = hi[idim]+1; + + // CPML stretching parameters (same staggering/indexing as sigma and sigma_star) + kappa [idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); + kappa_star [idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); + alpha [idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); + alpha_star [idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); + kappa [idim].m_lo = lo[idim]; + kappa [idim].m_hi = hi[idim]+1; + kappa_star [idim].m_lo = lo[idim]; + kappa_star [idim].m_hi = hi[idim]+1; + alpha [idim].m_lo = lo[idim]; + alpha [idim].m_hi = hi[idim]+1; + alpha_star [idim].m_lo = lo[idim]; + alpha_star [idim].m_hi = hi[idim]+1; } + m_ncell = ncell; + m_cpml_kappa_max = WarpX::cpml_kappa_max; + m_cpml_alpha_max = WarpX::cpml_alpha_max; + m_cpml_grading_order = WarpX::cpml_grading_order; + m_cpml_alpha_grading_order = WarpX::cpml_alpha_grading_order; + 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])); @@ -208,8 +279,7 @@ void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell, int olo = std::max(slo, dlo-ncell[idim]); int ohi = std::min(shi, dlo-1); if (ohi >= olo) { - FillLo(sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], + FillLo(*this, idim, olo, ohi, dlo, fac[idim], v_sigma_sb); } @@ -218,8 +288,7 @@ void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell, olo = std::max(slo, dlo); ohi = std::min(shi, dhi); if (ohi >= olo) { - FillZero(sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], + FillZero(*this, idim, olo, ohi); } #endif @@ -228,8 +297,7 @@ void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell, olo = std::max(slo, dhi+1); ohi = std::min(shi, dhi+ncell[idim]); if (ohi >= olo) { - FillHi(sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], + FillHi(*this, idim, olo, ohi, dhi, fac[idim], v_sigma_sb); } } @@ -306,8 +374,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int Box looverlap = lobox & box; if (looverlap.ok()) { - FillLo(sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], + FillLo(*this, idim, looverlap.smallEnd(idim), looverlap.bigEnd(idim), grid_box.smallEnd(idim), fac[idim], v_sigma_sb); } @@ -319,8 +386,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int #endif Box hioverlap = hibox & box; if (hioverlap.ok()) { - FillHi(sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], + FillHi(*this, idim, hioverlap.smallEnd(idim), hioverlap.bigEnd(idim), grid_box.bigEnd(idim), fac[idim], v_sigma_sb); } @@ -343,8 +409,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int "SigmaBox::SigmaBox(): side_side_edges, how did this happen?" ); - FillZero(sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], + FillZero(*this, idim, overlap.smallEnd(idim), overlap.bigEnd(idim)); } @@ -355,8 +420,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int Box lobox = amrex::adjCellLo(grid_box, idim, ncell[idim]); Box looverlap = lobox.grow(jdim,ncell[jdim]).grow(kdim,ncell[kdim]) & box; if (looverlap.ok()) { - FillLo(sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], + FillLo(*this, idim, looverlap.smallEnd(idim), looverlap.bigEnd(idim), grid_box.smallEnd(idim), fac[idim], v_sigma_sb); } @@ -364,8 +428,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int Box hibox = amrex::adjCellHi(grid_box, idim, ncell[idim]); Box hioverlap = hibox.grow(jdim,ncell[jdim]).grow(kdim,ncell[kdim]) & box; if (hioverlap.ok()) { - FillHi(sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], + FillHi(*this, idim, hioverlap.smallEnd(idim), hioverlap.bigEnd(idim), grid_box.bigEnd(idim), fac[idim], v_sigma_sb); } @@ -392,8 +455,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int "SigmaBox::SigmaBox(): side_faces, how did this happen?" ); - FillZero(sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], + FillZero(*this, idim, overlap.smallEnd(idim), overlap.bigEnd(idim)); } #endif @@ -405,8 +467,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int const Box& lobox = amrex::adjCellLo(grid_box, idim, ncell[idim]); Box looverlap = lobox & box; if (looverlap.ok()) { - FillLo(sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], + FillLo(*this, idim, looverlap.smallEnd(idim), looverlap.bigEnd(idim), grid_box.smallEnd(idim), fac[idim], v_sigma_sb); } @@ -414,8 +475,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell[idim]); Box hioverlap = hibox & box; if (hioverlap.ok()) { - FillHi(sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], + FillHi(*this, idim, hioverlap.smallEnd(idim), hioverlap.bigEnd(idim), grid_box.bigEnd(idim), fac[idim], v_sigma_sb); } @@ -596,6 +656,18 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri m_ok = true; } + m_convolutional = (WarpX::pml_type == PMLType::Convolutional); + if (m_convolutional) { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::Yee && + grid_type != GridType::Collocated, + "warpx.pml_type = convolutional is only implemented for the (staggered) Yee FDTD solver"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!do_pml_dive_cleaning && !do_pml_divb_cleaning, + "warpx.pml_type = convolutional does not support divergence cleaning in the PML"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!do_moving_window, + "warpx.pml_type = convolutional does not support the moving window"); + } + // Define the number of guard cells in each direction, for E, B, and F IntVect nge = IntVect(AMREX_D_DECL(2, 2, 2)); IntVect ngb = IntVect(AMREX_D_DECL(2, 2, 2)); @@ -696,6 +768,24 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri WarpX::AllocInitMultiFab(pml_H_fp[2], ba_Hz, dm, 2, ngb, "pml_H_fp[z]"); #endif + if (m_convolutional) { + // CPML memory variables: 2 components each, no guard cells needed + // (they are only read/written at the point where they are defined) + WarpX::AllocInitMultiFab(pml_psiE_fp[0], ba_Ex, dm, 2, IntVect(0), "pml_psiE_fp[x]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiE_fp[1], ba_Ey, dm, 2, IntVect(0), "pml_psiE_fp[y]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiE_fp[2], ba_Ez, dm, 2, IntVect(0), "pml_psiE_fp[z]", 0.0_rt); +#ifdef WARPX_MAG_LLG + // with LLG, H (not B) is advanced in the PML + WarpX::AllocInitMultiFab(pml_psiB_fp[0], ba_Hx, dm, 2, IntVect(0), "pml_psiB_fp[x]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiB_fp[1], ba_Hy, dm, 2, IntVect(0), "pml_psiB_fp[y]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiB_fp[2], ba_Hz, dm, 2, IntVect(0), "pml_psiB_fp[z]", 0.0_rt); +#else + WarpX::AllocInitMultiFab(pml_psiB_fp[0], ba_Bx, dm, 2, IntVect(0), "pml_psiB_fp[x]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiB_fp[1], ba_By, dm, 2, IntVect(0), "pml_psiB_fp[y]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiB_fp[2], ba_Bz, dm, 2, IntVect(0), "pml_psiB_fp[z]", 0.0_rt); +#endif + } + if (WarpX::em_solver_medium == MediumForEM::Macroscopic) { // Allocating macroproperties in pml at cell-centers pml_eps_fp = std::make_unique(ba, dm, 1, nge); @@ -937,6 +1027,21 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri WarpX::AllocInitMultiFab(pml_H_cp[2], cba_Hz, cdm, 2, ngb, "pml_H_cp[z]"); #endif + if (m_convolutional) { + WarpX::AllocInitMultiFab(pml_psiE_cp[0], cba_Ex, cdm, 2, IntVect(0), "pml_psiE_cp[x]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiE_cp[1], cba_Ey, cdm, 2, IntVect(0), "pml_psiE_cp[y]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiE_cp[2], cba_Ez, cdm, 2, IntVect(0), "pml_psiE_cp[z]", 0.0_rt); +#ifdef WARPX_MAG_LLG + WarpX::AllocInitMultiFab(pml_psiB_cp[0], cba_Hx, cdm, 2, IntVect(0), "pml_psiB_cp[x]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiB_cp[1], cba_Hy, cdm, 2, IntVect(0), "pml_psiB_cp[y]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiB_cp[2], cba_Hz, cdm, 2, IntVect(0), "pml_psiB_cp[z]", 0.0_rt); +#else + WarpX::AllocInitMultiFab(pml_psiB_cp[0], cba_Bx, cdm, 2, IntVect(0), "pml_psiB_cp[x]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiB_cp[1], cba_By, cdm, 2, IntVect(0), "pml_psiB_cp[y]", 0.0_rt); + WarpX::AllocInitMultiFab(pml_psiB_cp[2], cba_Bz, cdm, 2, IntVect(0), "pml_psiB_cp[z]", 0.0_rt); +#endif + } + // Allocating macroproperties in pml at cell-centers if (WarpX::em_solver_medium == MediumForEM::Macroscopic) { pml_eps_cp = std::make_unique(cba, dm, 1, nge); @@ -1193,6 +1298,30 @@ PML::GetPEC_fp (int comp) return pml_PEC_fp[comp].get(); } +std::array +PML::GetPsiE_fp () +{ + return {pml_psiE_fp[0].get(), pml_psiE_fp[1].get(), pml_psiE_fp[2].get()}; +} + +std::array +PML::GetPsiB_fp () +{ + return {pml_psiB_fp[0].get(), pml_psiB_fp[1].get(), pml_psiB_fp[2].get()}; +} + +std::array +PML::GetPsiE_cp () +{ + return {pml_psiE_cp[0].get(), pml_psiE_cp[1].get(), pml_psiE_cp[2].get()}; +} + +std::array +PML::GetPsiB_cp () +{ + return {pml_psiB_cp[0].get(), pml_psiB_cp[1].get(), pml_psiB_cp[2].get()}; +} + std::array PML::GetB_fp () { @@ -1647,6 +1776,14 @@ PML::CheckPoint (const std::string& dir) const VisMF::AsyncWrite(*pml_H_fp[1], dir+"_Hy_fp"); VisMF::AsyncWrite(*pml_H_fp[2], dir+"_Hz_fp"); #endif + if (m_convolutional) { + VisMF::AsyncWrite(*pml_psiE_fp[0], dir+"_psiEx_fp"); + VisMF::AsyncWrite(*pml_psiE_fp[1], dir+"_psiEy_fp"); + VisMF::AsyncWrite(*pml_psiE_fp[2], dir+"_psiEz_fp"); + VisMF::AsyncWrite(*pml_psiB_fp[0], dir+"_psiBx_fp"); + VisMF::AsyncWrite(*pml_psiB_fp[1], dir+"_psiBy_fp"); + VisMF::AsyncWrite(*pml_psiB_fp[2], dir+"_psiBz_fp"); + } } if (pml_E_cp[0]) @@ -1662,6 +1799,14 @@ PML::CheckPoint (const std::string& dir) const VisMF::AsyncWrite(*pml_H_cp[1], dir+"_Hy_cp"); VisMF::AsyncWrite(*pml_H_cp[2], dir+"_Hz_cp"); #endif + if (m_convolutional) { + VisMF::AsyncWrite(*pml_psiE_cp[0], dir+"_psiEx_cp"); + VisMF::AsyncWrite(*pml_psiE_cp[1], dir+"_psiEy_cp"); + VisMF::AsyncWrite(*pml_psiE_cp[2], dir+"_psiEz_cp"); + VisMF::AsyncWrite(*pml_psiB_cp[0], dir+"_psiBx_cp"); + VisMF::AsyncWrite(*pml_psiB_cp[1], dir+"_psiBy_cp"); + VisMF::AsyncWrite(*pml_psiB_cp[2], dir+"_psiBz_cp"); + } } } @@ -1681,6 +1826,14 @@ PML::Restart (const std::string& dir) VisMF::Read(*pml_H_fp[1], dir+"_Hy_fp"); VisMF::Read(*pml_H_fp[2], dir+"_Hz_fp"); #endif + if (m_convolutional) { + VisMF::Read(*pml_psiE_fp[0], dir+"_psiEx_fp"); + VisMF::Read(*pml_psiE_fp[1], dir+"_psiEy_fp"); + VisMF::Read(*pml_psiE_fp[2], dir+"_psiEz_fp"); + VisMF::Read(*pml_psiB_fp[0], dir+"_psiBx_fp"); + VisMF::Read(*pml_psiB_fp[1], dir+"_psiBy_fp"); + VisMF::Read(*pml_psiB_fp[2], dir+"_psiBz_fp"); + } } if (pml_E_cp[0]) @@ -1696,6 +1849,14 @@ PML::Restart (const std::string& dir) VisMF::Read(*pml_H_cp[1], dir+"_Hy_cp"); VisMF::Read(*pml_H_cp[2], dir+"_Hz_cp"); #endif + if (m_convolutional) { + VisMF::Read(*pml_psiE_cp[0], dir+"_psiEx_cp"); + VisMF::Read(*pml_psiE_cp[1], dir+"_psiEy_cp"); + VisMF::Read(*pml_psiE_cp[2], dir+"_psiEz_cp"); + VisMF::Read(*pml_psiB_cp[0], dir+"_psiBx_cp"); + VisMF::Read(*pml_psiB_cp[1], dir+"_psiBy_cp"); + VisMF::Read(*pml_psiB_cp[2], dir+"_psiBz_cp"); + } } } diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp index a09fd042e..cee4dd62b 100644 --- a/Source/BoundaryConditions/WarpXEvolvePML.cpp +++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp @@ -67,7 +67,9 @@ WarpX::DampPML (const int lev, PatchType patch_type) dt[lev]); } #endif - if (pml[lev]) { + // The convolutional PML absorbs through the stretched derivatives in the + // field updates; no separate damping step is applied + if (pml[lev] && !pml[lev]->isConvolutional()) { DampPML_Cartesian (lev, patch_type); } } diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 698ec4573..6303c0f71 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -497,9 +497,9 @@ WarpX::OneStep_nosub (Real cur_time) #ifndef WARPX_DIM_RZ if (WarpX::em_solver_medium == MediumForEM::Macroscopic) { //evolveM is not applicable to vacuum if (mag_time_scheme_order==1){ - MacroscopicEvolveHM(0.5*dt[0]); // we now have M^{n+1/2} and H^{n+1/2} + MacroscopicEvolveHM(0.5*dt[0], DtType::FirstHalf); // we now have M^{n+1/2} and H^{n+1/2} } else if (mag_time_scheme_order==2){ - MacroscopicEvolveHM_2nd(0.5*dt[0]); // we now have M^{n+1/2} and H^{n+1/2} + MacroscopicEvolveHM_2nd(0.5*dt[0], DtType::FirstHalf); // we now have M^{n+1/2} and H^{n+1/2} } else { amrex::Abort("unsupported mag_time_scheme_order for M field"); } @@ -557,9 +557,9 @@ WarpX::OneStep_nosub (Real cur_time) #ifndef WARPX_DIM_RZ if (WarpX::em_solver_medium == MediumForEM::Macroscopic) { if (mag_time_scheme_order==1){ - MacroscopicEvolveHM(0.5*dt[0]); // we now have M^{n+1} and H^{n+1} + MacroscopicEvolveHM(0.5*dt[0], DtType::SecondHalf); // we now have M^{n+1} and H^{n+1} } else if (mag_time_scheme_order==2){ - MacroscopicEvolveHM_2nd(0.5*dt[0]); // we now have M^{n+1} and H^{n+1} + MacroscopicEvolveHM_2nd(0.5*dt[0], DtType::SecondHalf); // we now have M^{n+1} and H^{n+1} } else { amrex::Abort("unsupported mag_time_scheme_order for M field"); } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt index c95e646ba..4d68ec0b1 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt +++ b/Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt @@ -3,6 +3,7 @@ target_sources(WarpX ComputeDivE.cpp EvolveB.cpp EvolveBPML.cpp + EvolveCPML.cpp EvolveE.cpp EvolveEPML.cpp EvolveF.cpp diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index c1d24d74c..498105709 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -12,6 +12,7 @@ #include "FiniteDifferenceSolver_fwd.H" #include "BoundaryConditions/PML_fwd.H" +#include "Evolve/WarpXDtType.H" #include "MacroscopicProperties/MacroscopicProperties_fwd.H" #include @@ -198,6 +199,46 @@ class FiniteDifferenceSolver amrex::MultiFab* const mu_mf, amrex::MultiFab* const sigma_mf); + /** Convolutional PML (warpx.pml_type = convolutional) updates (Yee only). + * The total field is stored in component 0 of the PML MultiFabs; psi* are + * the CPML memory variables (see BoundaryConditions/CPML_K.H). */ + /** B is advanced in two half steps per time step (FirstHalf/SecondHalf). + * To be consistent with leapfrog, the memory variable psiB is updated once + * per full step (on the SecondHalf, with the coefficients of the full step) + * and the same psiB is used by both halves. */ + void EvolveBCPML ( std::array< amrex::MultiFab*, 3 > Bfield, + std::array< amrex::MultiFab*, 3 > const Efield, + std::array< amrex::MultiFab*, 3 > const psiB, + MultiSigmaBox const& sigba, + amrex::Real const dt, DtType const dt_type ); + + void EvolveECPML ( std::array< amrex::MultiFab*, 3 > Efield, + std::array< amrex::MultiFab*, 3 > const Bfield, + std::array< amrex::MultiFab*, 3 > const Jfield, + std::array< amrex::MultiFab*, 3 > const psiE, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles ); + + /** BorHfield is B (H = B/mu computed with the local mu) or, with LLG, H */ + void MacroscopicEvolveECPML ( std::array< amrex::MultiFab*, 3 > Efield, + std::array< amrex::MultiFab*, 3 > const BorHfield, + std::array< amrex::MultiFab*, 3 > const Jfield, + std::array< amrex::MultiFab*, 3 > const psiE, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const eps_mf, + amrex::MultiFab* const mu_mf, + amrex::MultiFab* const sigma_mf ); + +#ifdef WARPX_MAG_LLG + void EvolveHCPML ( std::array< amrex::MultiFab*, 3 > Hfield, + std::array< amrex::MultiFab*, 3 > const Efield, + std::array< amrex::MultiFab*, 3 > const psiH, + MultiSigmaBox const& sigba, + amrex::Real const dt, DtType const dt_type ); +#endif + #ifndef WARPX_DIM_RZ #ifdef WARPX_MAG_LLG void EvolveHPML ( std::array< amrex::MultiFab*, 3 > Hfield, @@ -398,6 +439,19 @@ class FiniteDifferenceSolver amrex::MultiFab* const sigma_mf); + template< typename T_MacroAlgo > + void MacroscopicEvolveECPMLCartesian ( + std::array< amrex::MultiFab*, 3 > Efield, + std::array< amrex::MultiFab*, 3 > const BorHfield, + std::array< amrex::MultiFab*, 3 > const Jfield, + std::array< amrex::MultiFab*, 3 > const psiE, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const eps_mf, + amrex::MultiFab* const mu_mf, + amrex::MultiFab* const sigma_mf); + #ifdef WARPX_MAG_LLG template< typename T_Algo > void EvolveHPMLCartesian ( diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package index 0f96e6bb5..3d374cfe9 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package +++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package @@ -14,6 +14,7 @@ CEXE_sources += EvolveHPML.cpp #endif CEXE_sources += EvolveBPML.cpp +CEXE_sources += EvolveCPML.cpp CEXE_sources += EvolveEPML.cpp CEXE_sources += EvolveFPML.cpp CEXE_sources += MacroscopicEvolveEPML.cpp diff --git a/Source/FieldSolver/WarpXExternalEMFields.cpp b/Source/FieldSolver/WarpXExternalEMFields.cpp index 4b7d19e38..7ce721264 100644 --- a/Source/FieldSolver/WarpXExternalEMFields.cpp +++ b/Source/FieldSolver/WarpXExternalEMFields.cpp @@ -51,6 +51,14 @@ WarpX::ApplyExternalFieldExcitationOnGrid (int const externalfieldtype, DtType a Eyfield_flag_parser->compile<3>(), Ezfield_flag_parser->compile<3>(), lev, a_dt_type ); + if (pml[lev]->isConvolutional()) { + // With the CPML the total field is in component 0 and component 1 + // must stay zero (the excitation above was applied to every component) + for (int idir = 0; idir < 3; ++idir) { + amrex::MultiFab* mf = pml[lev]->GetE_fp(idir); + mf->setVal(0._rt, 1, mf->nComp()-1, mf->nGrowVect()); + } + } } } if (externalfieldtype == ExternalFieldType::AllExternal || externalfieldtype == ExternalFieldType::BfieldExternal) { diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index aa7e6a768..ea58340eb 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -865,7 +865,17 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt, DtType a_dt_typ // Evolve B field in PML cells if (do_pml && pml[lev]->ok()) { - if (patch_type == PatchType::fine) { + if (pml[lev]->isConvolutional()) { + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->EvolveBCPML( + pml[lev]->GetB_fp(), pml[lev]->GetE_fp(), pml[lev]->GetPsiB_fp(), + pml[lev]->GetMultiSigmaBox_fp(), a_dt, a_dt_type); + } else { + m_fdtd_solver_cp[lev]->EvolveBCPML( + pml[lev]->GetB_cp(), pml[lev]->GetE_cp(), pml[lev]->GetPsiB_cp(), + pml[lev]->GetMultiSigmaBox_cp(), a_dt, a_dt_type); + } + } else if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->EvolveBPML( pml[lev]->GetB_fp(), pml[lev]->GetE_fp(), a_dt, WarpX::do_dive_cleaning); } else { @@ -915,7 +925,29 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) } // Evolve E field in PML cells - if (do_pml && pml[lev]->ok()) { + if (do_pml && pml[lev]->ok() && pml[lev]->isConvolutional()) { + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->EvolveECPML( + pml[lev]->GetE_fp(), +#ifdef WARPX_MAG_LLG + pml[lev]->GetH_fp(), +#else + pml[lev]->GetB_fp(), +#endif + pml[lev]->Getj_fp(), pml[lev]->GetPsiE_fp(), + pml[lev]->GetMultiSigmaBox_fp(), a_dt, pml_has_particles ); + } else { + m_fdtd_solver_cp[lev]->EvolveECPML( + pml[lev]->GetE_cp(), +#ifdef WARPX_MAG_LLG + pml[lev]->GetH_cp(), +#else + pml[lev]->GetB_cp(), +#endif + pml[lev]->Getj_cp(), pml[lev]->GetPsiE_cp(), + pml[lev]->GetMultiSigmaBox_cp(), a_dt, pml_has_particles ); + } + } else if (do_pml && pml[lev]->ok()) { if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->EvolveEPML( pml[lev]->GetE_fp(), @@ -1094,7 +1126,41 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { m_macroscopic_properties); // Evolve E field in PML cells if (do_pml && pml[lev]->ok()) { - if (patch_type == PatchType::fine) { + if (pml[lev]->isConvolutional()) { + // CPML: same material update (alpha/beta with local eps, sigma; H = B/mu + // with local mu) as in the interior, applied to the stretched curl + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->MacroscopicEvolveECPML( + pml[lev]->GetE_fp(), +#ifndef WARPX_MAG_LLG + pml[lev]->GetB_fp(), +#else + pml[lev]->GetH_fp(), +#endif + pml[lev]->Getj_fp(), pml[lev]->GetPsiE_fp(), + pml[lev]->GetMultiSigmaBox_fp(), + a_dt, pml_has_particles, + m_macroscopic_properties, + pml[lev]->Geteps_fp(), + pml[lev]->Getmu_fp(), + pml[lev]->Getsigma_fp() ); + } else { + m_fdtd_solver_cp[lev]->MacroscopicEvolveECPML( + pml[lev]->GetE_cp(), +#ifndef WARPX_MAG_LLG + pml[lev]->GetB_cp(), +#else + pml[lev]->GetH_cp(), +#endif + pml[lev]->Getj_cp(), pml[lev]->GetPsiE_cp(), + pml[lev]->GetMultiSigmaBox_cp(), + a_dt, pml_has_particles, + m_macroscopic_properties, + pml[lev]->Geteps_cp(), + pml[lev]->Getmu_cp(), + pml[lev]->Getsigma_cp() ); + } + } else if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->MacroscopicEvolveEPML( pml[lev]->GetE_fp(), #ifndef WARPX_MAG_LLG @@ -1135,9 +1201,24 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { amrex::MultiFab* PECz_PEC = pml[lev]->GetPEC_fp(2); // 0 ghost cells or nge? + // Split-field PML: both split components are zeroed on the PEC edges. + // CPML: component 0 holds the total field (component 1 is zero), so the + // same multiply enforces E_tan = 0 on the PEC mask extended into the PML. MultiFab::Multiply(*Ex_PEC, *PECx_PEC, 0, 0, 2, 0); MultiFab::Multiply(*Ey_PEC, *PECy_PEC, 0, 0, 2, 0); MultiFab::Multiply(*Ez_PEC, *PECz_PEC, 0, 0, 2, 0); + + if (pml[lev]->isConvolutional()) { + // Also reset the CPML memory variables on PEC edges, so that no + // stale convolution history is carried by masked E components + const auto psiE = pml[lev]->GetPsiE_fp(); + MultiFab::Multiply(*psiE[0], *PECx_PEC, 0, 0, 1, 0); + MultiFab::Multiply(*psiE[0], *PECx_PEC, 0, 1, 1, 0); + MultiFab::Multiply(*psiE[1], *PECy_PEC, 0, 0, 1, 0); + MultiFab::Multiply(*psiE[1], *PECy_PEC, 0, 1, 1, 0); + MultiFab::Multiply(*psiE[2], *PECz_PEC, 0, 0, 1, 0); + MultiFab::Multiply(*psiE[2], *PECz_PEC, 0, 1, 1, 0); + } } } @@ -1149,25 +1230,25 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { #ifdef WARPX_MAG_LLG // define WarpX::MacroscopicEvolveHM void -WarpX::MacroscopicEvolveHM (amrex::Real a_dt) +WarpX::MacroscopicEvolveHM (amrex::Real a_dt, DtType a_dt_type) { for (int lev = 0; lev <= finest_level; ++lev ) { - MacroscopicEvolveHM(lev, a_dt); + MacroscopicEvolveHM(lev, a_dt, a_dt_type); } } void -WarpX::MacroscopicEvolveHM (int lev, amrex::Real a_dt) { +WarpX::MacroscopicEvolveHM (int lev, amrex::Real a_dt, DtType a_dt_type) { WARPX_PROFILE("WarpX::MacroscopicEvolveHM()"); - MacroscopicEvolveHM(lev, PatchType::fine, a_dt); + MacroscopicEvolveHM(lev, PatchType::fine, a_dt, a_dt_type); if (lev > 0) { amrex::Abort("Macroscopic EvolveHM is not implemented for lev>0, yet."); } } void -WarpX::MacroscopicEvolveHM (int lev, PatchType patch_type, amrex::Real a_dt) { +WarpX::MacroscopicEvolveHM (int lev, PatchType patch_type, amrex::Real a_dt, DtType a_dt_type) { // Evolve H field in regular cells if (patch_type == PatchType::fine) { @@ -1179,7 +1260,17 @@ WarpX::MacroscopicEvolveHM (int lev, PatchType patch_type, amrex::Real a_dt) { } // Evolve H field in PML cells - if (do_pml && pml[lev]->ok()) { + if (do_pml && pml[lev]->ok() && pml[lev]->isConvolutional()) { + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->EvolveHCPML( + pml[lev]->GetH_fp(), pml[lev]->GetE_fp(), pml[lev]->GetPsiB_fp(), + pml[lev]->GetMultiSigmaBox_fp(), a_dt, a_dt_type); + } else { + m_fdtd_solver_cp[lev]->EvolveHCPML( + pml[lev]->GetH_cp(), pml[lev]->GetE_cp(), pml[lev]->GetPsiB_cp(), + pml[lev]->GetMultiSigmaBox_cp(), a_dt, a_dt_type); + } + } else if (do_pml && pml[lev]->ok()) { if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->EvolveHPML( pml[lev]->GetH_fp(), pml[lev]->GetE_fp(), a_dt, WarpX::do_dive_cleaning); @@ -1192,25 +1283,25 @@ WarpX::MacroscopicEvolveHM (int lev, PatchType patch_type, amrex::Real a_dt) { // define WarpX::MacroscopicEvolveHM_2nd void -WarpX::MacroscopicEvolveHM_2nd (amrex::Real a_dt) +WarpX::MacroscopicEvolveHM_2nd (amrex::Real a_dt, DtType a_dt_type) { for (int lev = 0; lev <= finest_level; ++lev ) { - MacroscopicEvolveHM_2nd(lev, a_dt); + MacroscopicEvolveHM_2nd(lev, a_dt, a_dt_type); } } void -WarpX::MacroscopicEvolveHM_2nd (int lev, amrex::Real a_dt) { +WarpX::MacroscopicEvolveHM_2nd (int lev, amrex::Real a_dt, DtType a_dt_type) { WARPX_PROFILE("WarpX::MacroscopicEvolveHM_2nd()"); - MacroscopicEvolveHM_2nd(lev, PatchType::fine, a_dt); + MacroscopicEvolveHM_2nd(lev, PatchType::fine, a_dt, a_dt_type); if (lev > 0) { amrex::Abort("Macroscopic EvolveHM_2nd is not implemented for lev>0, yet."); } } void -WarpX::MacroscopicEvolveHM_2nd (int lev, PatchType patch_type, amrex::Real a_dt) { +WarpX::MacroscopicEvolveHM_2nd (int lev, PatchType patch_type, amrex::Real a_dt, DtType a_dt_type) { // Evolve H field in regular cells if (patch_type == PatchType::fine) { @@ -1222,7 +1313,17 @@ WarpX::MacroscopicEvolveHM_2nd (int lev, PatchType patch_type, amrex::Real a_dt) } // Evolve H field in PML cells - if (do_pml && pml[lev]->ok()) { + if (do_pml && pml[lev]->ok() && pml[lev]->isConvolutional()) { + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->EvolveHCPML( + pml[lev]->GetH_fp(), pml[lev]->GetE_fp(), pml[lev]->GetPsiB_fp(), + pml[lev]->GetMultiSigmaBox_fp(), a_dt, a_dt_type); + } else { + m_fdtd_solver_cp[lev]->EvolveHCPML( + pml[lev]->GetH_cp(), pml[lev]->GetE_cp(), pml[lev]->GetPsiB_cp(), + pml[lev]->GetMultiSigmaBox_cp(), a_dt, a_dt_type); + } + } else if (do_pml && pml[lev]->ok()) { if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->EvolveHPML( pml[lev]->GetH_fp(), pml[lev]->GetE_fp(), a_dt, WarpX::do_dive_cleaning ); diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index 0e5de537c..7b6b26e03 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -50,6 +50,19 @@ struct MacroscopicSolverAlgo { }; }; +/** + * \brief struct to select the PML formulation used with the FDTD solver + * Split: split-field (Berenger) PML (default, original implementation) + * Convolutional: convolutional PML (CPML, Roden & Gedney 2000) with the + * complex-frequency-shifted stretching s = kappa + sigma/(alpha + i omega) + */ +struct PMLType { + enum { + Split = 0, + Convolutional = 1 + }; +}; + struct GridType { enum { Collocated = 0, diff --git a/Source/WarpX.H b/Source/WarpX.H index 1cfe69217..d11799958 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -146,6 +146,17 @@ public: static std::string H_bias_ext_grid_s; #endif static int ApplyExcitationInPML; + + //! PML formulation: PMLType::Split (default) or PMLType::Convolutional (warpx.pml_type) + static int pml_type; + //! CPML: maximum of the real coordinate-stretching factor kappa (warpx.cpml_kappa_max, default 1) + static amrex::Real cpml_kappa_max; + //! CPML: maximum of the CFS parameter alpha, as a rate alpha_SI/eps0 in 1/s (warpx.cpml_alpha_max, default 0) + static amrex::Real cpml_alpha_max; + //! CPML: polynomial grading order of kappa (warpx.cpml_grading_order, default 2, same as sigma) + static amrex::Real cpml_grading_order; + //! CPML: polynomial grading order of alpha (warpx.cpml_alpha_grading_order, default 1) + static amrex::Real cpml_alpha_grading_order; //! Whether to apply the effect of an externally-defined electric field static bool add_external_E_field; //! Whether to apply the effect of an externally-defined magnetic field @@ -818,13 +829,13 @@ public: void MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real dt); #ifdef WARPX_MAG_LLG - void MacroscopicEvolveHM ( amrex::Real dt); - void MacroscopicEvolveHM (int lev, amrex::Real dt); - void MacroscopicEvolveHM (int lev, PatchType patch_type, amrex::Real dt); + void MacroscopicEvolveHM ( amrex::Real dt, DtType dt_type = DtType::Full); + void MacroscopicEvolveHM (int lev, amrex::Real dt, DtType dt_type = DtType::Full); + void MacroscopicEvolveHM (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type = DtType::Full); - void MacroscopicEvolveHM_2nd ( amrex::Real dt); - void MacroscopicEvolveHM_2nd (int lev, amrex::Real dt); - void MacroscopicEvolveHM_2nd (int lev, PatchType patch_type, amrex::Real dt); + void MacroscopicEvolveHM_2nd ( amrex::Real dt, DtType dt_type = DtType::Full); + void MacroscopicEvolveHM_2nd (int lev, amrex::Real dt, DtType dt_type = DtType::Full); + void MacroscopicEvolveHM_2nd (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type = DtType::Full); #endif /** apply QED correction on electric field diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 11a95f6c5..95a8e4d53 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -96,6 +96,11 @@ std::string WarpX::authors = ""; std::string WarpX::B_ext_grid_s = "default"; std::string WarpX::E_ext_grid_s = "default"; int WarpX::ApplyExcitationInPML = 0; // default is 0 +int WarpX::pml_type = PMLType::Split; +amrex::Real WarpX::cpml_kappa_max = 1.0_rt; +amrex::Real WarpX::cpml_alpha_max = 0.0_rt; +amrex::Real WarpX::cpml_grading_order = 2.0_rt; +amrex::Real WarpX::cpml_alpha_grading_order = 1.0_rt; // default type is constant std::string WarpX::B_excitation_grid_s = "default"; @@ -978,6 +983,32 @@ 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); + + // PML formulation: "split" (default) or "convolutional" (CPML) + { + std::string pml_type_str = "split"; + pp_warpx.query("pml_type", pml_type_str); + std::transform(pml_type_str.begin(), pml_type_str.end(), pml_type_str.begin(), ::tolower); + if (pml_type_str == "split") { + pml_type = PMLType::Split; + } else if (pml_type_str == "convolutional" || pml_type_str == "cpml") { + pml_type = PMLType::Convolutional; + } else { + WARPX_ABORT_WITH_MESSAGE("warpx.pml_type must be 'split' or 'convolutional'"); + } + utils::parser::queryWithParser(pp_warpx, "cpml_kappa_max", cpml_kappa_max); + utils::parser::queryWithParser(pp_warpx, "cpml_alpha_max", cpml_alpha_max); + utils::parser::queryWithParser(pp_warpx, "cpml_grading_order", cpml_grading_order); + utils::parser::queryWithParser(pp_warpx, "cpml_alpha_grading_order", cpml_alpha_grading_order); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(cpml_kappa_max >= 1._rt, + "warpx.cpml_kappa_max must be >= 1"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(cpml_alpha_max >= 0._rt, + "warpx.cpml_alpha_max must be >= 0"); + if (pml_type == PMLType::Convolutional) { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(do_pml_j_damping == 0, + "warpx.do_pml_j_damping is not used/supported with warpx.pml_type = convolutional"); + } + } // 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); From 3040831d971b39c597a409f4f293ad74d3738b27 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Fri, 30 Oct 2026 14:50:37 -0700 Subject: [PATCH 2/3] add files --- Source/BoundaryConditions/CPML_K.H | 138 ++++++ .../FiniteDifferenceSolver/EvolveCPML.cpp | 426 ++++++++++++++++++ 2 files changed, 564 insertions(+) create mode 100644 Source/BoundaryConditions/CPML_K.H create mode 100644 Source/FieldSolver/FiniteDifferenceSolver/EvolveCPML.cpp diff --git a/Source/BoundaryConditions/CPML_K.H b/Source/BoundaryConditions/CPML_K.H new file mode 100644 index 000000000..7919b3331 --- /dev/null +++ b/Source/BoundaryConditions/CPML_K.H @@ -0,0 +1,138 @@ +/* Convolutional PML (CPML) helpers. + * + * The CPML uses the complex-frequency-shifted coordinate stretching + * + * s_d = kappa_d + sigma_d / (alpha_d + i omega) (d = x, y, z) + * + * (sigma_d and alpha_d are rates, i.e. sigma_SI/eps0 and alpha_SI/eps0), which turns + * each spatial derivative into + * + * d/dd -> (1/kappa_d) d/dd + psi_d + * + * where psi_d is the recursive convolution (Roden & Gedney, 2000) + * + * psi_d^{new} = b_d psi_d^{old} + a_d (dF/dd), + * b_d = exp(-(sigma_d/kappa_d + alpha_d) dt), + * a_d = sigma_d / (kappa_d (sigma_d + kappa_d alpha_d)) (b_d - 1). + * + * Because the stretching only acts on the spatial derivatives, the material + * response (eps, mu, sigma of the medium) is untouched: the E update keeps the + * same alpha/beta (Lax-Wendroff or backward-Euler) coefficients as in the + * interior, evaluated with the local eps(x) and sigma(x). This is what makes the + * CPML reflectionless for arbitrary (and spatially varying) materials continued + * into the PML, in contrast to formulations whose damping must be matched to eps. + */ +#ifndef WARPX_CPML_K_H_ +#define WARPX_CPML_K_H_ + +#include "BoundaryConditions/PML.H" + +#include +#include +#include +#include +#include + +#include + +/** Device-copyable view of the CPML 1D profiles of one SigmaBox, indexed by the + * physical direction (0=x, 1=y, 2=z). Directions that are not resolved by the + * simulation (y in 2D XZ) have null pointers and act as "no stretching". */ +struct CPMLCoefs +{ + const amrex::Real* sig [3] = {nullptr, nullptr, nullptr}; // nodal positions + const amrex::Real* sigs[3] = {nullptr, nullptr, nullptr}; // cell-centered positions + const amrex::Real* kap [3] = {nullptr, nullptr, nullptr}; + const amrex::Real* kaps[3] = {nullptr, nullptr, nullptr}; + const amrex::Real* alp [3] = {nullptr, nullptr, nullptr}; + const amrex::Real* alps[3] = {nullptr, nullptr, nullptr}; + int lo[3] = {0, 0, 0}; + + CPMLCoefs () = default; + + explicit CPMLCoefs (SigmaBox const& sb) + { +#if defined(WARPX_DIM_3D) + constexpr int dirs[AMREX_SPACEDIM] = {0, 1, 2}; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + constexpr int dirs[AMREX_SPACEDIM] = {0, 2}; +#else + constexpr int dirs[AMREX_SPACEDIM] = {2}; +#endif + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + const int d = dirs[idim]; + sig [d] = sb.sigma[idim].data(); + sigs[d] = sb.sigma_star[idim].data(); + kap [d] = sb.kappa[idim].data(); + kaps[d] = sb.kappa_star[idim].data(); + alp [d] = sb.alpha[idim].data(); + alps[d] = sb.alpha_star[idim].data(); + lo [d] = sb.sigma[idim].lo(); + } + } +}; + +/** Staggering of a field component, indexed by physical direction (1 = nodal) */ +inline amrex::GpuArray cpml_stag (amrex::IndexType const& ixt) +{ + amrex::IntVect const iv = ixt.toIntVect(); +#if defined(WARPX_DIM_3D) + return {iv[0], iv[1], iv[2]}; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + return {iv[0], 0, iv[1]}; +#else + return {0, 0, iv[0]}; +#endif +} + +/** Grid index along physical direction dir */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +int cpml_index (int dir, int i, int j, int k) noexcept +{ +#if defined(WARPX_DIM_3D) + return (dir == 0) ? i : ((dir == 1) ? j : k); +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + amrex::ignore_unused(k); + return (dir == 0) ? i : j; +#else + amrex::ignore_unused(dir, j, k); + return i; +#endif +} + +/** Coefficients (b, a, 1/kappa) of the CPML recursive convolution along direction + * dir, at grid index ii, for a field with nodal (1) or cell-centered (0) staggering + * along dir, for a time step dt. */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void cpml_coef (CPMLCoefs const& c, int dir, int nodal, int ii, amrex::Real dt, + amrex::Real& b, amrex::Real& a, amrex::Real& inv_kappa) noexcept +{ + using namespace amrex::literals; + if (c.sig[dir] == nullptr) { + b = 1._rt; a = 0._rt; inv_kappa = 1._rt; + return; + } + const int n = ii - c.lo[dir]; + const amrex::Real s = nodal ? c.sig[dir][n] : c.sigs[dir][n]; + const amrex::Real kp = nodal ? c.kap[dir][n] : c.kaps[dir][n]; + const amrex::Real al = nodal ? c.alp[dir][n] : c.alps[dir][n]; + inv_kappa = 1._rt/kp; + b = std::exp(-(s*inv_kappa + al)*dt); + const amrex::Real denom = kp*(s + kp*al); + a = (denom > 0._rt) ? s/denom*(b - 1._rt) : 0._rt; +} + +/** Stretched derivative: if update_psi, first advances the memory variable psi + * (in place) over dt_psi; returns (1/kappa) dF + psi */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real cpml_stretched_derivative (CPMLCoefs const& c, int dir, int nodal, + int i, int j, int k, amrex::Real dt_psi, bool update_psi, + amrex::Real dF, amrex::Real& psi) noexcept +{ + amrex::Real b, a, ik; + cpml_coef(c, dir, nodal, cpml_index(dir, i, j, k), dt_psi, b, a, ik); + if (update_psi) { psi = b*psi + a*dF; } + return ik*dF + psi; +} + +#endif // WARPX_CPML_K_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveCPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveCPML.cpp new file mode 100644 index 000000000..62b1c789c --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveCPML.cpp @@ -0,0 +1,426 @@ +/* Convolutional PML (CPML) field updates for the Cartesian Yee FDTD solver. + * + * These are used when warpx.pml_type = convolutional. The PML field MultiFabs keep + * the same layout as for the split-field PML (2 or 3 components), so that the + * existing PML <-> domain exchange, FillBoundary, checkpointing, PEC handling, etc. + * can be reused unchanged: with the CPML the *total* field is stored in component 0 + * and component 1 is identically zero (PML::Exchange sums the components, and sets + * component 1 to zero in the guard cells it fills from the regular grid). + * + * See BoundaryConditions/CPML_K.H for the formulation. + */ +#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H" + +#include "BoundaryConditions/CPML_K.H" +#include "BoundaryConditions/PML.H" +#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H" +#ifndef WARPX_DIM_RZ +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" +# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H" +#endif +#include "Utils/TextMsg.H" +#include "Utils/WarpXAlgorithmSelection.H" +#include "Utils/WarpXConst.H" +#include "WarpX.H" + +#include + +#include +#include +#include +#include +#include + +#include + +using namespace amrex; + +namespace { + +#ifndef WARPX_DIM_RZ + void check_cpml_algo (int fdtd_algo, short grid_type) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + fdtd_algo == ElectromagneticSolverAlgo::Yee && grid_type != GridType::Collocated, + "Convolutional PML is only implemented for the staggered Yee solver"); + } + + /** B (or H) update in the CPML: F -= fac * dt * curl_s(E). + * fac = 1 for B, fac = 1/mu0 for H (LLG). + * + * B is advanced in two half steps (FirstHalf with E^n, SecondHalf with E^{n+1}), + * which together form the leapfrog step B^{n-1/2} -> B^{n+1/2} driven by E^n. + * The memory variable is therefore updated once per leapfrog step, at the start + * of the SecondHalf (psi^{n+1} = b psi^n + a dE^{n+1}, with the coefficients of + * the full step 2*dt), and that same psi is used again by the following FirstHalf + * (E, hence dE, is unchanged in between). This makes the CPML recursion + * consistent with leapfrog; updating psi at every half step with half-step + * coefficients is noticeably more reflective. + * A DtType::Full push updates psi with dt and applies it over dt. */ + void EvolveBorHCPMLYee ( + std::array< amrex::MultiFab*, 3 > Bfield, + std::array< amrex::MultiFab*, 3 > const Efield, + std::array< amrex::MultiFab*, 3 > const psiB, + MultiSigmaBox const& sigba, + amrex::Real const a_dt, DtType const dt_type, amrex::Real const fac, + amrex::Real const * const coefs_x_dev, int const n_coefs_x, + amrex::Real const * const coefs_y_dev, int const n_coefs_y, + amrex::Real const * const coefs_z_dev, int const n_coefs_z) + { + using Yee = CartesianYeeAlgorithm; + + // time step applied to B, and time step of the psi recursion + amrex::Real const dt = a_dt; + bool const update_psi = (dt_type != DtType::FirstHalf); + amrex::Real const dt_psi = (dt_type == DtType::SecondHalf) ? 2._rt*a_dt : a_dt; + + auto const Bx_stag = cpml_stag(Bfield[0]->ixType()); + auto const By_stag = cpml_stag(Bfield[1]->ixType()); + auto const Bz_stag = cpml_stag(Bfield[2]->ixType()); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Bfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + Array4 const& Bx = Bfield[0]->array(mfi); + Array4 const& By = Bfield[1]->array(mfi); + Array4 const& Bz = Bfield[2]->array(mfi); + Array4 const& Ex = Efield[0]->array(mfi); + Array4 const& Ey = Efield[1]->array(mfi); + Array4 const& Ez = Efield[2]->array(mfi); + Array4 const& Px = psiB[0]->array(mfi); + Array4 const& Py = psiB[1]->array(mfi); + Array4 const& Pz = psiB[2]->array(mfi); + + CPMLCoefs const cc(sigba[mfi]); + + Real const * const AMREX_RESTRICT coefs_x = coefs_x_dev; + Real const * const AMREX_RESTRICT coefs_y = coefs_y_dev; + Real const * const AMREX_RESTRICT coefs_z = coefs_z_dev; + + Box const& tbx = mfi.tilebox(Bfield[0]->ixType().toIntVect()); + Box const& tby = mfi.tilebox(Bfield[1]->ixType().toIntVect()); + Box const& tbz = mfi.tilebox(Bfield[2]->ixType().toIntVect()); + + amrex::ParallelFor(tbx, tby, tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + // curl_x E = dEz/dy - dEy/dz ; psi comps: xy=0 (d/dy), xz=1 (d/dz) + Real const dEz_dy = Yee::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, 0); + Real const dEy_dz = Yee::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, 0); + Real const sy = cpml_stretched_derivative(cc, 1, Bx_stag[1], i, j, k, dt_psi, update_psi, dEz_dy, Px(i,j,k,0)); + Real const sz = cpml_stretched_derivative(cc, 2, Bx_stag[2], i, j, k, dt_psi, update_psi, dEy_dz, Px(i,j,k,1)); + Bx(i, j, k, 0) -= fac * dt * (sy - sz); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + // curl_y E = dEx/dz - dEz/dx ; psi comps: yz=0 (d/dz), yx=1 (d/dx) + Real const dEx_dz = Yee::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, 0); + Real const dEz_dx = Yee::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, 0); + Real const sz = cpml_stretched_derivative(cc, 2, By_stag[2], i, j, k, dt_psi, update_psi, dEx_dz, Py(i,j,k,0)); + Real const sx = cpml_stretched_derivative(cc, 0, By_stag[0], i, j, k, dt_psi, update_psi, dEz_dx, Py(i,j,k,1)); + By(i, j, k, 0) -= fac * dt * (sz - sx); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + // curl_z E = dEy/dx - dEx/dy ; psi comps: zx=0 (d/dx), zy=1 (d/dy) + Real const dEy_dx = Yee::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, 0); + Real const dEx_dy = Yee::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, 0); + Real const sx = cpml_stretched_derivative(cc, 0, Bz_stag[0], i, j, k, dt_psi, update_psi, dEy_dx, Pz(i,j,k,0)); + Real const sy = cpml_stretched_derivative(cc, 1, Bz_stag[1], i, j, k, dt_psi, update_psi, dEx_dy, Pz(i,j,k,1)); + Bz(i, j, k, 0) -= fac * dt * (sx - sy); + } + ); + } + } +#endif // WARPX_DIM_RZ +} + +void FiniteDifferenceSolver::EvolveBCPML ( + std::array< amrex::MultiFab*, 3 > Bfield, + std::array< amrex::MultiFab*, 3 > const Efield, + std::array< amrex::MultiFab*, 3 > const psiB, + MultiSigmaBox const& sigba, + amrex::Real const dt, DtType const dt_type) +{ +#ifdef WARPX_DIM_RZ + amrex::ignore_unused(Bfield, Efield, psiB, sigba, dt, dt_type); + WARPX_ABORT_WITH_MESSAGE("PML are not implemented in cylindrical geometry."); +#else + check_cpml_algo(m_fdtd_algo, m_grid_type); + EvolveBorHCPMLYee(Bfield, Efield, psiB, sigba, dt, dt_type, 1._rt, + m_stencil_coefs_x.dataPtr(), static_cast(m_stencil_coefs_x.size()), + m_stencil_coefs_y.dataPtr(), static_cast(m_stencil_coefs_y.size()), + m_stencil_coefs_z.dataPtr(), static_cast(m_stencil_coefs_z.size())); +#endif +} + +#ifdef WARPX_MAG_LLG +void FiniteDifferenceSolver::EvolveHCPML ( + std::array< amrex::MultiFab*, 3 > Hfield, + std::array< amrex::MultiFab*, 3 > const Efield, + std::array< amrex::MultiFab*, 3 > const psiH, + MultiSigmaBox const& sigba, + amrex::Real const dt, DtType const dt_type) +{ +#ifdef WARPX_DIM_RZ + amrex::ignore_unused(Hfield, Efield, psiH, sigba, dt, dt_type); + WARPX_ABORT_WITH_MESSAGE("PML are not implemented in cylindrical geometry."); +#else + check_cpml_algo(m_fdtd_algo, m_grid_type); + // As in EvolveHPML, the PML region is treated as vacuum for the H update (mu = mu0) + EvolveBorHCPMLYee(Hfield, Efield, psiH, sigba, dt, dt_type, 1._rt/PhysConst::mu0, + m_stencil_coefs_x.dataPtr(), static_cast(m_stencil_coefs_x.size()), + m_stencil_coefs_y.dataPtr(), static_cast(m_stencil_coefs_y.size()), + m_stencil_coefs_z.dataPtr(), static_cast(m_stencil_coefs_z.size())); +#endif +} +#endif + +void FiniteDifferenceSolver::EvolveECPML ( + std::array< amrex::MultiFab*, 3 > Efield, + std::array< amrex::MultiFab*, 3 > const Bfield, + std::array< amrex::MultiFab*, 3 > const Jfield, + std::array< amrex::MultiFab*, 3 > const psiE, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles) +{ +#ifdef WARPX_DIM_RZ + amrex::ignore_unused(Efield, Bfield, Jfield, psiE, sigba, dt, pml_has_particles); + WARPX_ABORT_WITH_MESSAGE("PML are not implemented in cylindrical geometry."); +#else + check_cpml_algo(m_fdtd_algo, m_grid_type); + using Yee = CartesianYeeAlgorithm; + + // dE/dt = c^2 curl_s(B) - mu0 c^2 J (with LLG, Bfield holds H: dE/dt = curl_s(H)/eps0 - J/eps0) + Real c2 = PhysConst::c * PhysConst::c; +#ifdef WARPX_MAG_LLG + c2 *= PhysConst::mu0; +#endif + Real const jfac = PhysConst::c * PhysConst::c * PhysConst::mu0; // 1/eps0 + + auto const Ex_stag = cpml_stag(Efield[0]->ixType()); + auto const Ey_stag = cpml_stag(Efield[1]->ixType()); + auto const Ez_stag = cpml_stag(Efield[2]->ixType()); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Efield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + Array4 const& Ex = Efield[0]->array(mfi); + Array4 const& Ey = Efield[1]->array(mfi); + Array4 const& Ez = Efield[2]->array(mfi); + Array4 const& Bx = Bfield[0]->array(mfi); + Array4 const& By = Bfield[1]->array(mfi); + Array4 const& Bz = Bfield[2]->array(mfi); + Array4 const& Px = psiE[0]->array(mfi); + Array4 const& Py = psiE[1]->array(mfi); + Array4 const& Pz = psiE[2]->array(mfi); + + CPMLCoefs const cc(sigba[mfi]); + + Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + int const n_coefs_x = static_cast(m_stencil_coefs_x.size()); + Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = static_cast(m_stencil_coefs_y.size()); + Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = static_cast(m_stencil_coefs_z.size()); + + Box const& tex = mfi.tilebox(Efield[0]->ixType().toIntVect()); + Box const& tey = mfi.tilebox(Efield[1]->ixType().toIntVect()); + Box const& tez = mfi.tilebox(Efield[2]->ixType().toIntVect()); + + amrex::ParallelFor(tex, tey, tez, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Real const dBz_dy = Yee::DownwardDy(Bz, coefs_y, n_coefs_y, i, j, k, 0); + Real const dBy_dz = Yee::DownwardDz(By, coefs_z, n_coefs_z, i, j, k, 0); + Real const sy = cpml_stretched_derivative(cc, 1, Ex_stag[1], i, j, k, dt, true, dBz_dy, Px(i,j,k,0)); + Real const sz = cpml_stretched_derivative(cc, 2, Ex_stag[2], i, j, k, dt, true, dBy_dz, Px(i,j,k,1)); + Ex(i, j, k, 0) += c2 * dt * (sy - sz); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Real const dBx_dz = Yee::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k, 0); + Real const dBz_dx = Yee::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k, 0); + Real const sz = cpml_stretched_derivative(cc, 2, Ey_stag[2], i, j, k, dt, true, dBx_dz, Py(i,j,k,0)); + Real const sx = cpml_stretched_derivative(cc, 0, Ey_stag[0], i, j, k, dt, true, dBz_dx, Py(i,j,k,1)); + Ey(i, j, k, 0) += c2 * dt * (sz - sx); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Real const dBy_dx = Yee::DownwardDx(By, coefs_x, n_coefs_x, i, j, k, 0); + Real const dBx_dy = Yee::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k, 0); + Real const sx = cpml_stretched_derivative(cc, 0, Ez_stag[0], i, j, k, dt, true, dBy_dx, Pz(i,j,k,0)); + Real const sy = cpml_stretched_derivative(cc, 1, Ez_stag[1], i, j, k, dt, true, dBx_dy, Pz(i,j,k,1)); + Ez(i, j, k, 0) += c2 * dt * (sx - sy); + } + ); + + if (pml_has_particles) { + Array4 const& Jx = Jfield[0]->array(mfi); + Array4 const& Jy = Jfield[1]->array(mfi); + Array4 const& Jz = Jfield[2]->array(mfi); + amrex::ParallelFor(tex, tey, tez, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ Ex(i,j,k,0) -= jfac * dt * Jx(i,j,k); }, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ Ey(i,j,k,0) -= jfac * dt * Jy(i,j,k); }, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ Ez(i,j,k,0) -= jfac * dt * Jz(i,j,k); } + ); + } + } +#endif +} + +void FiniteDifferenceSolver::MacroscopicEvolveECPML ( + std::array< amrex::MultiFab*, 3 > Efield, + std::array< amrex::MultiFab*, 3 > const BorHfield, + std::array< amrex::MultiFab*, 3 > const Jfield, + std::array< amrex::MultiFab*, 3 > const psiE, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const eps_mf, + amrex::MultiFab* const mu_mf, + amrex::MultiFab* const sigma_mf) +{ +#ifdef WARPX_DIM_RZ + amrex::ignore_unused(Efield, BorHfield, Jfield, psiE, sigba, dt, pml_has_particles, + macroscopic_properties, eps_mf, mu_mf, sigma_mf); + WARPX_ABORT_WITH_MESSAGE("PML are not implemented in cylindrical geometry."); +#else + check_cpml_algo(m_fdtd_algo, m_grid_type); + if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) { + MacroscopicEvolveECPMLCartesian( + Efield, BorHfield, Jfield, psiE, sigba, dt, pml_has_particles, + macroscopic_properties, eps_mf, mu_mf, sigma_mf); + } else if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) { + MacroscopicEvolveECPMLCartesian( + Efield, BorHfield, Jfield, psiE, sigba, dt, pml_has_particles, + macroscopic_properties, eps_mf, mu_mf, sigma_mf); + } else { + WARPX_ABORT_WITH_MESSAGE("MacroscopicEvolveECPML: unknown macroscopic algorithm"); + } +#endif +} + +#ifndef WARPX_DIM_RZ +template +void FiniteDifferenceSolver::MacroscopicEvolveECPMLCartesian ( + std::array< amrex::MultiFab*, 3 > Efield, + std::array< amrex::MultiFab*, 3 > const BorHfield, + std::array< amrex::MultiFab*, 3 > const Jfield, + std::array< amrex::MultiFab*, 3 > const psiE, + MultiSigmaBox const& sigba, + amrex::Real const dt, bool pml_has_particles, + std::unique_ptr const& macroscopic_properties, + amrex::MultiFab* const eps_mf, + amrex::MultiFab* const mu_mf, + amrex::MultiFab* const sigma_mf) +{ + using Yee = CartesianYeeAlgorithm; +#ifdef WARPX_MAG_LLG + amrex::ignore_unused(mu_mf); +#endif + + // Index types used to interpolate the (cell-centered) material properties to + // the Ex, Ey, Ez locations, exactly as in the interior and split-field PML updates + amrex::GpuArray const& sigma_stag = macroscopic_properties->sigma_IndexType; + amrex::GpuArray const& epsilon_stag = macroscopic_properties->epsilon_IndexType; + amrex::GpuArray const& macro_cr = macroscopic_properties->macro_cr_ratio; + amrex::GpuArray const& Ex_mstag = macroscopic_properties->Ex_IndexType; + amrex::GpuArray const& Ey_mstag = macroscopic_properties->Ey_IndexType; + amrex::GpuArray const& Ez_mstag = macroscopic_properties->Ez_IndexType; + + auto const Ex_stag = cpml_stag(Efield[0]->ixType()); + auto const Ey_stag = cpml_stag(Efield[1]->ixType()); + auto const Ez_stag = cpml_stag(Efield[2]->ixType()); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Efield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + Array4 const& Ex = Efield[0]->array(mfi); + Array4 const& Ey = Efield[1]->array(mfi); + Array4 const& Ez = Efield[2]->array(mfi); + Array4 const& Px = psiE[0]->array(mfi); + Array4 const& Py = psiE[1]->array(mfi); + Array4 const& Pz = psiE[2]->array(mfi); + + amrex::Array4 const& sigma_arr = sigma_mf->array(mfi); + amrex::Array4 const& eps_arr = eps_mf->array(mfi); + +#ifndef WARPX_MAG_LLG + // H = B/mu, with mu the local permeability in the PML + amrex::Array4 const& mu_arr = mu_mf->array(mfi); + FieldAccessorMacroscopic const Hx(BorHfield[0]->const_array(mfi), mu_arr); + FieldAccessorMacroscopic const Hy(BorHfield[1]->const_array(mfi), mu_arr); + FieldAccessorMacroscopic const Hz(BorHfield[2]->const_array(mfi), mu_arr); +#else + Array4 const& Hx = BorHfield[0]->array(mfi); + Array4 const& Hy = BorHfield[1]->array(mfi); + Array4 const& Hz = BorHfield[2]->array(mfi); +#endif + + CPMLCoefs const cc(sigba[mfi]); + + Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); + int const n_coefs_x = static_cast(m_stencil_coefs_x.size()); + Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = static_cast(m_stencil_coefs_y.size()); + Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = static_cast(m_stencil_coefs_z.size()); + + Box const& tex = mfi.tilebox(Efield[0]->ixType().toIntVect()); + Box const& tey = mfi.tilebox(Efield[1]->ixType().toIntVect()); + Box const& tez = mfi.tilebox(Efield[2]->ixType().toIntVect()); + + const int scomp = 0; + const bool has_j = pml_has_particles; + Array4 Jx, Jy, Jz; + if (has_j) { + Jx = Jfield[0]->const_array(mfi); + Jy = Jfield[1]->const_array(mfi); + Jz = Jfield[2]->const_array(mfi); + } + + amrex::ParallelFor(tex, tey, tez, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Real const sig = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, Ex_mstag, macro_cr, i, j, k, scomp); + Real const eps = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, Ex_mstag, macro_cr, i, j, k, scomp); + Real const alpha = T_MacroAlgo::alpha(sig, eps, dt); + Real const beta = T_MacroAlgo::beta (sig, eps, dt); + + Real const dHz_dy = Yee::DownwardDy(Hz, coefs_y, n_coefs_y, i, j, k, 0); + Real const dHy_dz = Yee::DownwardDz(Hy, coefs_z, n_coefs_z, i, j, k, 0); + Real const sy = cpml_stretched_derivative(cc, 1, Ex_stag[1], i, j, k, dt, true, dHz_dy, Px(i,j,k,0)); + Real const sz = cpml_stretched_derivative(cc, 2, Ex_stag[2], i, j, k, dt, true, dHy_dz, Px(i,j,k,1)); + Real const jx = has_j ? Jx(i,j,k) : 0._rt; + Ex(i, j, k, 0) = alpha * Ex(i, j, k, 0) + beta * (sy - sz - jx); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Real const sig = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, Ey_mstag, macro_cr, i, j, k, scomp); + Real const eps = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, Ey_mstag, macro_cr, i, j, k, scomp); + Real const alpha = T_MacroAlgo::alpha(sig, eps, dt); + Real const beta = T_MacroAlgo::beta (sig, eps, dt); + + Real const dHx_dz = Yee::DownwardDz(Hx, coefs_z, n_coefs_z, i, j, k, 0); + Real const dHz_dx = Yee::DownwardDx(Hz, coefs_x, n_coefs_x, i, j, k, 0); + Real const sz = cpml_stretched_derivative(cc, 2, Ey_stag[2], i, j, k, dt, true, dHx_dz, Py(i,j,k,0)); + Real const sx = cpml_stretched_derivative(cc, 0, Ey_stag[0], i, j, k, dt, true, dHz_dx, Py(i,j,k,1)); + Real const jy = has_j ? Jy(i,j,k) : 0._rt; + Ey(i, j, k, 0) = alpha * Ey(i, j, k, 0) + beta * (sz - sx - jy); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Real const sig = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, Ez_mstag, macro_cr, i, j, k, scomp); + Real const eps = ablastr::coarsen::sample::Interp(eps_arr, epsilon_stag, Ez_mstag, macro_cr, i, j, k, scomp); + Real const alpha = T_MacroAlgo::alpha(sig, eps, dt); + Real const beta = T_MacroAlgo::beta (sig, eps, dt); + + Real const dHy_dx = Yee::DownwardDx(Hy, coefs_x, n_coefs_x, i, j, k, 0); + Real const dHx_dy = Yee::DownwardDy(Hx, coefs_y, n_coefs_y, i, j, k, 0); + Real const sx = cpml_stretched_derivative(cc, 0, Ez_stag[0], i, j, k, dt, true, dHy_dx, Pz(i,j,k,0)); + Real const sy = cpml_stretched_derivative(cc, 1, Ez_stag[1], i, j, k, dt, true, dHx_dy, Pz(i,j,k,1)); + Real const jz = has_j ? Jz(i,j,k) : 0._rt; + Ez(i, j, k, 0) = alpha * Ez(i, j, k, 0) + beta * (sx - sy - jz); + } + ); + } +} +#endif // WARPX_DIM_RZ From 50b7c477390d307227ceb37db41533769e166232 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Wed, 4 Nov 2026 07:19:34 -0800 Subject: [PATCH 3/3] bugfix: main-grid E is now masked immediately after its update, before the exchange. The mask after the excitation is still applied, so masked edges stay exactly zero. new feature: algo.extend_PEC_mask_into_pml (default 1, same behaviour as before). Set it to 0 to have no PEC mask in the PML. --- Docs/source/usage/parameters.rst | 6 ++++++ Source/BoundaryConditions/PML.cpp | 14 ++++++++++---- Source/FieldSolver/WarpXPushFieldsEM.cpp | 18 +++++++++++++++++- Source/WarpX.H | 3 +++ Source/WarpX.cpp | 2 ++ 5 files changed, 38 insertions(+), 5 deletions(-) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 7427ad612..49112c0a0 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -404,6 +404,12 @@ Additional PML parameters typically lowers the reflection further; a too-large ``alpha_max`` degrades absorption of the low frequencies. +* ``algo.extend_PEC_mask_into_pml`` (`int`; default: 1) + Only used with ``algo.use_PEC_mask = 1``. If 1, the PEC mask built from the numpy (GDS) + layer is extended into the PML (the mask values on the domain edge are continued through + the PML), for both ``warpx.pml_type = split`` and ``convolutional``. If 0, no PEC mask is + applied in the PML. + * ``warpx.cpml_kappa_max`` (`float`; default: 1) CPML only. Maximum of the real stretching factor :math:`\kappa`, graded as :math:`\kappa(d) = 1 + (\kappa_{max}-1)\,d^m` with :math:`d \in [0,1]` the normalized diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index a76ed5d09..59e9a7026 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -827,11 +827,17 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri pml_PEC_fp[1]->setVal(1.); pml_PEC_fp[2]->setVal(1.); - macroscopic_properties->InitializePECFromSigma(pml_sigma_fp.get(), pml_PEC_fp[0].get(), pml_PEC_fp[1].get(), - macroscopic_properties->m_npy_k_index); - if (!macroscopic_properties->m_sigma_npy_filename2.empty()) { + // Extend the PEC layer(s) into the PML (algo.extend_PEC_mask_into_pml, default on). + // pml_sigma_fp was filled from the numpy mask with the edge values of the domain + // continued into the PML, so the PEC edges in the PML continue the interior ones. + // If disabled, the PEC mask stays 1 (no PEC) everywhere in the PML. + if (WarpX::extend_PEC_mask_into_pml) { macroscopic_properties->InitializePECFromSigma(pml_sigma_fp.get(), pml_PEC_fp[0].get(), pml_PEC_fp[1].get(), - macroscopic_properties->m_npy_k_index2); + macroscopic_properties->m_npy_k_index); + if (!macroscopic_properties->m_sigma_npy_filename2.empty()) { + macroscopic_properties->InitializePECFromSigma(pml_sigma_fp.get(), pml_PEC_fp[0].get(), pml_PEC_fp[1].get(), + macroscopic_properties->m_npy_k_index2); + } } // no need for sigma anymore diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index ea58340eb..223163d45 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -1124,6 +1124,22 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { #endif current_fp[lev], m_edge_lengths[lev], a_dt, m_macroscopic_properties); + + // Apply the PEC mask to the interior E right after the update, *before* the + // E guard cells and the PML are filled (FillBoundaryE -> PML::Exchange). + // Previously the mask was only applied after FillBoundaryE (in OneStep), so the + // guard cells of the interior grids, and the PML copy of the interface plane + // (which the PML receives from the interior), kept unmasked tangential E on the + // PEC edges until the next exchange. The following B update then saw + // inconsistent E on the two sides of the PML interface (and of interior box + // boundaries) along the PEC layer. The mask is applied again after the + // excitation in OneStep, so masked edges stay exactly zero. + if (use_PEC_mask) { + MultiFab::Multiply(*Efield_fp[lev][0], *PEC_fp[lev][0], 0, 0, 1, 0); + MultiFab::Multiply(*Efield_fp[lev][1], *PEC_fp[lev][1], 0, 0, 1, 0); + MultiFab::Multiply(*Efield_fp[lev][2], *PEC_fp[lev][2], 0, 0, 1, 0); + } + // Evolve E field in PML cells if (do_pml && pml[lev]->ok()) { if (pml[lev]->isConvolutional()) { @@ -1192,7 +1208,7 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { pml[lev]->Getsigma_cp() ); } - if (use_PEC_mask) { + if (use_PEC_mask && extend_PEC_mask_into_pml) { amrex::MultiFab* Ex_PEC = pml[lev]->GetE_fp(0); amrex::MultiFab* Ey_PEC = pml[lev]->GetE_fp(1); amrex::MultiFab* Ez_PEC = pml[lev]->GetE_fp(2); diff --git a/Source/WarpX.H b/Source/WarpX.H index d11799958..c5aa50c5e 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -330,6 +330,9 @@ public: static int yee_coupled_solver_algo; static int use_PEC_mask; + //! With use_PEC_mask: extend the PEC mask of the numpy/GDS layer into the PML + //! (algo.extend_PEC_mask_into_pml, default 1). If 0, the PML has no PEC mask. + static int extend_PEC_mask_into_pml; static int use_lumped_inductor; static int use_lumped_resistor; static int use_lumped_capacitor; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 95a8e4d53..c49e65a9f 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -229,6 +229,7 @@ amrex::Vector WarpX::particle_boundary_lo(AMREX_SPACEDIM,P amrex::Vector WarpX::particle_boundary_hi(AMREX_SPACEDIM,ParticleBoundaryType::Absorbing); int WarpX::yee_coupled_solver_algo; int WarpX::use_PEC_mask = 0; +int WarpX::extend_PEC_mask_into_pml = 1; int WarpX::use_lumped_inductor = 0; int WarpX::use_lumped_resistor = 0; int WarpX::use_lumped_capacitor = 0; @@ -1302,6 +1303,7 @@ WarpX::ReadParameters () yee_coupled_solver_algo = GetAlgorithmInteger(pp_algo, "yee_coupled_solver"); pp_algo.query("use_PEC_mask",use_PEC_mask); + pp_algo.query("extend_PEC_mask_into_pml",extend_PEC_mask_into_pml); pp_algo.query("use_lumped_inductor",use_lumped_inductor); // Load balancing parameters