Enforce PEC boundaries and the GDS PEC mask after external field excitation

Examples/Tests/circuits/inputs_regression_roger_pec

@@ -134,12 +134,11 @@ 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) * (x < w_port/2) * (x > -w_port/2) * (z < -Lz/2 + h_port)"
 
 # Diagnostics
 

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()) {

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);

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.);
     }

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

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);
 }
 

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);