diff --git a/Examples/Tests/circuits/inputs_regression_roger_pec b/Examples/Tests/circuits/inputs_regression_roger_pec index 9c8d45672..06547f9a0 100644 --- a/Examples/Tests/circuits/inputs_regression_roger_pec +++ b/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 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);