From 1c90d79d56d1d4077fa809c474e40463f72a7f21 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 28 Jul 2026 23:22:50 +0000 Subject: [PATCH 1/2] extended external force to be able to access particle properties --- examples/external_fields/pgen.hpp | 3 ++- src/global/traits/archetypes.h | 27 ++++++++++++++++++++------- src/kernels/pushers/sr.hpp | 9 +++++---- tests/archetypes/pgen.cpp | 3 ++- tests/global/traits_archetypes.cpp | 6 +++--- tests/global/traits_policies.cpp | 2 +- tests/kernels/ext_force.cpp | 9 ++++++--- 7 files changed, 39 insertions(+), 20 deletions(-) diff --git a/examples/external_fields/pgen.hpp b/examples/external_fields/pgen.hpp index ac4f6faed..30fb2473d 100644 --- a/examples/external_fields/pgen.hpp +++ b/examples/external_fields/pgen.hpp @@ -40,7 +40,8 @@ namespace user { */ // f_ext: external force-field (acceleration): - Inline auto fx1(const coord_t&) const -> real_t { + Inline auto fx1(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const + -> real_t { return (sp % 2u == 0u) ? -HALF : HALF; } diff --git a/src/global/traits/archetypes.h b/src/global/traits/archetypes.h index e170170fb..c5582e74b 100644 --- a/src/global/traits/archetypes.h +++ b/src/global/traits/archetypes.h @@ -4,7 +4,7 @@ * @implements * - EnrgDistClass<> - checks if a class can be used as an energy distribution * - SpatialDistClass<> - checks if a class can be used as a spatial distribution - * - traits::fieldsetter::HasFx1, ::HasFx2, ::HasFx3 - checks for F functions in field setter class + * - traits::fieldsetter::HasFx1, ::HasFx2, ::HasFx3 - checks for particle-aware F functions fx*(x, particles, p) * - traits::fieldsetter::HasEx1, ::HasEx2, ::HasEx3 - checks for E functions in field setter class * - traits::fieldsetter::HasBx1, ::HasBx2, ::HasBx3 - checks for B functions in field setter class * - traits::fieldsetter::HasDx1, ::HasDx2, ::HasDx3 - checks for D functions in field setter class @@ -26,6 +26,10 @@ #include +namespace ntt { + struct ParticleArrays; +} // namespace ntt + template concept EnrgDistClass = requires(const ED& edist, const coord_t& x_Ph, @@ -50,18 +54,27 @@ concept SpatialDistClass = SimpleSpatialDistClass or namespace traits::fieldsetter { template - concept HasFx1 = requires(const T& t, const coord_t& x_Ph) { - { t.fx1(x_Ph) } -> std::convertible_to; + concept HasFx1 = requires(const T& t, + const coord_t& x_Ph, + const ntt::ParticleArrays& prtls, + prtlidx_t p) { + { t.fx1(x_Ph, prtls, p) } -> std::convertible_to; }; template - concept HasFx2 = requires(const T& t, const coord_t& x_Ph) { - { t.fx2(x_Ph) } -> std::convertible_to; + concept HasFx2 = requires(const T& t, + const coord_t& x_Ph, + const ntt::ParticleArrays& prtls, + prtlidx_t p) { + { t.fx2(x_Ph, prtls, p) } -> std::convertible_to; }; template - concept HasFx3 = requires(const T& t, const coord_t& x_Ph) { - { t.fx3(x_Ph) } -> std::convertible_to; + concept HasFx3 = requires(const T& t, + const coord_t& x_Ph, + const ntt::ParticleArrays& prtls, + prtlidx_t p) { + { t.fx3(x_Ph, prtls, p) } -> std::convertible_to; }; template diff --git a/src/kernels/pushers/sr.hpp b/src/kernels/pushers/sr.hpp index 78b5c4916..bc7d79f08 100644 --- a/src/kernels/pushers/sr.hpp +++ b/src/kernels/pushers/sr.hpp @@ -245,7 +245,7 @@ namespace kernel::sr { // compute the external force either user-provided or from the atmosphere model if constexpr (HasExtForce or Atm) { - getExternalForce(xp_Cd, xp_Ph, external_force_Cart); + getExternalForce(xp_Cd, xp_Ph, p, external_force_Cart); } if (ctx.pusher_flags & ParticlePusher::GCA) { @@ -1425,18 +1425,19 @@ namespace kernel::sr { Inline void getExternalForce(const coord_t& xp_Cd, const coord_t& xp_Ph, + [[maybe_unused]] prtlidx_t p, vec_t& external_force_Cart) const requires(Atm or HasExtForce) { real_t f_x1 = ZERO, f_x2 = ZERO, f_x3 = ZERO; if constexpr (HasExtFx1) { - f_x1 = policies.external_fields_policy.fx1(xp_Ph); + f_x1 = policies.external_fields_policy.fx1(xp_Ph, particles, p); } if constexpr (HasExtFx2) { - f_x2 = policies.external_fields_policy.fx2(xp_Ph); + f_x2 = policies.external_fields_policy.fx2(xp_Ph, particles, p); } if constexpr (HasExtFx3) { - f_x3 = policies.external_fields_policy.fx3(xp_Ph); + f_x3 = policies.external_fields_policy.fx3(xp_Ph, particles, p); } if constexpr (Atm) { if constexpr (D == Dim::_1D or D == Dim::_2D or D == Dim::_3D) { diff --git a/tests/archetypes/pgen.cpp b/tests/archetypes/pgen.cpp index d95e9cf9c..fb428493b 100644 --- a/tests/archetypes/pgen.cpp +++ b/tests/archetypes/pgen.cpp @@ -25,7 +25,8 @@ struct CustomFieldsetter { template struct ExtForce { - Inline auto fx1(const coord_t&) const -> real_t { + Inline auto fx1(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const + -> real_t { return ZERO; } diff --git a/tests/global/traits_archetypes.cpp b/tests/global/traits_archetypes.cpp index d72545102..575b7abe5 100644 --- a/tests/global/traits_archetypes.cpp +++ b/tests/global/traits_archetypes.cpp @@ -89,15 +89,15 @@ struct WithDx1Dx2Dx3 { }; struct WithFx1Fx2Fx3 { - real_t fx1(const coord_t&) const { + real_t fx1(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const { return ZERO; } - real_t fx2(const coord_t&) const { + real_t fx2(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const { return ZERO; } - real_t fx3(const coord_t&) const { + real_t fx3(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const { return ZERO; } }; diff --git a/tests/global/traits_policies.cpp b/tests/global/traits_policies.cpp index 5f5a83f9f..82ceb7e12 100644 --- a/tests/global/traits_policies.cpp +++ b/tests/global/traits_policies.cpp @@ -105,7 +105,7 @@ static_assert(not EmissionPolicyClass); // --- ExtFieldsPolicyClass --- struct WithFx1 { - real_t fx1(const coord_t&) const { + real_t fx1(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const { return ZERO; } }; diff --git a/tests/kernels/ext_force.cpp b/tests/kernels/ext_force.cpp index 69d773d10..d9e8203a5 100644 --- a/tests/kernels/ext_force.cpp +++ b/tests/kernels/ext_force.cpp @@ -54,15 +54,18 @@ void put_value(array_t& arr, T v, prtlidx_t p) { struct Force { Force(real_t force) : force { force } {} - Inline auto fx1(const coord_t&) const -> real_t { + Inline auto fx1(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const + -> real_t { return force * math::sin(ONE) * math::sin(ONE); } - Inline auto fx2(const coord_t&) const -> real_t { + Inline auto fx2(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const + -> real_t { return force * math::sin(ONE) * math::cos(ONE); } - Inline auto fx3(const coord_t&) const -> real_t { + Inline auto fx3(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const + -> real_t { return force * math::cos(ONE); } From 87bf6ad608d0550a625b630f866f95cb39106b42 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 16 Sep 2026 20:27:38 +0000 Subject: [PATCH 2/2] refactor particle-dependent external force via `HasExtFx*WithIndex` to keep backwards compatible `HasExtFx*` --- examples/external_fields/pgen.hpp | 3 +- src/global/traits/archetypes.h | 46 +++++++++------- src/global/traits/policies.h | 11 ++-- src/kernels/pushers/sr.hpp | 36 +++++++++---- tests/archetypes/pgen.cpp | 3 +- tests/global/traits_archetypes.cpp | 34 ++++++++++-- tests/global/traits_policies.cpp | 9 +++- tests/kernels/ext_force.cpp | 84 +++++++++++++++++++++++------- 8 files changed, 163 insertions(+), 63 deletions(-) diff --git a/examples/external_fields/pgen.hpp b/examples/external_fields/pgen.hpp index 30fb2473d..ac4f6faed 100644 --- a/examples/external_fields/pgen.hpp +++ b/examples/external_fields/pgen.hpp @@ -40,8 +40,7 @@ namespace user { */ // f_ext: external force-field (acceleration): - Inline auto fx1(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const - -> real_t { + Inline auto fx1(const coord_t&) const -> real_t { return (sp % 2u == 0u) ? -HALF : HALF; } diff --git a/src/global/traits/archetypes.h b/src/global/traits/archetypes.h index c5582e74b..0dfced45c 100644 --- a/src/global/traits/archetypes.h +++ b/src/global/traits/archetypes.h @@ -4,7 +4,9 @@ * @implements * - EnrgDistClass<> - checks if a class can be used as an energy distribution * - SpatialDistClass<> - checks if a class can be used as a spatial distribution - * - traits::fieldsetter::HasFx1, ::HasFx2, ::HasFx3 - checks for particle-aware F functions fx*(x, particles, p) + * - traits::fieldsetter::HasFx1, ::HasFx2, ::HasFx3 - checks for F functions in field setter class + * - traits::fieldsetter::HasFx1WithIndex, ::HasFx2WithIndex, ::HasFx3WithIndex + * - checks for F functions taking an additional particle index * - traits::fieldsetter::HasEx1, ::HasEx2, ::HasEx3 - checks for E functions in field setter class * - traits::fieldsetter::HasBx1, ::HasBx2, ::HasBx3 - checks for B functions in field setter class * - traits::fieldsetter::HasDx1, ::HasDx2, ::HasDx3 - checks for D functions in field setter class @@ -26,10 +28,6 @@ #include -namespace ntt { - struct ParticleArrays; -} // namespace ntt - template concept EnrgDistClass = requires(const ED& edist, const coord_t& x_Ph, @@ -54,27 +52,35 @@ concept SpatialDistClass = SimpleSpatialDistClass or namespace traits::fieldsetter { template - concept HasFx1 = requires(const T& t, - const coord_t& x_Ph, - const ntt::ParticleArrays& prtls, - prtlidx_t p) { - { t.fx1(x_Ph, prtls, p) } -> std::convertible_to; + concept HasFx1 = requires(const T& t, const coord_t& x_Ph) { + { t.fx1(x_Ph) } -> std::convertible_to; + }; + + template + concept HasFx2 = requires(const T& t, const coord_t& x_Ph) { + { t.fx2(x_Ph) } -> std::convertible_to; + }; + + template + concept HasFx3 = requires(const T& t, const coord_t& x_Ph) { + { t.fx3(x_Ph) } -> std::convertible_to; + }; + + // force components which additionally depend on the particle index + // (e.g., to access per-particle properties captured by the setter itself) + template + concept HasFx1WithIndex = requires(const T& t, const coord_t& x_Ph, prtlidx_t p) { + { t.fx1(x_Ph, p) } -> std::convertible_to; }; template - concept HasFx2 = requires(const T& t, - const coord_t& x_Ph, - const ntt::ParticleArrays& prtls, - prtlidx_t p) { - { t.fx2(x_Ph, prtls, p) } -> std::convertible_to; + concept HasFx2WithIndex = requires(const T& t, const coord_t& x_Ph, prtlidx_t p) { + { t.fx2(x_Ph, p) } -> std::convertible_to; }; template - concept HasFx3 = requires(const T& t, - const coord_t& x_Ph, - const ntt::ParticleArrays& prtls, - prtlidx_t p) { - { t.fx3(x_Ph, prtls, p) } -> std::convertible_to; + concept HasFx3WithIndex = requires(const T& t, const coord_t& x_Ph, prtlidx_t p) { + { t.fx3(x_Ph, p) } -> std::convertible_to; }; template diff --git a/src/global/traits/policies.h b/src/global/traits/policies.h index 7f6d98985..a317fd687 100644 --- a/src/global/traits/policies.h +++ b/src/global/traits/policies.h @@ -109,10 +109,13 @@ namespace traits::extfields { template concept ExtFieldsPolicyClass = (::traits::fieldsetter::HasFx1 or ::traits::fieldsetter::HasFx2 or - ::traits::fieldsetter::HasFx3 or ::traits::fieldsetter::HasEx1 or - ::traits::fieldsetter::HasEx2 or ::traits::fieldsetter::HasEx3 or - ::traits::fieldsetter::HasBx1 or ::traits::fieldsetter::HasBx2 or - ::traits::fieldsetter::HasBx3) or + ::traits::fieldsetter::HasFx3 or + ::traits::fieldsetter::HasFx1WithIndex or + ::traits::fieldsetter::HasFx2WithIndex or + ::traits::fieldsetter::HasFx3WithIndex or + ::traits::fieldsetter::HasEx1 or ::traits::fieldsetter::HasEx2 or + ::traits::fieldsetter::HasEx3 or ::traits::fieldsetter::HasBx1 or + ::traits::fieldsetter::HasBx2 or ::traits::fieldsetter::HasBx3) or ::traits::extfields::IsNoPolicy; namespace traits::custom_prtl_update { diff --git a/src/kernels/pushers/sr.hpp b/src/kernels/pushers/sr.hpp index bc7d79f08..6d691d3df 100644 --- a/src/kernels/pushers/sr.hpp +++ b/src/kernels/pushers/sr.hpp @@ -67,11 +67,19 @@ namespace kernel::sr { using F = typename P::ExternalFieldsPolicy; static constexpr auto Atm = P::ApplyAtmosphere; - static constexpr auto D = M::Dim; - static constexpr auto HasExtFx1 = ::traits::fieldsetter::HasFx1; - static constexpr auto HasExtFx2 = ::traits::fieldsetter::HasFx2; - static constexpr auto HasExtFx3 = ::traits::fieldsetter::HasFx3; - static constexpr auto HasExtForce = HasExtFx1 or HasExtFx2 or HasExtFx3; + static constexpr auto D = M::Dim; + static constexpr auto HasExtFx1 = ::traits::fieldsetter::HasFx1; + static constexpr auto HasExtFx2 = ::traits::fieldsetter::HasFx2; + static constexpr auto HasExtFx3 = ::traits::fieldsetter::HasFx3; + static constexpr auto HasExtFx1WithIndex = + ::traits::fieldsetter::HasFx1WithIndex; + static constexpr auto HasExtFx2WithIndex = + ::traits::fieldsetter::HasFx2WithIndex; + static constexpr auto HasExtFx3WithIndex = + ::traits::fieldsetter::HasFx3WithIndex; + static constexpr auto HasExtForce = HasExtFx1 or HasExtFx2 or HasExtFx3 or + HasExtFx1WithIndex or + HasExtFx2WithIndex or HasExtFx3WithIndex; static constexpr auto HasExtEx1 = ::traits::fieldsetter::HasEx1; static constexpr auto HasExtEx2 = ::traits::fieldsetter::HasEx2; static constexpr auto HasExtEx3 = ::traits::fieldsetter::HasEx3; @@ -1430,14 +1438,20 @@ namespace kernel::sr { requires(Atm or HasExtForce) { real_t f_x1 = ZERO, f_x2 = ZERO, f_x3 = ZERO; - if constexpr (HasExtFx1) { - f_x1 = policies.external_fields_policy.fx1(xp_Ph, particles, p); + if constexpr (HasExtFx1WithIndex) { + f_x1 = policies.external_fields_policy.fx1(xp_Ph, p); + } else if constexpr (HasExtFx1) { + f_x1 = policies.external_fields_policy.fx1(xp_Ph); } - if constexpr (HasExtFx2) { - f_x2 = policies.external_fields_policy.fx2(xp_Ph, particles, p); + if constexpr (HasExtFx2WithIndex) { + f_x2 = policies.external_fields_policy.fx2(xp_Ph, p); + } else if constexpr (HasExtFx2) { + f_x2 = policies.external_fields_policy.fx2(xp_Ph); } - if constexpr (HasExtFx3) { - f_x3 = policies.external_fields_policy.fx3(xp_Ph, particles, p); + if constexpr (HasExtFx3WithIndex) { + f_x3 = policies.external_fields_policy.fx3(xp_Ph, p); + } else if constexpr (HasExtFx3) { + f_x3 = policies.external_fields_policy.fx3(xp_Ph); } if constexpr (Atm) { if constexpr (D == Dim::_1D or D == Dim::_2D or D == Dim::_3D) { diff --git a/tests/archetypes/pgen.cpp b/tests/archetypes/pgen.cpp index fb428493b..d95e9cf9c 100644 --- a/tests/archetypes/pgen.cpp +++ b/tests/archetypes/pgen.cpp @@ -25,8 +25,7 @@ struct CustomFieldsetter { template struct ExtForce { - Inline auto fx1(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const - -> real_t { + Inline auto fx1(const coord_t&) const -> real_t { return ZERO; } diff --git a/tests/global/traits_archetypes.cpp b/tests/global/traits_archetypes.cpp index 575b7abe5..70ce3e9e2 100644 --- a/tests/global/traits_archetypes.cpp +++ b/tests/global/traits_archetypes.cpp @@ -89,15 +89,29 @@ struct WithDx1Dx2Dx3 { }; struct WithFx1Fx2Fx3 { - real_t fx1(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const { + real_t fx1(const coord_t&) const { return ZERO; } - real_t fx2(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const { + real_t fx2(const coord_t&) const { return ZERO; } - real_t fx3(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const { + real_t fx3(const coord_t&) const { + return ZERO; + } +}; + +struct WithIndexedFx1Fx2Fx3 { + real_t fx1(const coord_t&, prtlidx_t) const { + return ZERO; + } + + real_t fx2(const coord_t&, prtlidx_t) const { + return ZERO; + } + + real_t fx3(const coord_t&, prtlidx_t) const { return ZERO; } }; @@ -135,6 +149,20 @@ static_assert(traits::fieldsetter::HasFx1); static_assert(traits::fieldsetter::HasFx2); static_assert(traits::fieldsetter::HasFx3); +// the indexed variants are a separate, opt-in signature: neither form +// satisfies the other's trait +static_assert( + traits::fieldsetter::HasFx1WithIndex); +static_assert( + traits::fieldsetter::HasFx2WithIndex); +static_assert( + traits::fieldsetter::HasFx3WithIndex); +static_assert( + not traits::fieldsetter::HasFx1WithIndex); +static_assert( + not traits::fieldsetter::HasFx1); +static_assert(not traits::fieldsetter::HasFx1WithIndex); + // conditional variants require 3-arg signature returning Kokkos::pair static_assert( traits::fieldsetter::HasConditionalEx1); diff --git a/tests/global/traits_policies.cpp b/tests/global/traits_policies.cpp index 82ceb7e12..7f212dab3 100644 --- a/tests/global/traits_policies.cpp +++ b/tests/global/traits_policies.cpp @@ -105,7 +105,13 @@ static_assert(not EmissionPolicyClass); // --- ExtFieldsPolicyClass --- struct WithFx1 { - real_t fx1(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const { + real_t fx1(const coord_t&) const { + return ZERO; + } +}; + +struct WithIndexedFx1 { + real_t fx1(const coord_t&, prtlidx_t) const { return ZERO; } }; @@ -113,6 +119,7 @@ struct WithFx1 { struct Empty {}; static_assert(ExtFieldsPolicyClass); +static_assert(ExtFieldsPolicyClass); static_assert(not ExtFieldsPolicyClass); // --- CustomParticleUpdatePolicyClass with a real updater --- diff --git a/tests/kernels/ext_force.cpp b/tests/kernels/ext_force.cpp index d9e8203a5..0a705f19b 100644 --- a/tests/kernels/ext_force.cpp +++ b/tests/kernels/ext_force.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include using namespace ntt; @@ -51,21 +52,25 @@ void put_value(array_t& arr, T v, prtlidx_t p) { Kokkos::deep_copy(arr, h); } +// per-particle weight given to the two test particles; the indexed force +// setter reads it back from the particle arrays and scales the force by it +auto weight_of(prtlidx_t p) -> real_t { + return (p == 0) ? ONE : TWO; +} + +// force setter depending only on the position struct Force { Force(real_t force) : force { force } {} - Inline auto fx1(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const - -> real_t { + Inline auto fx1(const coord_t&) const -> real_t { return force * math::sin(ONE) * math::sin(ONE); } - Inline auto fx2(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const - -> real_t { + Inline auto fx2(const coord_t&) const -> real_t { return force * math::sin(ONE) * math::cos(ONE); } - Inline auto fx3(const coord_t&, const ntt::ParticleArrays&, prtlidx_t) const - -> real_t { + Inline auto fx3(const coord_t&) const -> real_t { return force * math::cos(ONE); } @@ -73,7 +78,31 @@ struct Force { const real_t force; }; -template +// force setter carrying its own copy of the particle arrays and reading +// per-particle properties through the particle index +struct IndexedForce { + IndexedForce(real_t force, const ParticleArrays& prtls) + : force { force } + , prtls { prtls } {} + + Inline auto fx1(const coord_t&, prtlidx_t p) const -> real_t { + return prtls.weight(p) * force * math::sin(ONE) * math::sin(ONE); + } + + Inline auto fx2(const coord_t&, prtlidx_t p) const -> real_t { + return prtls.weight(p) * force * math::sin(ONE) * math::cos(ONE); + } + + Inline auto fx3(const coord_t&, prtlidx_t p) const -> real_t { + return prtls.weight(p) * force * math::cos(ONE); + } + +private: + const real_t force; + const ParticleArrays prtls; +}; + +template void testPusher(const std::vector& res) { static_assert(M::Dim == 3); raise::ErrorIf(res.size() != M::Dim, "res.size() != M::Dim", HERE); @@ -144,6 +173,7 @@ void testPusher(const std::vector& res) { put_value(ux1, ux1_0, 0); put_value(ux2, ux2_0, 0); put_value(ux3, ux3_0, 0); + put_value(weight, weight_of(0), 0); put_value(tag, ParticleTag::alive, 0); put_value(i1, (int)(x1_0), 1); @@ -155,12 +185,11 @@ void testPusher(const std::vector& res) { put_value(ux1, -ux1_0, 1); put_value(ux2, -ux2_0, 1); put_value(ux3, -ux3_0, 1); + put_value(weight, weight_of(1), 1); put_value(tag, ParticleTag::alive, 1); const real_t eps = std::is_same_v ? 1e-4 : 1e-6; - const auto ext_force = Force { f_mag }; - static plog::RollingFileAppender file_appender( "pusher_log.csv"); plog::init(plog::verbose, &file_appender); @@ -189,8 +218,18 @@ void testPusher(const std::vector& res) { pusher_arrays.ux2 = ux2; pusher_arrays.ux3 = ux3; pusher_arrays.phi = phi; + pusher_arrays.weight = weight; pusher_arrays.tag = tag; + using ext_force_t = std::conditional_t; + const auto ext_force = [&]() -> ext_force_t { + if constexpr (Indexed) { + return { f_mag, pusher_arrays }; + } else { + return { f_mag }; + } + }(); + const auto pusher_policy = ::kernel::sr::PusherPolicy, ::traits::emission::NoPolicy_t, @@ -256,12 +295,16 @@ void testPusher(const std::vector& res) { ux2_(1), ux3_(1)); + // the indexed setter scales the force by the particle weight + const real_t scale_0 = Indexed ? weight_of(0) : ONE; + const real_t scale_1 = Indexed ? weight_of(1) : ONE; + { - const real_t ux1_expect = ux1_0 + (time + dt) * f_mag * std::sin(ONE) * - std::sin(ONE); - const real_t ux2_expect = ux2_0 + (time + dt) * f_mag * std::sin(ONE) * - std::cos(ONE); - const real_t ux3_expect = ux3_0 + (time + dt) * f_mag * std::cos(ONE); + const real_t f = (time + dt) * f_mag * scale_0; + + const real_t ux1_expect = ux1_0 + f * std::sin(ONE) * std::sin(ONE); + const real_t ux2_expect = ux2_0 + f * std::sin(ONE) * std::cos(ONE); + const real_t ux3_expect = ux3_0 + f * std::cos(ONE); check_value(t, ux1_(0), ux1_expect, eps, "Particle #1 ux1"); check_value(t, ux2_(0), ux2_expect, eps, "Particle #1 ux2"); @@ -269,11 +312,11 @@ void testPusher(const std::vector& res) { } { - const real_t ux1_expect = -ux1_0 + (time + dt) * f_mag * std::sin(ONE) * - std::sin(ONE); - const real_t ux2_expect = -ux2_0 + (time + dt) * f_mag * std::sin(ONE) * - std::cos(ONE); - const real_t ux3_expect = -ux3_0 + (time + dt) * f_mag * std::cos(ONE); + const real_t f = (time + dt) * f_mag * scale_1; + + const real_t ux1_expect = -ux1_0 + f * std::sin(ONE) * std::sin(ONE); + const real_t ux2_expect = -ux2_0 + f * std::sin(ONE) * std::cos(ONE); + const real_t ux3_expect = -ux3_0 + f * std::cos(ONE); check_value(t, ux1_(1), ux1_expect, eps, "Particle #2 ux1"); check_value(t, ux2_(1), ux2_expect, eps, "Particle #2 ux2"); @@ -288,7 +331,8 @@ auto main(int argc, char* argv[]) -> int { try { using namespace ntt; - testPusher>({ 10, 10, 10 }); + testPusher, false>({ 10, 10, 10 }); + testPusher, true>({ 10, 10, 10 }); } catch (std::exception& e) { std::cerr << e.what() << '\n';