From c017cefbb8fe20d69a4f4207a2ecab4dfc2a03a4 Mon Sep 17 00:00:00 2001 From: AnujKankani Date: Tue, 23 Dec 2025 11:22:45 -0600 Subject: [PATCH 1/3] weights compatible with 1.4rc --- src/archetypes/spatial_dist.h | 27 +++++++++++++++++++-------- src/engines/grpic.hpp | 4 ++-- src/kernels/injectors.hpp | 32 ++++++++++++++++++-------------- 3 files changed, 39 insertions(+), 24 deletions(-) diff --git a/src/archetypes/spatial_dist.h b/src/archetypes/spatial_dist.h index 68477208..c8f3c8a3 100644 --- a/src/archetypes/spatial_dist.h +++ b/src/archetypes/spatial_dist.h @@ -41,8 +41,11 @@ namespace arch { struct Uniform : public SpatialDistribution { Uniform(const M& metric) : SpatialDistribution { metric } {} - Inline auto operator()(const coord_t&) const -> real_t { - return ONE; + Inline auto operator()(const coord_t&, + real_t& nppc_distribution, + real_t& weight_distribution) const { + nppc_distribution = ONE; + weight_distribution = ONE; } }; @@ -66,7 +69,9 @@ namespace arch { , target_density { target_density } , target_max_density { target_max_density } {} - Inline auto operator()(const coord_t& x_Ph) const -> real_t { + Inline auto operator()(const coord_t& x_Ph, + real_t& nppc_distribution, + real_t& weight_distribution) const { coord_t x_Cd { ZERO }; metric.template convert(x_Ph, x_Cd); real_t dens { ZERO }; @@ -86,9 +91,11 @@ namespace arch { } const auto target = target_density(x_Ph); if (0.9 * target > dens) { - return (target - dens) / target_max_density; + nppc_distribution = (target - dens) / target_max_density; + weight_distribution = ONE; } else { - return ZERO; + nppc_distribution = ZERO; + weight_distribution = ONE; } } }; @@ -110,7 +117,9 @@ namespace arch { , idx { idx } , target_density { target_density } {} - Inline auto operator()(const coord_t& x_Ph) const -> real_t { + Inline auto operator()(const coord_t& x_Ph, + real_t& nppc_distribution, + real_t& weight_distribution) const { coord_t x_Cd { ZERO }; metric.template convert(x_Ph, x_Cd); real_t dens { ZERO }; @@ -129,9 +138,11 @@ namespace arch { raise::KernelError(HERE, "Invalid dimension"); } if (0.9 * target_density > dens) { - return (target_density - dens) / target_density; + nppc_distribution = (target_density - dens) / target_density; + weight_distribution = ONE; } else { - return ZERO; + nppc_distribution = ZERO; + weight_distribution = ONE; } } }; diff --git a/src/engines/grpic.hpp b/src/engines/grpic.hpp index 35e544a9..4d9d89d9 100644 --- a/src/engines/grpic.hpp +++ b/src/engines/grpic.hpp @@ -94,9 +94,9 @@ namespace ntt { void step_forward(timer::Timers& timers, domain_t& dom) override { const auto fieldsolver_enabled = m_params.template get( - "algorithms.toggles.fieldsolver"); + "algorithms.fieldsolver.enable"); const auto deposit_enabled = m_params.template get( - "algorithms.toggles.deposit"); + "algorithms.deposit.enable"); const auto clear_interval = m_params.template get( "particles.clear_interval"); diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index cddb1883..c51f811f 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -615,8 +615,11 @@ namespace kernel { return idx_h(); } - Inline auto injected_ppc(const coord_t& x_Ph) const -> npart_t { - const auto ppc_real = ppc0 * spatial_dist(x_Ph); + Inline auto injected_ppc(const coord_t& x_Ph, + real_t& ppc_dist, + real_t& weight_dist) const -> npart_t { + spatial_dist(x_Ph, ppc_dist, weight_dist); + const auto ppc_real = ppc0 * ppc_dist; auto ppc = static_cast(ppc_real); auto rand_gen = random_pool.get_state(); if (Random(rand_gen) < (ppc_real - static_cast(ppc))) { @@ -682,15 +685,15 @@ namespace kernel { coord_t x_Cd { i1_ + HALF }; coord_t x_Ph { ZERO }; metric.template convert(x_Cd, x_Ph); - - const auto ppc = injected_ppc(x_Ph); + auto ppc_dist { ZERO }, weight_dist { ONE }; + const auto ppc = injected_ppc(x_Ph, ppc_dist, weight_dist); if (ppc == 0) { return; } - auto weight = ONE; + auto weight = ONE * weight_dist; if constexpr (M::CoordType != Coord::Cart) { - weight = metric.sqrt_det_h({ i1_ + HALF }) * inv_V0; + weight = metric.sqrt_det_h({ i1_ + HALF }) * inv_V0 * weight_dist; } for (auto p { 0u }; p < ppc; ++p) { const auto index = Kokkos::atomic_fetch_add(&idx(), 1); @@ -732,15 +735,16 @@ namespace kernel { x_Cd_[2] = ZERO; } metric.template convert(x_Cd, x_Ph); - - const auto ppc = injected_ppc(x_Ph); + + auto ppc_dist { ZERO }, weight_dist { ONE }; + const auto ppc = injected_ppc(x_Ph, ppc_dist, weight_dist); if (ppc == 0) { return; } - auto weight = ONE; + auto weight = ONE * weight_dist; if constexpr (M::CoordType != Coord::Cart) { - weight = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF }) * inv_V0; + weight = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF }) * inv_V0 * weight_dist; } for (auto p { 0u }; p < ppc; ++p) { const auto index = Kokkos::atomic_fetch_add(&idx(), 1); @@ -796,16 +800,16 @@ namespace kernel { coord_t x_Cd { i1_ + HALF, i2_ + HALF, i3_ + HALF }; coord_t x_Ph { ZERO }; metric.template convert(x_Cd, x_Ph); - - const auto ppc = injected_ppc(x_Ph); + auto ppc_dist { ZERO }, weight_dist { ONE }; + const auto ppc = injected_ppc(x_Ph, ppc_dist, weight_dist); if (ppc == 0) { return; } - auto weight = ONE; + auto weight = ONE * weight_dist; if constexpr (M::CoordType != Coord::Cart) { weight = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF, i3_ + HALF }) * - inv_V0; + inv_V0 * weight_dist; } for (auto p { 0u }; p < ppc; ++p) { const auto index = Kokkos::atomic_fetch_add(&idx(), 1); From 17471c59e125775fb2a9ef147417d62e77ac37c2 Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Wed, 14 Jan 2026 16:09:01 -0500 Subject: [PATCH 2/3] small-angle approximation for polar area in GR, replicates SR --- src/metrics/kerr_schild.h | 20 +++++++++++++++----- src/metrics/qkerr_schild.h | 28 +++++++++++++++++++++------- 2 files changed, 36 insertions(+), 12 deletions(-) diff --git a/src/metrics/kerr_schild.h b/src/metrics/kerr_schild.h index 14441f99..45a8d1cf 100644 --- a/src/metrics/kerr_schild.h +++ b/src/metrics/kerr_schild.h @@ -85,7 +85,8 @@ namespace metric { , dphi { (x3_max - x3_min) / nx3 } , dr_inv { ONE / dr } , dtheta_inv { ONE / dtheta } - , dphi_inv { ONE / dphi } { + , dphi_inv { ONE / dphi } + , small_angle { HALF * dtheta < constant::SMALL_ANGLE } { set_dxMin(find_dxMin()); } @@ -443,12 +444,21 @@ namespace metric { /** * differential area at the pole (used in axisymmetric solvers) * @param x1 radial coordinate along the axis (code units) + * @note uses small-angle approximation when the resolution is too high */ Inline auto polar_area(real_t x1) const -> real_t { - return dr * (SQR(x1 * dr + x1_min) + SQR(a)) * - math::sqrt(ONE + TWO * (x1 * dr + x1_min) / - (SQR(x1 * dr + x1_min) + SQR(a))) * - (ONE - math::cos(HALF * dtheta)); + if (small_angle) { + return dr * (SQR(x1 * dr + x1_min) + SQR(a)) * + math::sqrt(ONE + TWO * (x1 * dr + x1_min) / + (SQR(x1 * dr + x1_min) + SQR(a))) * + (static_cast(48) - SQR(dtheta)) * SQR(dtheta) / + static_cast(384); + } else { + return dr * (SQR(x1 * dr + x1_min) + SQR(a)) * + math::sqrt(ONE + TWO * (x1 * dr + x1_min) / + (SQR(x1 * dr + x1_min) + SQR(a))) * + (ONE - math::cos(HALF * dtheta)); + } } /** diff --git a/src/metrics/qkerr_schild.h b/src/metrics/qkerr_schild.h index 9b608f35..90f366ef 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(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()); } @@ -503,15 +505,27 @@ namespace metric { * differential area at the pole (used in axisymmetric solvers) * @note approximate solution for the polar area * @param x1 radial coordinate along the axis (code units) + * @note uses small-angle approximation when the resolution is too high */ Inline auto polar_area(real_t x1) const -> real_t { if constexpr (D != Dim::_1D) { - 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 (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))); + } } } From 690385d5c9317dc6455a887c683226b9c55aa42e Mon Sep 17 00:00:00 2001 From: Alisa Galishnikova <55898700+alisagk@users.noreply.github.com> Date: Wed, 14 Jan 2026 16:14:42 -0500 Subject: [PATCH 3/3] added definition of small_angle to ks metric --- src/metrics/kerr_schild.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/metrics/kerr_schild.h b/src/metrics/kerr_schild.h index 45a8d1cf..4df88cf4 100644 --- a/src/metrics/kerr_schild.h +++ b/src/metrics/kerr_schild.h @@ -40,6 +40,7 @@ namespace metric { const real_t dr, dtheta, dphi; const real_t dr_inv, dtheta_inv, dphi_inv; + const bool small_angle; Inline auto Delta(real_t r) const -> real_t { return SQR(r) - TWO * r + SQR(a);