From f6d3098f28a91fb0e75bc3ccf1ce0f8172efb58b Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Mon, 7 Jul 2025 15:31:38 -0400 Subject: [PATCH] > SR and GR: bug in deposit: wrong jx3 in 2D > GR: bug in gr range_with_axis_BCs() > GR: symmetric BC for axes > GR: small angle check for qks metric --- src/engines/grpic.hpp | 2 +- src/kernels/currents_deposit.hpp | 2 +- src/kernels/fields_bcs.hpp | 8 ++++++-- src/metrics/qkerr_schild.h | 29 ++++++++++++++++++++++------- 4 files changed, 30 insertions(+), 11 deletions(-) diff --git a/src/engines/grpic.hpp b/src/engines/grpic.hpp index 168f2c73..cd04838d 100644 --- a/src/engines/grpic.hpp +++ b/src/engines/grpic.hpp @@ -840,7 +840,7 @@ namespace ntt { if constexpr (M::Dim == Dim::_2D) { if (domain.mesh.flds_bc_in({ 0, +1 }) == FldsBC::AXIS) { range = CreateRangePolicy( - { domain.mesh.i_min(in::x1) - 1, domain.mesh.i_min(in::x2) }, + { domain.mesh.i_min(in::x1), domain.mesh.i_min(in::x2) }, { domain.mesh.i_max(in::x1), domain.mesh.i_max(in::x2) + 1 }); } } else if constexpr (M::Dim == Dim::_3D) { diff --git a/src/kernels/currents_deposit.hpp b/src/kernels/currents_deposit.hpp index 98d00a9b..e0a72e56 100644 --- a/src/kernels/currents_deposit.hpp +++ b/src/kernels/currents_deposit.hpp @@ -232,7 +232,7 @@ namespace kernel { cur::jx3) += Fx3_1 * (ONE - Wx1_1) * (ONE - Wx2_1); J_acc(i1_prev(p) + N_GHOSTS + 1, i2_prev(p) + N_GHOSTS, - cur::jx3) += Fx3_1 * Wx1_2 * (ONE - Wx2_1); + cur::jx3) += Fx3_1 * Wx1_1 * (ONE - Wx2_1); J_acc(i1_prev(p) + N_GHOSTS, i2_prev(p) + N_GHOSTS + 1, cur::jx3) += Fx3_1 * (ONE - Wx1_1) * Wx2_1; diff --git a/src/kernels/fields_bcs.hpp b/src/kernels/fields_bcs.hpp index 04e63052..4891688a 100644 --- a/src/kernels/fields_bcs.hpp +++ b/src/kernels/fields_bcs.hpp @@ -843,21 +843,25 @@ namespace kernel::bc { if constexpr (not P) { if (setE) { Fld(i1, i_edge - 1, em::ex2) = -Fld(i1, i_edge, em::ex2); - Fld(i1, i_edge, em::ex3) = ZERO; + Fld(i1, i_edge, em::ex3) = ZERO; + Fld(i1, i_edge - 1, em::ex3) = Fld(i1, i_edge + 1, em::ex3); } if (setB) { Fld(i1, i_edge - 1, em::bx1) = Fld(i1, i_edge, em::bx1); Fld(i1, i_edge, em::bx2) = ZERO; + Fld(i1, i_edge - 1, em::bx2) = - Fld(i1, i_edge + 1, em::bx2); Fld(i1, i_edge - 1, em::bx3) = Fld(i1, i_edge, em::bx3); } } else { if (setE) { Fld(i1, i_edge, em::ex2) = -Fld(i1, i_edge - 1, em::ex2); - Fld(i1, i_edge, em::ex3) = ZERO; + Fld(i1, i_edge, em::ex3) = ZERO; + Fld(i1, i_edge + 1, em::ex3) = Fld(i1, i_edge - 1, em::ex3); } if (setB) { Fld(i1, i_edge, em::bx1) = Fld(i1, i_edge - 1, em::bx1); Fld(i1, i_edge, em::bx2) = ZERO; + Fld(i1, i_edge + 1, em::bx2) = - Fld(i1, i_edge - 1, em::bx2); Fld(i1, i_edge, em::bx3) = Fld(i1, i_edge - 1, em::bx3); } } diff --git a/src/metrics/qkerr_schild.h b/src/metrics/qkerr_schild.h index efe44611..6ff114b4 100644 --- a/src/metrics/qkerr_schild.h +++ b/src/metrics/qkerr_schild.h @@ -39,6 +39,7 @@ namespace metric { const real_t chi_min, eta_min, phi_min; const real_t dchi, deta, dphi; const real_t dchi_inv, deta_inv, dphi_inv; + const bool small_angle; Inline auto Delta(const real_t& r) const -> real_t { return SQR(r) - TWO * r + SQR(a); @@ -89,7 +90,8 @@ namespace metric { , dphi { (x3_max - phi_min) / nx3 } , dchi_inv { ONE / dchi } , deta_inv { ONE / deta } - , dphi_inv { ONE / dphi } { + , dphi_inv { ONE / dphi } + , small_angle { eta2theta(HALF * deta) < constant::SMALL_ANGLE } { set_dxMin(find_dxMin()); } @@ -484,12 +486,25 @@ namespace metric { * @param x1 radial coordinate along the axis (code units) */ Inline auto polar_area(const real_t& x1) const -> real_t { - return dchi * math::exp(x1 * dchi + chi_min) * - (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a)) * - math::sqrt( - ONE + TWO * (r0 + math::exp(x1 * dchi + chi_min)) / - (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a))) * - (ONE - math::cos(eta2theta(HALF * deta))); + if constexpr (D != Dim::_1D) { + if (small_angle) { + const real_t dtheta = eta2theta(HALF * deta); + return dchi * math::exp(x1 * dchi + chi_min) * + (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a)) * + math::sqrt( + ONE + TWO * (r0 + math::exp(x1 * dchi + chi_min)) / + (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a))) * + (static_cast(48) - SQR(dtheta)) * SQR(dtheta) / + static_cast(384); + } else { + return dchi * math::exp(x1 * dchi + chi_min) * + (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a)) * + math::sqrt( + ONE + TWO * (r0 + math::exp(x1 * dchi + chi_min)) / + (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a))) * + (ONE - math::cos(eta2theta(HALF * deta))); + } + } } /**