From 5e728e3df135a70e99c3a7c797a7c49b5aa88992 Mon Sep 17 00:00:00 2001 From: Justin Angus Date: Sat, 8 Feb 2025 12:02:21 -0800 Subject: [PATCH 1/4] making BC for J at PEC and PMC (symmetry) boundaries consistent with BC for E as required by energy conservation. --- .../WarpXFieldBoundaries.cpp | 16 +- Source/BoundaryConditions/WarpX_PEC.H | 12 +- Source/BoundaryConditions/WarpX_PEC.cpp | 155 ++++++++---------- 3 files changed, 82 insertions(+), 101 deletions(-) diff --git a/Source/BoundaryConditions/WarpXFieldBoundaries.cpp b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp index 6217eb04a33..fd24e73787e 100644 --- a/Source/BoundaryConditions/WarpXFieldBoundaries.cpp +++ b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp @@ -266,14 +266,12 @@ void WarpX::ApplyBfieldBoundary (const int lev, PatchType patch_type, DtType a_d void WarpX::ApplyRhofieldBoundary (const int lev, MultiFab* rho, PatchType patch_type) { - if (::isAnyBoundary(particle_boundary_lo, particle_boundary_hi) || - ::isAnyBoundary(particle_boundary_lo, particle_boundary_hi) || - ::isAnyBoundary(field_boundary_lo, field_boundary_hi) || + if (::isAnyBoundary(field_boundary_lo, field_boundary_hi) || ::isAnyBoundary(field_boundary_lo, field_boundary_hi)) { - PEC::ApplyReflectiveBoundarytoRhofield(rho, + // This routine handles rho at both PEC and PMC boundaries + PEC::ApplyBoundarytoRhofield(rho, field_boundary_lo, field_boundary_hi, - particle_boundary_lo, particle_boundary_hi, Geom(lev), lev, patch_type, ref_ratio); } } @@ -282,14 +280,12 @@ void WarpX::ApplyJfieldBoundary (const int lev, amrex::MultiFab* Jx, amrex::MultiFab* Jy, amrex::MultiFab* Jz, PatchType patch_type) { - if (::isAnyBoundary(particle_boundary_lo, particle_boundary_hi) || - ::isAnyBoundary(particle_boundary_lo, particle_boundary_hi) || - ::isAnyBoundary(field_boundary_lo, field_boundary_hi) || + if (::isAnyBoundary(field_boundary_lo, field_boundary_hi) || ::isAnyBoundary(field_boundary_lo, field_boundary_hi)) { - PEC::ApplyReflectiveBoundarytoJfield(Jx, Jy, Jz, + // This routine handles J at both PEC and PMC boundaries + PEC::ApplyBoundarytoJfield(Jx, Jy, Jz, field_boundary_lo, field_boundary_hi, - particle_boundary_lo, particle_boundary_hi, Geom(lev), lev, patch_type, ref_ratio); } } diff --git a/Source/BoundaryConditions/WarpX_PEC.H b/Source/BoundaryConditions/WarpX_PEC.H index e3fd804b62c..df2a2bf7019 100644 --- a/Source/BoundaryConditions/WarpX_PEC.H +++ b/Source/BoundaryConditions/WarpX_PEC.H @@ -61,7 +61,7 @@ namespace PEC { bool split_pml_field = false); /** - * \brief Reflects charge density deposited over the PEC boundary back into + * \brief Reflects charge density deposited over the PEC or PMC boundary back into * the simulation domain. * * \param[in,out] rho Multifab containing the charge density @@ -74,17 +74,15 @@ namespace PEC { * \param[in] patch_type coarse or fine * \param[in] ref_ratios vector containing the refinement ratios of the refinement levels */ - void ApplyReflectiveBoundarytoRhofield( + void ApplyBoundarytoRhofield( amrex::MultiFab* rho, const amrex::Array& field_boundary_lo, const amrex::Array& field_boundary_hi, - const amrex::Array& particle_boundary_lo, - const amrex::Array& particle_boundary_hi, const amrex::Geometry& geom, int lev, PatchType patch_type, const amrex::Vector& ref_ratios); /** - * \brief Reflects current density deposited over the PEC boundary back into + * \brief Reflects current density deposited over the PEC or PMC boundary back into * the simulation domain. * * \param[in,out] Jx, Jy, Jz Multifabs containing the current density @@ -97,13 +95,11 @@ namespace PEC { * \param[in] patch_type coarse or fine * \param[in] ref_ratios vector containing the refinement ratios of the refinement levels */ - void ApplyReflectiveBoundarytoJfield( + void ApplyBoundarytoJfield( amrex::MultiFab* Jx, amrex::MultiFab* Jy, amrex::MultiFab* Jz, const amrex::Array& field_boundary_lo, const amrex::Array& field_boundary_hi, - const amrex::Array& particle_boundary_lo, - const amrex::Array& particle_boundary_hi, const amrex::Geometry& geom, int lev, PatchType patch_type, const amrex::Vector& ref_ratios); diff --git a/Source/BoundaryConditions/WarpX_PEC.cpp b/Source/BoundaryConditions/WarpX_PEC.cpp index a3b75791582..16b070d82e8 100644 --- a/Source/BoundaryConditions/WarpX_PEC.cpp +++ b/Source/BoundaryConditions/WarpX_PEC.cpp @@ -152,8 +152,8 @@ namespace const bool is_tangent_to_PEC = (icomp != idim); #endif if (isPECBoundary) { - // grid point ijk_vec is ig number of points pass the - // domain boundary in direction, idim + // grid point ijk_vec is ig number of points past the + // domain boundary in idim direction const int ig = ::get_cell_count_to_boundary( dom_lo, dom_hi, ijk_vec, is_nodal, idim, iside); @@ -333,25 +333,36 @@ namespace /** - * \brief Sets the rho or J field value in cells close to and on reflecting particle boundary - * or PEC field boundary. The charge/current density deposited - * in the guard cells are either reflected - * back into the simulation domain (if a reflecting particle - * boundary is used), or the opposite charge/current density is deposited - * back in the domain to capture the effect of an image charge. - * The charge/current density on the reflecting boundary is set to 0 while values - * in the guard cells are set equal (and opposite) to their mirror - * location inside the domain - representing image charges - in the - * normal (tangential) direction. + * \brief Sets the rho or J field value near PEC and PMC boundaries by appropriately + * folding in the rho or J deposited into the ghost cells. + * Note that energy conservation demands that the boundary conditions for J be + * consistent with those for E used to push the particles. If E_i is odd, then + * J_i deposited to ghost cells should be subtracted from its mirror location + * inside the domain. If E_i is even, then J_i deposited to the ghost cells should + * be added to its mirror location inside the domain. If this is not done, + * then artificial heating can occur at the boundaries. + * PEC: The charge and current density parallel to the boundary deposited to the + * guard cells is subtracted from its mirror location inside the domain. + * The current density normal to the boundary deposited to the guard cells is + * added to its mirror location inside the domain. + * This can be interpreted as depositing the charge/current from the oppositely + * charged mirror particle inside the PEC. It can also be interpreted as using + * a lower order shape factor for the particles near the boundary. + * PMC: The charge and current density parallel to the boundary deposited to the + * guard cells is added to its mirror location inside the domain. + * The current density normal to the boundary deposited to the guard cells is + * subtracted from its mirror location inside the domain. Since PMC is a + * symmetry boundary, this charge and current comes from the physical + * mirror particle with the same charge on the other side of the boundary. * - * \param[in] n index of the MultiFab component being updated - * \param[in] ijk_vec indices along the x(i), y(j), z(k) of the rho Array4 - * \param[in out] field field data to be updated - * \param[in] mirrorfac mirror cell is given by mirrorfac - ijk_vec - * \param[in] psign Whether the field value should be flipped across the boundary - * \param[in] is_reflective Whether the given particle boundary is reflecting or field boundary is pec - * \param[in] tangent_to_bndy Whether a given direction is perpendicular to the boundary - * \param[in] fabbox multifab box including ghost cells + * \param[in] n index of the MultiFab component being updated + * \param[in] ijk_vec indices along the x(i), y(j), z(k) of the rho Array4 + * \param[in out] field field data to be updated + * \param[in] mirrorfac mirror cell is given by mirrorfac - ijk_vec + * \param[in] psign Whether the field value should be flipped across the boundary + * \param[in] is_pec_pmc_bndy Whether the field boundary is pec or pmc + * \param[in] tangent_to_bndy Whether a given direction is perpendicular to the boundary + * \param[in] fabbox multifab box including ghost cells */ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void SetRhoOrJfieldFromPEC (const int n, @@ -359,7 +370,7 @@ namespace amrex::Array4 const& field, amrex::GpuArray, AMREX_SPACEDIM> const& mirrorfac, amrex::GpuArray, AMREX_SPACEDIM> const& psign, - amrex::GpuArray, AMREX_SPACEDIM> const& is_reflective, + amrex::GpuArray, AMREX_SPACEDIM> const& is_pec_pmc_bndy, amrex::GpuArray const& tangent_to_bndy, amrex::Box const& fabbox) { @@ -370,14 +381,14 @@ namespace { for (int iside = 0; iside < 2; ++iside) { - if (!is_reflective[idim][iside]) { continue; } + if (!is_pec_pmc_bndy[idim][iside]) { continue; } // Get the mirror guard cell index amrex::IntVect iv_mirror = ijk_vec; iv_mirror[idim] = mirrorfac[idim][iside] - ijk_vec[idim]; // Update the cell if the mirror guard cell exists - if (ijk_vec == iv_mirror && is_reflective[idim][iside] == 1) { + if (ijk_vec == iv_mirror && psign[idim][iside] == -1) { field(ijk_vec,n) = 0._rt; } else if (fabbox.contains(iv_mirror)) { // Note that this includes the cells on the boundary for PMC @@ -391,7 +402,7 @@ namespace { for (int iside = 0; iside < 2; ++iside) { - if (!is_reflective[idim][iside]) { continue; } + if (!is_pec_pmc_bndy[idim][iside]) { continue; } amrex::IntVect iv_mirror = ijk_vec; iv_mirror[idim] = mirrorfac[idim][iside] - ijk_vec[idim]; @@ -633,12 +644,10 @@ PEC::ApplyPECtoBfield ( * location inside the domain - representing image charges. **/ void -PEC::ApplyReflectiveBoundarytoRhofield ( +PEC::ApplyBoundarytoRhofield ( amrex::MultiFab* rho, const amrex::Array& field_boundary_lo, const amrex::Array& field_boundary_hi, - const amrex::Array& particle_boundary_lo, - const amrex::Array& particle_boundary_hi, const amrex::Geometry& geom, const int lev, PatchType patch_type, const amrex::Vector& ref_ratios) { @@ -658,34 +667,28 @@ PEC::ApplyReflectiveBoundarytoRhofield ( // cells for boundaries that are NOT PEC amrex::Box grown_domain_box = domain_box; - amrex::GpuArray, AMREX_SPACEDIM> is_reflective; + amrex::GpuArray, AMREX_SPACEDIM> is_pec_pmc_bndy; amrex::GpuArray is_tangent_to_bndy; amrex::GpuArray, AMREX_SPACEDIM> psign; amrex::GpuArray, AMREX_SPACEDIM> mirrorfac; for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { - is_reflective[idim][0] = ( particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting) - || ( particle_boundary_lo[idim] == ParticleBoundaryType::Thermal) - || ( field_boundary_lo[idim] == FieldBoundaryType::PEC); - if (field_boundary_lo[idim] == FieldBoundaryType::PMC) { is_reflective[idim][0] = 2; } - is_reflective[idim][1] = ( particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting) - || ( particle_boundary_hi[idim] == ParticleBoundaryType::Thermal) - || ( field_boundary_hi[idim] == FieldBoundaryType::PEC); - if (field_boundary_hi[idim] == FieldBoundaryType::PMC) { is_reflective[idim][1] = 2; } - if (!is_reflective[idim][0]) { grown_domain_box.growLo(idim, ng_fieldgather[idim]); } - if (!is_reflective[idim][1]) { grown_domain_box.growHi(idim, ng_fieldgather[idim]); } + is_pec_pmc_bndy[idim][0] = ( field_boundary_lo[idim] == FieldBoundaryType::PEC + || field_boundary_lo[idim] == FieldBoundaryType::PMC ); + is_pec_pmc_bndy[idim][1] = ( field_boundary_hi[idim] == FieldBoundaryType::PEC + || field_boundary_lo[idim] == FieldBoundaryType::PMC ); + // What are the below lines for? Why grow box to include ghosts only if not reflective? + // What does E and B BC routines do if the boundary is not a PEC or PMC boundary? + if (!is_pec_pmc_bndy[idim][0]) { grown_domain_box.growLo(idim, ng_fieldgather[idim]); } + if (!is_pec_pmc_bndy[idim][1]) { grown_domain_box.growHi(idim, ng_fieldgather[idim]); } // rho values inside guard cells are updated the same as tangential // components of the current density is_tangent_to_bndy[idim] = true; - psign[idim][0] = ((particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting) - ||(particle_boundary_lo[idim] == ParticleBoundaryType::Thermal) - ||(field_boundary_lo[idim] == FieldBoundaryType::PMC)) - ? 1._rt : -1._rt; - psign[idim][1] = ((particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting) - ||(particle_boundary_hi[idim] == ParticleBoundaryType::Thermal) - ||(field_boundary_hi[idim] == FieldBoundaryType::PMC)) - ? 1._rt : -1._rt; + // Default is -1 for PEC boundary + psign[idim][0] = (field_boundary_lo[idim] == FieldBoundaryType::PMC ? 1._rt : -1._rt); + psign[idim][1] = (field_boundary_hi[idim] == FieldBoundaryType::PMC ? 1._rt : -1._rt); + mirrorfac[idim][0] = 2*domain_lo[idim] - (1 - rho_nodal[idim]); mirrorfac[idim][1] = 2*domain_hi[idim] + (1 - rho_nodal[idim]); } @@ -701,6 +704,7 @@ PEC::ApplyReflectiveBoundarytoRhofield ( // If grown_domain_box contains fabbox it means there are no PEC // boundaries to handle so continue to next box + // Is this because it is grown above? Don't follow. if (grown_domain_box.contains(fabbox)) { continue; } // Extract field data @@ -714,7 +718,7 @@ PEC::ApplyReflectiveBoundarytoRhofield ( const amrex::IntVect iv(AMREX_D_DECL(i,j,k)); ::SetRhoOrJfieldFromPEC( - n, iv, rho_array, mirrorfac, psign, is_reflective, + n, iv, rho_array, mirrorfac, psign, is_pec_pmc_bndy, is_tangent_to_bndy, fabbox ); }); @@ -723,12 +727,10 @@ PEC::ApplyReflectiveBoundarytoRhofield ( void -PEC::ApplyReflectiveBoundarytoJfield( +PEC::ApplyBoundarytoJfield( amrex::MultiFab* Jx, amrex::MultiFab* Jy, amrex::MultiFab* Jz, const amrex::Array& field_boundary_lo, const amrex::Array& field_boundary_hi, - const amrex::Array& particle_boundary_lo, - const amrex::Array& particle_boundary_hi, const amrex::Geometry& geom, const int lev, PatchType patch_type, const amrex::Vector& ref_ratios) { @@ -758,23 +760,19 @@ PEC::ApplyReflectiveBoundarytoJfield( // directions of the current density multifab const amrex::IntVect ng_fieldgather = Jx->nGrowVect(); - amrex::GpuArray, AMREX_SPACEDIM> is_reflective; + amrex::GpuArray, AMREX_SPACEDIM> is_pec_pmc_bndy; amrex::GpuArray, 3> is_tangent_to_bndy; amrex::GpuArray, AMREX_SPACEDIM>, 3> psign; amrex::GpuArray, AMREX_SPACEDIM>, 3> mirrorfac; for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { - is_reflective[idim][0] = ( particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting) - || ( particle_boundary_lo[idim] == ParticleBoundaryType::Thermal) - || ( field_boundary_lo[idim] == FieldBoundaryType::PEC) - || ( field_boundary_lo[idim] == FieldBoundaryType::PMC); - if (field_boundary_lo[idim] == FieldBoundaryType::PMC) { is_reflective[idim][0] = 2; } - is_reflective[idim][1] = ( particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting) - || ( particle_boundary_hi[idim] == ParticleBoundaryType::Thermal) - || ( field_boundary_hi[idim] == FieldBoundaryType::PEC) - || ( field_boundary_hi[idim] == FieldBoundaryType::PMC); - if (field_boundary_hi[idim] == FieldBoundaryType::PMC) { is_reflective[idim][1] = 2; } - if (!is_reflective[idim][0]) { grown_domain_box.growLo(idim, ng_fieldgather[idim]); } - if (!is_reflective[idim][1]) { grown_domain_box.growHi(idim, ng_fieldgather[idim]); } + is_pec_pmc_bndy[idim][0] = ( field_boundary_lo[idim] == FieldBoundaryType::PEC + || field_boundary_lo[idim] == FieldBoundaryType::PMC ); + is_pec_pmc_bndy[idim][1] = ( field_boundary_hi[idim] == FieldBoundaryType::PEC + || field_boundary_hi[idim] == FieldBoundaryType::PMC ); + // What are the below lines for? Why grow box to include ghosts only if not reflective? + // What does E and B BC routines do if the boundary is not a PEC or PMC boundary? + if (!is_pec_pmc_bndy[idim][0]) { grown_domain_box.growLo(idim, ng_fieldgather[idim]); } + if (!is_pec_pmc_bndy[idim][1]) { grown_domain_box.growHi(idim, ng_fieldgather[idim]); } for (int icomp=0; icomp < 3; ++icomp) { // Set the psign value correctly for each current component for each @@ -793,24 +791,14 @@ PEC::ApplyReflectiveBoundarytoJfield( #endif if (is_tangent_to_bndy[icomp][idim]){ - psign[icomp][idim][0] = ( (particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting) - ||(particle_boundary_lo[idim] == ParticleBoundaryType::Thermal) - ||(field_boundary_lo[idim] == FieldBoundaryType::PMC)) - ? 1._rt : -1._rt; - psign[icomp][idim][1] = ( (particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting) - ||(particle_boundary_hi[idim] == ParticleBoundaryType::Thermal) - ||(field_boundary_hi[idim] == FieldBoundaryType::PMC)) - ? 1._rt : -1._rt; + // Default is -1 for PEC boundary + psign[icomp][idim][0] = ( field_boundary_lo[idim] == FieldBoundaryType::PMC ? 1._rt : -1._rt ); + psign[icomp][idim][1] = ( field_boundary_hi[idim] == FieldBoundaryType::PMC ? 1._rt : -1._rt ); } else { - psign[icomp][idim][0] = ( (particle_boundary_lo[idim] == ParticleBoundaryType::Reflecting) - ||(particle_boundary_lo[idim] == ParticleBoundaryType::Thermal) - ||(field_boundary_lo[idim] == FieldBoundaryType::PMC)) - ? -1._rt : 1._rt; - psign[icomp][idim][1] = ( (particle_boundary_hi[idim] == ParticleBoundaryType::Reflecting) - ||(particle_boundary_hi[idim] == ParticleBoundaryType::Thermal) - ||(field_boundary_hi[idim] == FieldBoundaryType::PMC)) - ? -1._rt : 1._rt; + // Default is +1 for PEC boundary + psign[icomp][idim][0] = ( field_boundary_lo[idim] == FieldBoundaryType::PMC ? -1._rt : 1._rt ); + psign[icomp][idim][1] = ( field_boundary_hi[idim] == FieldBoundaryType::PMC ? -1._rt : 1._rt ); } } // Set the correct mirror cell calculation for each current @@ -836,7 +824,8 @@ PEC::ApplyReflectiveBoundarytoJfield( Box const& fabbox = mfi.fabbox(); // If grown_domain_box contains fabbox it means there are no PEC - // boundaries to handle so continue to next box + // boundaries to handle so continue to next box. + // Is this because it is grown above? Don't follow. if (grown_domain_box.contains(fabbox)) { continue; } // Extract field data @@ -850,7 +839,7 @@ PEC::ApplyReflectiveBoundarytoJfield( const amrex::IntVect iv(AMREX_D_DECL(i,j,k)); ::SetRhoOrJfieldFromPEC( - n, iv, Jx_array, mirrorfac[0], psign[0], is_reflective, + n, iv, Jx_array, mirrorfac[0], psign[0], is_pec_pmc_bndy, is_tangent_to_bndy[0], fabbox ); }); @@ -881,7 +870,7 @@ PEC::ApplyReflectiveBoundarytoJfield( const amrex::IntVect iv(AMREX_D_DECL(i,j,k)); ::SetRhoOrJfieldFromPEC( - n, iv, Jy_array, mirrorfac[1], psign[1], is_reflective, + n, iv, Jy_array, mirrorfac[1], psign[1], is_pec_pmc_bndy, is_tangent_to_bndy[1], fabbox ); }); @@ -912,7 +901,7 @@ PEC::ApplyReflectiveBoundarytoJfield( const amrex::IntVect iv(AMREX_D_DECL(i,j,k)); ::SetRhoOrJfieldFromPEC( - n, iv, Jz_array, mirrorfac[2], psign[2], is_reflective, + n, iv, Jz_array, mirrorfac[2], psign[2], is_pec_pmc_bndy, is_tangent_to_bndy[2], fabbox ); }); From 4b9cfb83daea6ef22a3e731dd8feecc9dd1c14fe Mon Sep 17 00:00:00 2001 From: Justin Angus Date: Sun, 9 Feb 2025 08:10:28 -0800 Subject: [PATCH 2/4] updating briefs. --- Source/BoundaryConditions/WarpX_PEC.cpp | 44 ++++++++++++++++--------- 1 file changed, 28 insertions(+), 16 deletions(-) diff --git a/Source/BoundaryConditions/WarpX_PEC.cpp b/Source/BoundaryConditions/WarpX_PEC.cpp index 16b070d82e8..98eb86e6acf 100644 --- a/Source/BoundaryConditions/WarpX_PEC.cpp +++ b/Source/BoundaryConditions/WarpX_PEC.cpp @@ -632,16 +632,18 @@ PEC::ApplyPECtoBfield ( } } - /** - * \brief Sets the rho field value in cells close to and inside a PEC boundary. - * The charge density deposited in the guard cells are either reflected - * back into the simulation domain (if a reflecting particle - * boundary is used), or the opposite charge density is deposited - * back in the domain to capture the effect of an image charge. - * The charge density on the PEC boundary is set to 0 while values - * in the guard cells are set equal and opposite to their mirror - * location inside the domain - representing image charges. + * \brief Sets the rho field value in cells close to and inside PEC/PMC boundaries. + * The charge density deposited in the guard cells at PMC boundaries is added + * to the mirror location in the simulation domain to capture rho from image + * charge of the same sign on the other side of the symmetry boundary. + * The charge density deposited in the guard cells at PEC boundaries is subtracted + * from the mirror location in the simulation domain to capture rho from image + * charges of the opposite sign on the conductor boundary. + * For PEC boundaries, this method can also be interpretted as using a lower + * order shape function for particles near the boundary. + * Note that the boundary handling of rho at PEC and PMC boundaries is exactly + * the same as that for the current density parallel to the boundary. **/ void PEC::ApplyBoundarytoRhofield ( @@ -676,8 +678,6 @@ PEC::ApplyBoundarytoRhofield ( || field_boundary_lo[idim] == FieldBoundaryType::PMC ); is_pec_pmc_bndy[idim][1] = ( field_boundary_hi[idim] == FieldBoundaryType::PEC || field_boundary_lo[idim] == FieldBoundaryType::PMC ); - // What are the below lines for? Why grow box to include ghosts only if not reflective? - // What does E and B BC routines do if the boundary is not a PEC or PMC boundary? if (!is_pec_pmc_bndy[idim][0]) { grown_domain_box.growLo(idim, ng_fieldgather[idim]); } if (!is_pec_pmc_bndy[idim][1]) { grown_domain_box.growHi(idim, ng_fieldgather[idim]); } @@ -704,7 +704,6 @@ PEC::ApplyBoundarytoRhofield ( // If grown_domain_box contains fabbox it means there are no PEC // boundaries to handle so continue to next box - // Is this because it is grown above? Don't follow. if (grown_domain_box.contains(fabbox)) { continue; } // Extract field data @@ -725,7 +724,23 @@ PEC::ApplyBoundarytoRhofield ( } } - +/** + * \brief Sets the J field value in cells close to and inside PEC/PMC boundaries. + * The proper way to reflect current density deposited to ghost cells back + * into the simulation domain, as determined by conservation of energy, depends + * on the boundary condition used for the electric field used in the gather + * to push the particles. If E_i is odd, then J_i deposited to ghost cells is + * subtracted from its mirror location inside the domain. Likewise, if E_i is + * even, then J_i deposited in the ghost cells is added to its mirror location + * inside the domain. + * For PEC boundaries, E_parallel is odd and E_perp is even. + * For PMC boundaries, E_parallel is even and E_perp is odd. + * For PMC boundaries, this J comes from the mirror charge of the same sign + * on the other side of the symmetry boundary. + * For PEC boundaries, this J comes from the mirror charge of the opposite + * sign on the conductor boundary, and can also be interpretted as using a + * lower order shape function for particles near the boundary. + **/ void PEC::ApplyBoundarytoJfield( amrex::MultiFab* Jx, amrex::MultiFab* Jy, amrex::MultiFab* Jz, @@ -769,8 +784,6 @@ PEC::ApplyBoundarytoJfield( || field_boundary_lo[idim] == FieldBoundaryType::PMC ); is_pec_pmc_bndy[idim][1] = ( field_boundary_hi[idim] == FieldBoundaryType::PEC || field_boundary_hi[idim] == FieldBoundaryType::PMC ); - // What are the below lines for? Why grow box to include ghosts only if not reflective? - // What does E and B BC routines do if the boundary is not a PEC or PMC boundary? if (!is_pec_pmc_bndy[idim][0]) { grown_domain_box.growLo(idim, ng_fieldgather[idim]); } if (!is_pec_pmc_bndy[idim][1]) { grown_domain_box.growHi(idim, ng_fieldgather[idim]); } @@ -825,7 +838,6 @@ PEC::ApplyBoundarytoJfield( // If grown_domain_box contains fabbox it means there are no PEC // boundaries to handle so continue to next box. - // Is this because it is grown above? Don't follow. if (grown_domain_box.contains(fabbox)) { continue; } // Extract field data From e0a8dbbf21f09e3a8a2e5e5640eb55ee1283b00c Mon Sep 17 00:00:00 2001 From: Justin Angus Date: Fri, 14 Feb 2025 09:00:54 -0800 Subject: [PATCH 3/4] removing if condition --- Source/BoundaryConditions/WarpX_PEC.cpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/Source/BoundaryConditions/WarpX_PEC.cpp b/Source/BoundaryConditions/WarpX_PEC.cpp index 98eb86e6acf..112f2fed4e5 100644 --- a/Source/BoundaryConditions/WarpX_PEC.cpp +++ b/Source/BoundaryConditions/WarpX_PEC.cpp @@ -387,13 +387,9 @@ namespace amrex::IntVect iv_mirror = ijk_vec; iv_mirror[idim] = mirrorfac[idim][iside] - ijk_vec[idim]; - // Update the cell if the mirror guard cell exists - if (ijk_vec == iv_mirror && psign[idim][iside] == -1) { - field(ijk_vec,n) = 0._rt; - } else if (fabbox.contains(iv_mirror)) { - // Note that this includes the cells on the boundary for PMC - field(ijk_vec,n) += psign[idim][iside] * field(iv_mirror,n); - } + // Note that this includes the cells on the boundary + field(ijk_vec,n) += psign[idim][iside] * field(iv_mirror,n); + } } // 2) The guard cells are updated with the appropriate image From 917d0d8f4892dffb7df84477dc4c9d12ab940180 Mon Sep 17 00:00:00 2001 From: Justin Angus Date: Fri, 14 Feb 2025 09:41:31 -0800 Subject: [PATCH 4/4] fixing bug from last commit. --- Source/BoundaryConditions/WarpX_PEC.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/BoundaryConditions/WarpX_PEC.cpp b/Source/BoundaryConditions/WarpX_PEC.cpp index 112f2fed4e5..a76adc7db74 100644 --- a/Source/BoundaryConditions/WarpX_PEC.cpp +++ b/Source/BoundaryConditions/WarpX_PEC.cpp @@ -388,7 +388,9 @@ namespace iv_mirror[idim] = mirrorfac[idim][iside] - ijk_vec[idim]; // Note that this includes the cells on the boundary - field(ijk_vec,n) += psign[idim][iside] * field(iv_mirror,n); + if (fabbox.contains(iv_mirror)) { + field(ijk_vec,n) += psign[idim][iside] * field(iv_mirror,n); + } } }