From 848c748c974db677f2ba6000d879c18d2298237b Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 9 Nov 2026 13:16:59 -0800 Subject: [PATCH 1/3] External field excitations are applied in `OneStep_nosub` *after* `EvolveE`/`EvolveB` have already applied the field boundary conditions and the PEC mask. A hard (or soft) source whose flag region touches a PEC domain boundary, or a PEC layer from a GDS/numpy file, therefore overwrites them. For example, it can leave a nonzero tangential E on a PEC wall or along a metal trace. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This was the original motivation for #20. That PR is now stale (~1900 commits behind), and rebasing it as written would reintroduce soft-source double counting, because it moved the excitation inside the per-level/per-patch evolve routines. This PR fixes the issue without moving the excitation calls. At the end of `ApplyExternalFieldExcitationOnGrid(int externalfieldtype, DtType)`, the routine now re-enforces PEC on whatever was just excited: - **E** (`AllExternal`, `EfieldExternal`, `EfieldExternalPML`): - `ApplyEfieldBoundary(lev, PatchType::fine)`, i.e., domain PEC, including the split E fields in the PML - `ApplyPECMaskToEfield(lev, include_guard_cells=true)`, i.e., the GDS/numpy PEC mask on the grid (valid + guard cells, since the excitation also writes guard cells) and, if extended, in the PML - **B** (`AllExternal`, `BfieldExternal`): `PEC::ApplyPECtoBfield` only. The full `ApplyBfieldBoundary` is intentionally not called, because its Silver-Mueller update is not idempotent and must not run twice in a half step. Each step only runs when it applies: domain PEC when `PEC::isAnyBoundaryPEC()`, the mask when `use_PEC_mask` is on with the macroscopic solver, and only for fields whose excitation parser is active. Because the fix lives inside the excitation routine, it covers every call site: the t=0 `AllExternal` source, both B half-steps, E, E in the PML (`ApplyExcitationInPML`), and the LLG call in `SyncCurrentAndRho`. The existing `DtType` soft-source scaling is unchanged. New `WarpX::ApplyPECMaskToEfield(lev, include_guard_cells)`: - grid: `E *= PEC_fp` on the fine level, optionally including guard cells - PML (with `extend_PEC_mask_into_pml`): mask both split components; for the CPML, also reset the ψ_E memory variables on masked edges - no-op unless `use_PEC_mask` is on and `em_solver_medium == Macroscopic`. The mask is only initialized by `MacroscopicProperties`, so the old unconditional multiply in `OneStep` would have used an uninitialized mask with the vacuum solver. The two inline mask blocks (grid, then PML + CPML ψ) are replaced by a single `ApplyPECMaskToEfield(lev, false)` call after the PML update. It is still applied *before* `FillBoundaryE` / `PML::Exchange`, so guard cells and the PML interface plane see masked E. The PML E update does not read grid E, so applying the grid mask after the PML update instead of before is equivalent. The mask block after the E excitation is removed. It is now covered by `MacroscopicEvolveE` (no excitation) and by `ApplyExternalFieldExcitationOnGrid` (after excitation). The old block also only masked valid cells and did not touch the PML. `InitializePECFromSigma` only sets the valid region, so `PEC_fp` guard cells stayed 1. A `FillBoundary` on `PEC_fp[x,y,z]` now fills interior and periodic guard cells from the neighboring grids, so the mask can be applied to E guard cells. Guard cells outside the domain stay 1; the PML has its own mask. | Case | Before | After | |---|---|---| | Source touching a PEC domain boundary | excitation overwrites PEC | PEC enforced | | Source overlapping the GDS mask on the grid | masked (valid cells only); guard cells stale | masked, valid + guard cells | | Source overlapping the GDS mask in the PML (`ApplyExcitationInPML=1`) | excitation overwrites mask | masked; CPML ψ reset | | t=0 hard source overlapping PEC / mask | unmasked for the first B half-step | masked | - H / H_bias (LLG): there is no PEC boundary condition or mask for H, so a hard H source on a PEC wall or inside metal is not corrected. - B is not masked by the GDS PEC mask. - PML guard cells after a PML excitation: only valid PML cells are masked; guard cells are refreshed by the `FillBoundaryE` that follows `DampPML`. - Silver-Mueller boundaries are not re-applied after excitation, so a source on a Silver-Mueller boundary still overrides it for that step. - [x] Syntax check (`-fsyntax-only`, 3D, AMReX 25.07) with and without `WARPX_MAG_LLG` - [ ] Full build (GNUmake / CMake), with and without `USE_LLG` - [ ] GDS case with no excitation: results unchanged vs. `development` (up to roundoff) - [ ] Hard-source plane through a metal trace: tangential E on masked edges is exactly zero on the grid - [ ] Same with `ApplyExcitationInPML=1`, split-field PML and CPML: masked edges zero in the PML, ψ_E zero there - [ ] Hard-source plane intersecting a PEC domain boundary: tangential E on the wall stays zero - [ ] Multi-box run: masked E consistent across box boundaries (guard cells) - [ ] Existing regression tests unchanged Supersedes #20. --- .../WarpXFieldBoundaries.cpp | 39 +++++++++++++ Source/Evolve/WarpXEvolve.cpp | 15 +---- .../MacroscopicProperties.cpp | 7 +++ Source/FieldSolver/WarpXExternalEMFields.cpp | 36 ++++++++++++ Source/FieldSolver/WarpXPushFieldsEM.cpp | 58 +++++-------------- Source/WarpX.H | 12 ++++ 6 files changed, 110 insertions(+), 57 deletions(-) diff --git a/Source/BoundaryConditions/WarpXFieldBoundaries.cpp b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp index cd8e54229..1064aa424 100644 --- a/Source/BoundaryConditions/WarpXFieldBoundaries.cpp +++ b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp @@ -42,6 +42,45 @@ void WarpX::ApplyEfieldBoundary(const int lev, PatchType patch_type) } } +void WarpX::ApplyPECMaskToEfield (const int lev, const bool include_guard_cells) +{ + // The PEC mask is only initialized by the macroscopic solver (MacroscopicProperties) + if (!use_PEC_mask || em_solver_medium != MediumForEM::Macroscopic) { return; } + + // Grid: the guard cells of PEC_fp are filled at initialization + // (see MacroscopicProperties::InitData), so the mask can also be applied + // to the E guard cells, e.g., after an excitation has written into them. + for (int idim = 0; idim < 3; ++idim) { + const amrex::IntVect ng = include_guard_cells ? + amrex::min(Efield_fp[lev][idim]->nGrowVect(), PEC_fp[lev][idim]->nGrowVect()) : + amrex::IntVect::TheZeroVector(); + amrex::MultiFab::Multiply(*Efield_fp[lev][idim], *PEC_fp[lev][idim], 0, 0, 1, ng); + } + + // PML: mask extended into the PML (algo.extend_PEC_mask_into_pml) + if (do_pml && extend_PEC_mask_into_pml && pml[lev] && pml[lev]->ok()) { + const auto E_pml = pml[lev]->GetE_fp(); + const auto PEC_pml = pml[lev]->GetPEC_fp(); + + // 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. + for (int idim = 0; idim < 3; ++idim) { + amrex::MultiFab::Multiply(*E_pml[idim], *PEC_pml[idim], 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(); + for (int idim = 0; idim < 3; ++idim) { + amrex::MultiFab::Multiply(*psiE[idim], *PEC_pml[idim], 0, 0, 1, 0); + amrex::MultiFab::Multiply(*psiE[idim], *PEC_pml[idim], 0, 1, 1, 0); + } + } + } +} + void WarpX::ApplyBfieldBoundary (const int lev, PatchType patch_type, DtType a_dt_type) { if (PEC::isAnyBoundaryPEC()) { diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 348ae238f..ca5a6d2ba 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -540,18 +540,9 @@ WarpX::OneStep_nosub (Real cur_time) ApplyExternalFieldExcitationOnGrid(ExternalFieldType::EfieldExternalPML); } - if (use_PEC_mask) { - // multiply E field by PEC mask - amrex::MultiFab& Ex = *Efield_fp[0][0].get(); - amrex::MultiFab& Ey = *Efield_fp[0][1].get(); - amrex::MultiFab& Ez = *Efield_fp[0][2].get(); - amrex::MultiFab& PECx = *PEC_fp[0][0].get(); - amrex::MultiFab& PECy = *PEC_fp[0][1].get(); - amrex::MultiFab& PECz = *PEC_fp[0][2].get(); - MultiFab::Multiply(Ex, PECx, 0, 0, 1, 0); - MultiFab::Multiply(Ey, PECy, 0, 0, 1, 0); - MultiFab::Multiply(Ez, PECz, 0, 0, 1, 0); - } + // The PEC mask (use_PEC_mask) is applied in MacroscopicEvolveE, and again inside + // ApplyExternalFieldExcitationOnGrid after any E excitation (grid and PML), so + // masked edges are exactly zero at this point. EvolveF(0.5_rt * dt[0], DtType::SecondHalf); EvolveG(0.5_rt * dt[0], DtType::SecondHalf); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp index b2388c12a..2fcde2fc5 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp @@ -330,6 +330,13 @@ MacroscopicProperties::InitData () if (!m_sigma_npy_filename2.empty()) { InitializePECFromSigma(m_sigma_mf.get(), PECx, PECy, m_npy_k_index2); } + // InitializePECFromSigma only sets the valid region. Fill the guard cells from + // the neighboring grids so the mask can also be applied to the E guard cells + // (e.g., after an excitation, which writes into them). Guard cells outside the + // domain keep the value 1; the PML has its own mask. + PECx->FillBoundary(warpx.Geom(lev).periodicity()); + PECy->FillBoundary(warpx.Geom(lev).periodicity()); + PECz->FillBoundary(warpx.Geom(lev).periodicity()); // zero sigma m_sigma_mf->setVal(0.); } diff --git a/Source/FieldSolver/WarpXExternalEMFields.cpp b/Source/FieldSolver/WarpXExternalEMFields.cpp index 7ce721264..66637deac 100644 --- a/Source/FieldSolver/WarpXExternalEMFields.cpp +++ b/Source/FieldSolver/WarpXExternalEMFields.cpp @@ -1,5 +1,6 @@ #include "WarpX.H" #include "BoundaryConditions/PML.H" +#include "BoundaryConditions/WarpX_PEC.H" #include "Evolve/WarpXDtType.H" #include "Utils/WarpXConst.H" #include "Utils/WarpXUtil.H" @@ -106,6 +107,41 @@ WarpX::ApplyExternalFieldExcitationOnGrid (int const externalfieldtype, DtType a } #endif } // for loop over level + + // Re-apply the field boundary conditions, and the PEC mask, on any field that + // was just excited. The excitation is applied after EvolveE/EvolveB, which have + // already applied the boundary conditions and the PEC mask, so a hard (or soft) + // source whose flag region touches a PEC domain boundary or a PEC (GDS/numpy) + // layer would otherwise overwrite them (e.g., leave a nonzero tangential E on a + // PEC wall). Only the fine patch is excited, so only the fine patch is corrected. + const bool e_excited = + E_excitation_grid_s == "parse_e_excitation_grid_function" && + (externalfieldtype == ExternalFieldType::AllExternal || + externalfieldtype == ExternalFieldType::EfieldExternal || + externalfieldtype == ExternalFieldType::EfieldExternalPML); + const bool b_excited = + B_excitation_grid_s == "parse_b_excitation_grid_function" && + (externalfieldtype == ExternalFieldType::AllExternal || + externalfieldtype == ExternalFieldType::BfieldExternal); + for (int lev = 0; lev <= finest_level; ++lev) { + if (e_excited) { + if (PEC::isAnyBoundaryPEC()) { + // PEC on E, including the split E fields in the PML region + ApplyEfieldBoundary(lev, PatchType::fine); + } + // PEC mask on the grid (including the guard cells, which the excitation + // also writes) and, if extended, in the PML + const bool include_guard_cells = true; + ApplyPECMaskToEfield(lev, include_guard_cells); + } + if (b_excited && PEC::isAnyBoundaryPEC()) { + // PEC only: Silver-Mueller (in ApplyBfieldBoundary) advances B using E + // and must not be applied a second time within a half step + PEC::ApplyPECtoBfield( { get_pointer_Bfield_fp(lev, 0), + get_pointer_Bfield_fp(lev, 1), + get_pointer_Bfield_fp(lev, 2) }, lev, PatchType::fine); + } + } } void diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 223163d45..709811f1f 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -1125,21 +1125,6 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { 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()) { @@ -1207,38 +1192,21 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { pml[lev]->Getmu_cp(), pml[lev]->Getsigma_cp() ); } - - 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); - amrex::MultiFab* PECx_PEC = pml[lev]->GetPEC_fp(0); - amrex::MultiFab* PECy_PEC = pml[lev]->GetPEC_fp(1); - 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); - } - } - } + // Apply the PEC mask to the interior E (and, if extended, to the PML E and the + // CPML memory variables) right after the update, *before* the E guard cells and + // the PML are filled (FillBoundaryE -> PML::Exchange). If the mask were only + // applied after FillBoundaryE, the guard cells of the interior grids, and the PML + // copy of the interface plane (which the PML receives from the interior), would + // keep unmasked tangential E on the PEC edges until the next exchange, and the + // following B update would see 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 any E excitation (ApplyExternalFieldExcitationOnGrid), + // so masked edges stay exactly zero. + const bool include_guard_cells = false; + ApplyPECMaskToEfield(lev, include_guard_cells); + ApplyEfieldBoundary(lev, patch_type); } diff --git a/Source/WarpX.H b/Source/WarpX.H index f89a0e8b9..c08e696bc 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -931,6 +931,18 @@ public: void ApplyEfieldBoundary (const int lev, PatchType patch_type); void ApplyBfieldBoundary (const int lev, PatchType patch_type, DtType dt_type); + /** + * \brief With algo.use_PEC_mask (macroscopic solver only), zero the tangential E + * on the PEC edges of the numpy/GDS layer(s): on the fine patch of the grid and, + * with algo.extend_PEC_mask_into_pml, on the PML E fields (both split components) + * together with the CPML memory variables. + * + * \param[in] lev mesh refinement level + * \param[in] include_guard_cells also mask the guard cells of the grid E field + * (the PML is always masked on its valid region only) + */ + void ApplyPECMaskToEfield (const int lev, const bool include_guard_cells); + void DampPML (); void DampPML (const int lev); void DampPML (const int lev, PatchType patch_type); From 6a748478d069735ee63d5dbd0bac87eab2f909b3 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Tue, 10 Nov 2026 13:25:26 -0800 Subject: [PATCH 2/3] move excitation off PEC domain boundary, one cell inward. Regression benchmark will need to be changed --- Examples/Tests/circuits/inputs_regression_roger_pec | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Examples/Tests/circuits/inputs_regression_roger_pec b/Examples/Tests/circuits/inputs_regression_roger_pec index 9c8d45672..6d8c0c36e 100644 --- a/Examples/Tests/circuits/inputs_regression_roger_pec +++ b/Examples/Tests/circuits/inputs_regression_roger_pec @@ -134,12 +134,12 @@ warpx.Ez_excitation_flag_function(x,y,z) = " flag_hs * (z > -Lz/2 + th_si - tiny) * (z < -Lz/2 + th_si + th_nb + tiny) * (x > -w_line/2 - tiny) * (x < w_line/2 + tiny) + flag_hs * (z > -Lz/2 + th_si - tiny) * (z < -Lz/2 + th_si + th_nb + tiny) * (x > -w_line/2 - gap_cpw - w_gnd - tiny) * (x < -w_line/2 - gap_cpw + tiny) + flag_hs * (z > -Lz/2 + th_si - tiny) * (z < -Lz/2 + th_si + th_nb + tiny) * (x < +w_line/2 + gap_cpw + w_gnd + tiny) * (x > +w_line/2 + gap_cpw - tiny) -+ flag_hs * (x < w_port/2) * (x > -w_port/2) * (z < -Lz/2 + h_port) * (y > -Ly/2 - tiny) * (y < -Ly/2 + tiny)" ++ flag_hs * (x < w_port/2) * (x > -w_port/2) * (z < -Lz/2 + h_port) * (y > -Ly/2 + dy - tiny) * (y < -Ly/2 + dy + tiny)" warpx.Ex_excitation_grid_function(x,y,z,t) = "0." warpx.Ey_excitation_grid_function(x,y,z,t) = "0." warpx.Ez_excitation_grid_function(x,y,z,t) = "exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*c/wavelength*t) * (y > -Ly/2 - tiny) * (y < -Ly/2 + tiny)" -warpx.Ez_excitation_grid_function(x,y,z,t) = "sin(2*pi*c*t/wavelength) * (y > -Ly/2 - tiny) * (y < -Ly/2 + tiny)" +warpx.Ez_excitation_grid_function(x,y,z,t) = "sin(2*pi*c*t/wavelength) * (y > -Ly/2 - dy + tiny) * (y < -Ly/2 + dy + tiny)" # Diagnostics From 68cba94015ec510699ed0d1416b6a5e90a2f5e25 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Tue, 10 Nov 2026 13:31:44 -0800 Subject: [PATCH 3/3] restrict source area --- Examples/Tests/circuits/inputs_regression_roger_pec | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Examples/Tests/circuits/inputs_regression_roger_pec b/Examples/Tests/circuits/inputs_regression_roger_pec index 6d8c0c36e..06547f9a0 100644 --- a/Examples/Tests/circuits/inputs_regression_roger_pec +++ b/Examples/Tests/circuits/inputs_regression_roger_pec @@ -138,8 +138,7 @@ warpx.Ez_excitation_flag_function(x,y,z) = " warpx.Ex_excitation_grid_function(x,y,z,t) = "0." warpx.Ey_excitation_grid_function(x,y,z,t) = "0." -warpx.Ez_excitation_grid_function(x,y,z,t) = "exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*c/wavelength*t) * (y > -Ly/2 - tiny) * (y < -Ly/2 + tiny)" -warpx.Ez_excitation_grid_function(x,y,z,t) = "sin(2*pi*c*t/wavelength) * (y > -Ly/2 - dy + tiny) * (y < -Ly/2 + dy + tiny)" +warpx.Ez_excitation_grid_function(x,y,z,t) = "sin(2*pi*c*t/wavelength) * (y > -Ly/2 - dy + tiny) * (y < -Ly/2 + dy + tiny) * (x < w_port/2) * (x > -w_port/2) * (z < -Lz/2 + h_port)" # Diagnostics