From 57f706ab2dd8da816fddf4ca23a02d4c060ecf9b Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Sat, 21 Mar 2026 01:00:16 +0000 Subject: [PATCH 01/14] added option to define custom particle BCs in pgen --- src/archetypes/traits.h | 8 ++++++ src/engines/srpic/particle_pusher.h | 30 ++++++++++++++++------ src/kernels/digital_filter.hpp | 4 +-- src/kernels/particle_pusher_sr.hpp | 40 ++++++++++++++++++++++++++--- 4 files changed, 68 insertions(+), 14 deletions(-) diff --git a/src/archetypes/traits.h b/src/archetypes/traits.h index 5d11ba3ab..d95ecf355 100644 --- a/src/archetypes/traits.h +++ b/src/archetypes/traits.h @@ -85,6 +85,14 @@ namespace arch { pgen.EmissionPolicy(time, sp, domain); }; + template + concept HasCustomPrtlBC = requires(const PG& pgen, + simtime_t time, + spidx_t sp, + D& domain) { + pgen.CustomParticleBoundary(time, sp, domain); + }; + template concept HasInitPrtls = requires(PG& pgen, D& domain) { { pgen.InitPrtls(domain) } -> std::same_as; diff --git a/src/engines/srpic/particle_pusher.h b/src/engines/srpic/particle_pusher.h index 2a0848de4..b144050f5 100644 --- a/src/engines/srpic/particle_pusher.h +++ b/src/engines/srpic/particle_pusher.h @@ -36,18 +36,29 @@ namespace ntt { const M& metric, const PG& pgen, const F& external_fields) { + auto get_custom_prtl_bc = [&]() { + if constexpr (arch::traits::pgen::HasCustomPrtlBC) { + return pgen.CustomParticleBoundary( + pusher_params.time, pusher_params.species_index, domain); + } else { + return kernel::sr::NoCustomPrtlBC_t {}; + } + }; + const auto custom_prtl_bc = get_custom_prtl_bc(); + if (emission_policy_flag == EmissionType::NONE) { const auto no_emission = kernel::NoEmissionPolicy_t {}; Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, - no_emission)); + no_emission, + custom_prtl_bc)); } else if (emission_policy_flag == EmissionType::SYNCHROTRON) { const auto photon_species = params.get( "radiation.emission.synchrotron.photon_species"); @@ -77,13 +88,14 @@ namespace ntt { Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, - emission_policy)); + emission_policy, + custom_prtl_bc)); const auto n_inj = emission_policy.numbers_injected(); raise::ErrorIf(n_inj.size() != 1, "Synchrotron emission should only inject one species", @@ -120,13 +132,14 @@ namespace ntt { Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, - emission_policy)); + emission_policy, + custom_prtl_bc)); const auto n_inj = emission_policy.numbers_injected(); raise::ErrorIf(n_inj.size() != 1, "Compton emission should only inject one species", @@ -148,13 +161,14 @@ namespace ntt { Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, - emission_policy)); + emission_policy, + custom_prtl_bc)); const auto emitted_species = emission_policy.emitted_species_indices(); const auto n_inj = emission_policy.numbers_injected(); raise::ErrorIf(emitted_species.size() != n_inj.size(), diff --git a/src/kernels/digital_filter.hpp b/src/kernels/digital_filter.hpp index 5ac60327d..7f3e52f1f 100644 --- a/src/kernels/digital_filter.hpp +++ b/src/kernels/digital_filter.hpp @@ -93,8 +93,8 @@ namespace kernel { ? (boundaries[2].second == FldsBC::CONDUCTOR) : false } , i1_max { size_[0] + N_GHOSTS } - , i2_max { (short)D > 1 ? (size_[1] + N_GHOSTS) : 0 } - , i3_max { (short)D > 2 ? (size_[2] + N_GHOSTS) : 0 } {} + , i2_max { (short)D > 1 ? (size_[(short)D > 1 ? 1 : 0] + N_GHOSTS) : 0 } + , i3_max { (short)D > 2 ? (size_[(short)D > 2 ? 2 : 0] + N_GHOSTS) : 0 } {} Inline void operator()(index_t i1) const { if constexpr ((D == Dim::_1D) && (C == Coord::Cart)) { diff --git a/src/kernels/particle_pusher_sr.hpp b/src/kernels/particle_pusher_sr.hpp index b8ad572de..b13e9ee7e 100644 --- a/src/kernels/particle_pusher_sr.hpp +++ b/src/kernels/particle_pusher_sr.hpp @@ -101,11 +101,17 @@ namespace kernel::sr { array_t tag; }; + template + struct NoCustomPrtlBC_t { + template + Inline void operator()(index_t, int, bool, const PusherKernel&) const {} + }; + /** * @tparam M Metric * @tparam F Additional force */ - template + template > requires metric::traits::HasD && metric::traits::HasTransformXYZ && metric::traits::HasConvertXYZ && metric::traits::HasTransform_i && metric::traits::HasConvert_i @@ -120,12 +126,12 @@ namespace kernel::sr { std::is_same>::value, "Invalid emission policy E for Pusher_kernel"); - private: const ParticlePusherFlags pusher_flags; const RadiativeDragFlags radiative_drag_flags; const randacc_ndfield_t EB; + public: const spidx_t sp; array_t i1, i2, i3; array_t i1_prev, i2_prev, i3_prev; @@ -136,14 +142,17 @@ namespace kernel::sr { array_t weight; array_t tag; const M metric; + const int ni1, ni2, ni3; + + private: const F external_fields; const E emission_policy; + const PBC custom_prtl_bc; const simtime_t time; const real_t coeff, dt; - const int ni1, ni2, ni3; bool is_absorb_i1min { false }, is_absorb_i1max { false }; bool is_absorb_i2min { false }, is_absorb_i2max { false }; bool is_absorb_i3min { false }, is_absorb_i3max { false }; @@ -153,6 +162,9 @@ namespace kernel::sr { bool is_reflect_i1min { false }, is_reflect_i1max { false }; bool is_reflect_i2min { false }, is_reflect_i2max { false }; bool is_reflect_i3min { false }, is_reflect_i3max { false }; + bool is_custom_i1min { false }, is_custom_i1max { false }; + bool is_custom_i2min { false }, is_custom_i2max { false }; + bool is_custom_i3min { false }, is_custom_i3max { false }; bool is_axis_i2min { false }, is_axis_i2max { false }; // gca parameters @@ -169,7 +181,8 @@ namespace kernel::sr { const randacc_ndfield_t& EB, const M& metric, const F& external_fields, - const E& emission_policy) + const E& emission_policy, + const PBC& custom_prtl_bc) : pusher_flags { pusher_params.pusher_flags } , radiative_drag_flags { pusher_params.radiative_drag_flags } , EB { EB } @@ -195,6 +208,7 @@ namespace kernel::sr { , metric { metric } , external_fields { external_fields } , emission_policy { emission_policy } + , custom_prtl_bc { custom_prtl_bc } , time { pusher_params.time } , coeff { HALF * (pusher_params.charge / pusher_params.mass) * pusher_params.omegaB0 * pusher_params.dt } @@ -256,6 +270,8 @@ namespace kernel::sr { is_periodic_i1max = (pusher_params.boundaries[0].second == PrtlBC::PERIODIC); is_reflect_i1min = (pusher_params.boundaries[0].first == PrtlBC::REFLECT); is_reflect_i1max = (pusher_params.boundaries[0].second == PrtlBC::REFLECT); + is_custom_i1min = (pusher_params.boundaries[0].first == PrtlBC::CUSTOM); + is_custom_i1max = (pusher_params.boundaries[0].second == PrtlBC::CUSTOM); if constexpr ((D == Dim::_2D) || (D == Dim::_3D)) { raise::ErrorIf(pusher_params.boundaries.size() < 2, "pusher_params.boundaries defined incorrectly", @@ -269,6 +285,8 @@ namespace kernel::sr { is_periodic_i2max = (pusher_params.boundaries[1].second == PrtlBC::PERIODIC); is_reflect_i2min = (pusher_params.boundaries[1].first == PrtlBC::REFLECT); is_reflect_i2max = (pusher_params.boundaries[1].second == PrtlBC::REFLECT); + is_custom_i2min = (pusher_params.boundaries[1].first == PrtlBC::CUSTOM); + is_custom_i2max = (pusher_params.boundaries[1].second == PrtlBC::CUSTOM); is_axis_i2min = (pusher_params.boundaries[1].first == PrtlBC::AXIS); is_axis_i2max = (pusher_params.boundaries[1].second == PrtlBC::AXIS); } @@ -285,6 +303,8 @@ namespace kernel::sr { is_periodic_i3max = (pusher_params.boundaries[2].second == PrtlBC::PERIODIC); is_reflect_i3min = (pusher_params.boundaries[2].first == PrtlBC::REFLECT); is_reflect_i3max = (pusher_params.boundaries[2].second == PrtlBC::REFLECT); + is_custom_i3min = (pusher_params.boundaries[2].first == PrtlBC::CUSTOM); + is_custom_i3max = (pusher_params.boundaries[2].second == PrtlBC::CUSTOM); } } @@ -1315,6 +1335,8 @@ namespace kernel::sr { i1(p) = 0; dx1(p) = ONE - dx1(p); invert_vel = true; + } else if (is_custom_i1min) { + custom_prtl_bc(p, 1, true, *this); } } else if (i1(p) >= ni1) { if (is_periodic_i1max) { @@ -1326,6 +1348,8 @@ namespace kernel::sr { i1(p) = ni1 - 1; dx1(p) = ONE - dx1(p); invert_vel = true; + } else if (is_custom_i1max) { + custom_prtl_bc(p, 1, false, *this); } } if (invert_vel) { @@ -1357,6 +1381,8 @@ namespace kernel::sr { i2(p) = 0; dx2(p) = ONE - dx2(p); invert_vel = true; + } else if (is_custom_i2min) { + custom_prtl_bc(p, 2, true, *this); } else if (is_axis_i2min) { i2(p) = 0; dx2(p) = ONE - dx2(p); @@ -1371,6 +1397,8 @@ namespace kernel::sr { i2(p) = ni2 - 1; dx2(p) = ONE - dx2(p); invert_vel = true; + } else if (is_custom_i2max) { + custom_prtl_bc(p, 2, false, *this); } else if (is_axis_i2max) { i2(p) = ni2 - 1; dx2(p) = ONE - dx2(p); @@ -1405,6 +1433,8 @@ namespace kernel::sr { i3(p) = 0; dx3(p) = ONE - dx3(p); invert_vel = true; + } else if (is_custom_i3min) { + custom_prtl_bc(p, 3, true, *this); } } else if (i3(p) >= ni3) { if (is_periodic_i3max) { @@ -1416,6 +1446,8 @@ namespace kernel::sr { i3(p) = ni3 - 1; dx3(p) = ONE - dx3(p); invert_vel = true; + } else if (is_custom_i3max) { + custom_prtl_bc(p, 3, false, *this); } } if (invert_vel) { From 64641af2d6875193e293993513795fba26af7980 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 31 Mar 2026 13:31:45 +0000 Subject: [PATCH 02/14] rewrite to custom particle update --- src/archetypes/traits.h | 4 ++-- src/engines/srpic/particle_pusher.h | 28 ++++++++++++++-------------- src/kernels/particle_pusher_sr.hpp | 28 +++++++++++++++------------- 3 files changed, 31 insertions(+), 29 deletions(-) diff --git a/src/archetypes/traits.h b/src/archetypes/traits.h index a7004fbff..51e4eafd6 100644 --- a/src/archetypes/traits.h +++ b/src/archetypes/traits.h @@ -86,11 +86,11 @@ namespace arch { }; template - concept HasCustomPrtlBC = requires(const PG& pgen, + concept HasCustomPrtlUpdate = requires(const PG& pgen, simtime_t time, spidx_t sp, D& domain) { - pgen.CustomParticleBoundary(time, sp, domain); + pgen.CustomParticleUpdate(time, sp, domain); }; template diff --git a/src/engines/srpic/particle_pusher.h b/src/engines/srpic/particle_pusher.h index d0849c625..abd5f3a59 100644 --- a/src/engines/srpic/particle_pusher.h +++ b/src/engines/srpic/particle_pusher.h @@ -26,7 +26,7 @@ namespace ntt { namespace srpic { template - void CallPusher(Domain& domain, + void CallPusher_WithExternalFieldFlag(Domain& domain, const SimulationParams& params, const kernel::sr::PusherParams& pusher_params, kernel::sr::PusherArrays& pusher_arrays, @@ -36,29 +36,29 @@ namespace ntt { const M& metric, const PG& pgen, const F& external_fields) { - auto get_custom_prtl_bc = [&]() { - if constexpr (arch::traits::pgen::HasCustomPrtlBC) { - return pgen.CustomParticleBoundary( + auto get_custom_prtl_update = [&]() { + if constexpr (arch::traits::pgen::HasCustomPrtlUpdate) { + return pgen.CustomParticleUpdate( pusher_params.time, pusher_params.species_index, domain); } else { - return kernel::sr::NoCustomPrtlBC_t {}; + return kernel::sr::NoCustomPrtlUpdate_t {}; } }; - const auto custom_prtl_bc = get_custom_prtl_bc(); + const auto custom_prtl_update = get_custom_prtl_update(); if (emission_policy_flag == EmissionType::NONE) { const auto no_emission = kernel::NoEmissionPolicy_t {}; Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, no_emission, - custom_prtl_bc)); + custom_prtl_update)); } else if (emission_policy_flag == EmissionType::SYNCHROTRON) { const auto photon_species = params.get( "radiation.emission.synchrotron.photon_species"); @@ -88,14 +88,14 @@ namespace ntt { Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, emission_policy, - custom_prtl_bc)); + custom_prtl_update)); const auto n_inj = emission_policy.numbers_injected(); raise::ErrorIf(n_inj.size() != 1, "Synchrotron emission should only inject one species", @@ -132,14 +132,14 @@ namespace ntt { Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, emission_policy, - custom_prtl_bc)); + custom_prtl_update)); const auto n_inj = emission_policy.numbers_injected(); raise::ErrorIf(n_inj.size() != 1, "Compton emission should only inject one species", @@ -160,14 +160,14 @@ namespace ntt { Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, emission_policy, - custom_prtl_bc)); + custom_prtl_update)); const auto emitted_species = emission_policy.emitted_species_indices(); const auto n_inj = emission_policy.numbers_injected(); raise::ErrorIf(emitted_species.size() != n_inj.size(), diff --git a/src/kernels/particle_pusher_sr.hpp b/src/kernels/particle_pusher_sr.hpp index 27b5d0a18..18f7659ce 100644 --- a/src/kernels/particle_pusher_sr.hpp +++ b/src/kernels/particle_pusher_sr.hpp @@ -102,16 +102,16 @@ namespace kernel::sr { }; template - struct NoCustomPrtlBC_t { + struct NoCustomPrtlUpdate_t { template - Inline void operator()(index_t, int, bool, const PusherKernel&) const {} + Inline void operator()(index_t, coord_t&, const PusherKernel&) const {} }; /** * @tparam M Metric * @tparam F Additional force */ - template > + template > requires metric::traits::HasD && metric::traits::HasTransformXYZ && metric::traits::HasConvertXYZ && metric::traits::HasTransform_i && metric::traits::HasConvert_i @@ -133,6 +133,8 @@ namespace kernel::sr { public: const spidx_t sp; + simtime_t time; + const real_t dt; array_t i1, i2, i3; array_t i1_prev, i2_prev, i3_prev; array_t dx1, dx2, dx3; @@ -148,9 +150,9 @@ namespace kernel::sr { const F external_fields; const E emission_policy; - const PBC custom_prtl_bc; + const PUPD custom_prtl_update; - const real_t coeff, dt; + const real_t coeff; bool is_absorb_i1min { false }, is_absorb_i1max { false }; bool is_absorb_i2min { false }, is_absorb_i2max { false }; @@ -181,10 +183,11 @@ namespace kernel::sr { const M& metric, const F& external_fields, const E& emission_policy, - const PBC& custom_prtl_bc) + const PUPD& custom_prtl_update) : pusher_flags { pusher_params.pusher_flags } , radiative_drag_flags { pusher_params.radiative_drag_flags } , EB { EB } + , sp { pusher_arrays.sp } , i1 { pusher_arrays.i1 } , i2 { pusher_arrays.i2 } , i3 { pusher_arrays.i3 } @@ -206,7 +209,7 @@ namespace kernel::sr { , metric { metric } , external_fields { external_fields } , emission_policy { emission_policy } - , custom_prtl_bc { custom_prtl_bc } + , custom_prtl_update { custom_prtl_update } , time { pusher_params.time } , coeff { HALF * (pusher_params.charge / pusher_params.mass) * pusher_params.omegaB0 * pusher_params.dt } @@ -653,6 +656,11 @@ namespace kernel::sr { from_Xi_to_i_di(xp[2], i3(p), dx3(p)); } } + + // call custom particle position update + custom_prtl_update(p, xp, *this); + + // apply boundary conditions boundaryConditions(p, xp); } @@ -1375,7 +1383,6 @@ namespace kernel::sr { dx1(p) = ONE - dx1(p); invert_vel = true; } else if (is_custom_i1min) { - custom_prtl_bc(p, 1, true, *this); } } else if (i1(p) >= ni1) { if (is_periodic_i1max) { @@ -1388,7 +1395,6 @@ namespace kernel::sr { dx1(p) = ONE - dx1(p); invert_vel = true; } else if (is_custom_i1max) { - custom_prtl_bc(p, 1, false, *this); } } if (invert_vel) { @@ -1421,7 +1427,6 @@ namespace kernel::sr { dx2(p) = ONE - dx2(p); invert_vel = true; } else if (is_custom_i2min) { - custom_prtl_bc(p, 2, true, *this); } else if (is_axis_i2min) { i2(p) = 0; dx2(p) = ONE - dx2(p); @@ -1437,7 +1442,6 @@ namespace kernel::sr { dx2(p) = ONE - dx2(p); invert_vel = true; } else if (is_custom_i2max) { - custom_prtl_bc(p, 2, false, *this); } else if (is_axis_i2max) { i2(p) = ni2 - 1; dx2(p) = ONE - dx2(p); @@ -1473,7 +1477,6 @@ namespace kernel::sr { dx3(p) = ONE - dx3(p); invert_vel = true; } else if (is_custom_i3min) { - custom_prtl_bc(p, 3, true, *this); } } else if (i3(p) >= ni3) { if (is_periodic_i3max) { @@ -1486,7 +1489,6 @@ namespace kernel::sr { dx3(p) = ONE - dx3(p); invert_vel = true; } else if (is_custom_i3max) { - custom_prtl_bc(p, 3, false, *this); } } if (invert_vel) { From 077f9088fffc4b8b5aa651ef2cd3afd0ccd0f3a7 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 31 Mar 2026 20:08:10 +0000 Subject: [PATCH 03/14] remove redundant class default value --- src/kernels/particle_pusher_sr.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kernels/particle_pusher_sr.hpp b/src/kernels/particle_pusher_sr.hpp index 18f7659ce..3e1c80166 100644 --- a/src/kernels/particle_pusher_sr.hpp +++ b/src/kernels/particle_pusher_sr.hpp @@ -111,7 +111,7 @@ namespace kernel::sr { * @tparam M Metric * @tparam F Additional force */ - template > + template requires metric::traits::HasD && metric::traits::HasTransformXYZ && metric::traits::HasConvertXYZ && metric::traits::HasTransform_i && metric::traits::HasConvert_i From 211fa865f79c61880dac15b5f0f0a8bdc455b183 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 31 Mar 2026 21:02:49 +0000 Subject: [PATCH 04/14] example for `CustomParticleUpdate` --- .../particle_update.toml | 82 ++++++ .../examples/custom_particle_update/pgen.hpp | 234 ++++++++++++++++++ 2 files changed, 316 insertions(+) create mode 100644 pgens/examples/custom_particle_update/particle_update.toml create mode 100644 pgens/examples/custom_particle_update/pgen.hpp diff --git a/pgens/examples/custom_particle_update/particle_update.toml b/pgens/examples/custom_particle_update/particle_update.toml new file mode 100644 index 000000000..017ba82bd --- /dev/null +++ b/pgens/examples/custom_particle_update/particle_update.toml @@ -0,0 +1,82 @@ +[simulation] + name = "_reflect" + engine = "srpic" + runtime = 100.0 + +[grid] + resolution = [256, 256] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["CONDUCTOR", "CONDUCTOR"], ["PERIODIC"]] + particles = [["PERIODIC"], ["PERIODIC"]] + +[scales] + larmor0 = 1.0 + skindepth0 = 1.0 + +[algorithms] + current_filters = 8 + + [algorithms.timestep] + CFL = 0.7071 + +[particles] + ppc0 = 10.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e2 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e2 + +[setup] + x1_e = [0.9] + x2_e = [0.8] + x3_e = [0.0] + ux1_e = [3.162] + ux2_e = [2.0] + ux3_e = [0.0] + + x1_i = [0.1] + x2_i = [0.2] + x3_i = [0.0] + ux1_i = [-3.162] + ux2_i = [-2.0] + ux3_i = [0.0] + + temperature = 0.04 + temperature_gradient = 10.0 + +[output] + format = "BPFile" + interval_time = 0.04 + + [output.fields] + quantities = ["E"] + + [output.particles] + enable = true + stride = 1 + + [output.spectra] + enable = false + + [output.debug] + ghosts = false + +[checkpoint] + keep = 1 + +[diagnostics] + log_level = "WARNING" + blocking_timers = true diff --git a/pgens/examples/custom_particle_update/pgen.hpp b/pgens/examples/custom_particle_update/pgen.hpp new file mode 100644 index 000000000..2162a909b --- /dev/null +++ b/pgens/examples/custom_particle_update/pgen.hpp @@ -0,0 +1,234 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "arch/traits.h" + +#include "archetypes/particle_injector.h" +#include "archetypes/problem_generator.h" +#include "archetypes/energy_dist.h" +#include "framework/domain/domain.h" +#include "framework/domain/metadomain.h" + +#include + +/* -------------------------------------------------------------------------- */ +/* Local macros (same as in particle_pusher_sr.hpp) */ +/* -------------------------------------------------------------------------- */ +#define from_Xi_to_i(XI, I) \ + { \ + I = static_cast((XI + 1)) - 1; \ + } + +#define from_Xi_to_i_di(XI, I, DI) \ + { \ + from_Xi_to_i((XI), (I)); \ + DI = static_cast((XI)) - static_cast(I); \ + } + +#define i_di_to_Xi(I, DI) static_cast((I)) + static_cast((DI)) + + +namespace user { + using namespace ntt; + + template + struct PGen : public arch::ProblemGenerator { + + // compatibility traits for the problem generator + static constexpr auto engines { + arch::traits::pgen::compatible_with::value + }; + static constexpr auto metrics { + arch::traits::pgen::compatible_with::value + }; + static constexpr auto dimensions { + arch::traits::pgen::compatible_with::value + }; + + // for easy access to variables in the child class + using arch::ProblemGenerator::D; + using arch::ProblemGenerator::C; + using arch::ProblemGenerator::params; + + const Metadomain& global_domain; + + const real_t temperature, temperature_gradient; + + inline PGen(const SimulationParams& p, const Metadomain& global_domain) + : arch::ProblemGenerator { p } + , global_domain { global_domain } + , temperature { params.template get("setup.temperature", 0.0) } + , temperature_gradient { params.template get("setup.temperature_gradient", 0.0) } + {} + + inline void InitPrtls(Domain& local_domain) { + const auto empty = std::vector {}; + const auto x1_e = params.template get>("setup.x1_e", + empty); + const auto x2_e = params.template get>("setup.x2_e", + empty); + const auto x3_e = params.template get>("setup.x3_e", + empty); + const auto phi_e = params.template get>("setup.phi_e", + empty); + const auto ux1_e = params.template get>("setup.ux1_e", + empty); + const auto ux2_e = params.template get>("setup.ux2_e", + empty); + const auto ux3_e = params.template get>("setup.ux3_e", + empty); + + const auto x1_i = params.template get>("setup.x1_i", + empty); + const auto x2_i = params.template get>("setup.x2_i", + empty); + const auto x3_i = params.template get>("setup.x3_i", + empty); + const auto phi_i = params.template get>("setup.phi_i", + empty); + const auto ux1_i = params.template get>("setup.ux1_i", + empty); + const auto ux2_i = params.template get>("setup.ux2_i", + empty); + const auto ux3_i = params.template get>("setup.ux3_i", + empty); + std::map> data_e { + { "x1", x1_e }, + { "x2", x2_e }, + { "ux1", ux1_e }, + { "ux2", ux2_e }, + { "ux3", ux3_e } + }; + std::map> data_i { + { "x1", x1_i }, + { "x2", x2_i }, + { "ux1", ux1_i }, + { "ux2", ux2_i }, + { "ux3", ux3_i } + }; + if constexpr (M::CoordType == Coord::Cart or D == Dim::_3D) { + data_e["x3"] = x3_e; + data_i["x3"] = x3_i; + } else if constexpr (D == Dim::_2D) { + data_e["phi"] = phi_e; + data_i["phi"] = phi_i; + } + + arch::InjectGlobally(global_domain, local_domain, (spidx_t)1, data_e); + arch::InjectGlobally(global_domain, local_domain, (spidx_t)2, data_i); + } + + auto FixFieldsConst(const bc_in&, const em&) const -> std::pair { + return { ZERO, false }; + } + + struct CustomPrtlUpdate { + random_number_pool_t pool; + real_t temp_cold, temp_hot; + real_t xmin, xmax; + + template + Inline void operator()(index_t p, Coord& xp, PusherKernel& pusher) const { + + const auto x_Cd = static_cast(pusher.i1(p)) + + static_cast(pusher.dx1(p)); + const auto x_Ph = pusher.metric.template convert<1, Crd::Cd, Crd::XYZ>(x_Cd); + vec_t v {ZERO}; + const coord_t x_dummy { ZERO }; + + // step 1: calculate the particle 3 velocity + const real_t gamma_p = math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); + + const real_t beta_x_p = pusher.ux1(p)/gamma_p; + const real_t xp_prev = i_di_to_Xi(pusher.i1_prev(p), pusher.dx1_prev(p)); + + // Reflecting boundary that resamples velocity + if (x_Ph < xmin) { + arch::JuttnerSinge(v, temp_cold, pool); + + // step 2: calculate the time for the particle to reach the wall + const int delta_i1_to_wall = pusher.i1_prev(p); + const prtldx_t delta_dx1_to_wall = pusher.dx1_prev(p); + const real_t dx_to_wall = i_di_to_Xi(delta_i1_to_wall, delta_dx1_to_wall); + const real_t dt_to_wall = dx_to_wall / pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, beta_x_p); + + // update the velocity to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) + pusher.ux1(p) = math::abs(v[0]); + pusher.ux2(p) = v[1]; + pusher.ux3(p) = v[2]; + + // step 3: calculate remaining time after the collision + const real_t remaining_dt = pusher.dt - dt_to_wall; + const real_t remaining_dt_inv_energy = remaining_dt / math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); + + pusher.i1(p) = 0; + pusher.dx1(p) = pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, pusher.ux1(p))*remaining_dt_inv_energy ; + + // // alternative: simply kill the particle at the boundary without reflection + // pusher.ux1(p) = ZERO; + // pusher.ux2(p) = ZERO; + // pusher.ux3(p) = ZERO; + + // pusher.i1(p) = pusher.i1_prev(p); + // pusher.dx1(p) = pusher.dx1_prev(p); + + // Re-sync coordinate variable xp for subsequent boundary procedures + xp[0] = static_cast(pusher.i1(p)) + + static_cast(pusher.dx1(p)); + + + } else if (x_Ph > xmax) { + arch::JuttnerSinge(v, temp_hot, pool); + + // step 2: calculate the time for the particle to reach the piston + const int delta_i1_to_wall = pusher.ni1 - pusher.i1_prev(p); + const prtldx_t delta_dx1_to_wall = ONE - pusher.dx1_prev(p); + const real_t dx_to_wall = i_di_to_Xi(delta_i1_to_wall, delta_dx1_to_wall); + + const real_t dt_to_wall = dx_to_wall / pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, beta_x_p); + + pusher.ux1(p) = -math::abs(v[0]); + pusher.ux2(p) = v[1]; + pusher.ux3(p) = v[2]; + + // step 3: calculate remaining time after the collision + const real_t remaining_dt = pusher.dt - dt_to_wall; + const real_t remaining_dt_inv_energy = remaining_dt / math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); + + + pusher.i1(p) = pusher.ni1 - 1; + pusher.dx1(p) = ONE - pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, math::abs(pusher.ux1(p)))*remaining_dt_inv_energy; + + // // alternative: simply kill the particle at the boundary without reflection + // pusher.ux1(p) = ZERO; + // pusher.ux2(p) = ZERO; + // pusher.ux3(p) = ZERO; + + // pusher.i1(p) = pusher.i1_prev(p); + // pusher.dx1(p) = pusher.dx1_prev(p); + + // Re-sync coordinate variable xp for subsequent boundary procedures + xp[0] = static_cast(pusher.i1(p)) + + static_cast(pusher.dx1(p)); + + } + } + }; + + template + auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const { + return CustomPrtlUpdate { domain.random_pool(), + temperature / domain.species[sp].mass(), + temperature_gradient * temperature / domain.species[sp].mass(), + global_domain.mesh().extent(in::x1).first, // xmin + global_domain.mesh().extent(in::x1).second // xmax + }; + } + }; + +} // namespace user + +#endif // PROBLEM_GENERATOR_H \ No newline at end of file From 27761e1c28d5d373368a6e59a374765068c87406 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 31 Mar 2026 21:31:04 +0000 Subject: [PATCH 05/14] removed remaining is_custom_in if clauses --- src/kernels/particle_pusher_sr.hpp | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/src/kernels/particle_pusher_sr.hpp b/src/kernels/particle_pusher_sr.hpp index 3e1c80166..c15a89e0d 100644 --- a/src/kernels/particle_pusher_sr.hpp +++ b/src/kernels/particle_pusher_sr.hpp @@ -163,9 +163,6 @@ namespace kernel::sr { bool is_reflect_i1min { false }, is_reflect_i1max { false }; bool is_reflect_i2min { false }, is_reflect_i2max { false }; bool is_reflect_i3min { false }, is_reflect_i3max { false }; - bool is_custom_i1min { false }, is_custom_i1max { false }; - bool is_custom_i2min { false }, is_custom_i2max { false }; - bool is_custom_i3min { false }, is_custom_i3max { false }; bool is_axis_i2min { false }, is_axis_i2max { false }; // gca parameters @@ -271,8 +268,6 @@ namespace kernel::sr { is_periodic_i1max = (pusher_params.boundaries[0].second == PrtlBC::PERIODIC); is_reflect_i1min = (pusher_params.boundaries[0].first == PrtlBC::REFLECT); is_reflect_i1max = (pusher_params.boundaries[0].second == PrtlBC::REFLECT); - is_custom_i1min = (pusher_params.boundaries[0].first == PrtlBC::CUSTOM); - is_custom_i1max = (pusher_params.boundaries[0].second == PrtlBC::CUSTOM); if constexpr ((D == Dim::_2D) || (D == Dim::_3D)) { raise::ErrorIf(pusher_params.boundaries.size() < 2, "pusher_params.boundaries defined incorrectly", @@ -286,8 +281,6 @@ namespace kernel::sr { is_periodic_i2max = (pusher_params.boundaries[1].second == PrtlBC::PERIODIC); is_reflect_i2min = (pusher_params.boundaries[1].first == PrtlBC::REFLECT); is_reflect_i2max = (pusher_params.boundaries[1].second == PrtlBC::REFLECT); - is_custom_i2min = (pusher_params.boundaries[1].first == PrtlBC::CUSTOM); - is_custom_i2max = (pusher_params.boundaries[1].second == PrtlBC::CUSTOM); is_axis_i2min = (pusher_params.boundaries[1].first == PrtlBC::AXIS); is_axis_i2max = (pusher_params.boundaries[1].second == PrtlBC::AXIS); } @@ -304,8 +297,6 @@ namespace kernel::sr { is_periodic_i3max = (pusher_params.boundaries[2].second == PrtlBC::PERIODIC); is_reflect_i3min = (pusher_params.boundaries[2].first == PrtlBC::REFLECT); is_reflect_i3max = (pusher_params.boundaries[2].second == PrtlBC::REFLECT); - is_custom_i3min = (pusher_params.boundaries[2].first == PrtlBC::CUSTOM); - is_custom_i3max = (pusher_params.boundaries[2].second == PrtlBC::CUSTOM); } } @@ -1382,7 +1373,6 @@ namespace kernel::sr { i1(p) = 0; dx1(p) = ONE - dx1(p); invert_vel = true; - } else if (is_custom_i1min) { } } else if (i1(p) >= ni1) { if (is_periodic_i1max) { @@ -1394,7 +1384,6 @@ namespace kernel::sr { i1(p) = ni1 - 1; dx1(p) = ONE - dx1(p); invert_vel = true; - } else if (is_custom_i1max) { } } if (invert_vel) { @@ -1426,7 +1415,6 @@ namespace kernel::sr { i2(p) = 0; dx2(p) = ONE - dx2(p); invert_vel = true; - } else if (is_custom_i2min) { } else if (is_axis_i2min) { i2(p) = 0; dx2(p) = ONE - dx2(p); @@ -1441,7 +1429,6 @@ namespace kernel::sr { i2(p) = ni2 - 1; dx2(p) = ONE - dx2(p); invert_vel = true; - } else if (is_custom_i2max) { } else if (is_axis_i2max) { i2(p) = ni2 - 1; dx2(p) = ONE - dx2(p); @@ -1476,7 +1463,6 @@ namespace kernel::sr { i3(p) = 0; dx3(p) = ONE - dx3(p); invert_vel = true; - } else if (is_custom_i3min) { } } else if (i3(p) >= ni3) { if (is_periodic_i3max) { @@ -1488,7 +1474,6 @@ namespace kernel::sr { i3(p) = ni3 - 1; dx3(p) = ONE - dx3(p); invert_vel = true; - } else if (is_custom_i3max) { } } if (invert_vel) { From dedae54dc8cd91e08252b7ea2b5a67d0a2378e41 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 31 Mar 2026 22:26:34 +0000 Subject: [PATCH 06/14] bugfixes in example --- .../examples/custom_particle_update/pgen.hpp | 41 ++++++------------- 1 file changed, 13 insertions(+), 28 deletions(-) diff --git a/pgens/examples/custom_particle_update/pgen.hpp b/pgens/examples/custom_particle_update/pgen.hpp index 2162a909b..34310e037 100644 --- a/pgens/examples/custom_particle_update/pgen.hpp +++ b/pgens/examples/custom_particle_update/pgen.hpp @@ -142,14 +142,14 @@ namespace user { // step 1: calculate the particle 3 velocity const real_t gamma_p = math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); - const real_t beta_x_p = pusher.ux1(p)/gamma_p; + const real_t beta_x_p = math::abs(pusher.ux1(p))/gamma_p; const real_t xp_prev = i_di_to_Xi(pusher.i1_prev(p), pusher.dx1_prev(p)); // Reflecting boundary that resamples velocity if (x_Ph < xmin) { arch::JuttnerSinge(v, temp_cold, pool); - // step 2: calculate the time for the particle to reach the wall + // calculate the time for the particle to reach the wall const int delta_i1_to_wall = pusher.i1_prev(p); const prtldx_t delta_dx1_to_wall = pusher.dx1_prev(p); const real_t dx_to_wall = i_di_to_Xi(delta_i1_to_wall, delta_dx1_to_wall); @@ -160,36 +160,29 @@ namespace user { pusher.ux2(p) = v[1]; pusher.ux3(p) = v[2]; - // step 3: calculate remaining time after the collision + // calculate remaining time after the collision const real_t remaining_dt = pusher.dt - dt_to_wall; const real_t remaining_dt_inv_energy = remaining_dt / math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); - + + // update the position to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) pusher.i1(p) = 0; - pusher.dx1(p) = pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, pusher.ux1(p))*remaining_dt_inv_energy ; - - // // alternative: simply kill the particle at the boundary without reflection - // pusher.ux1(p) = ZERO; - // pusher.ux2(p) = ZERO; - // pusher.ux3(p) = ZERO; - - // pusher.i1(p) = pusher.i1_prev(p); - // pusher.dx1(p) = pusher.dx1_prev(p); + pusher.dx1(p) = pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, pusher.ux1(p))*remaining_dt_inv_energy; // Re-sync coordinate variable xp for subsequent boundary procedures xp[0] = static_cast(pusher.i1(p)) + static_cast(pusher.dx1(p)); - } else if (x_Ph > xmax) { arch::JuttnerSinge(v, temp_hot, pool); // step 2: calculate the time for the particle to reach the piston - const int delta_i1_to_wall = pusher.ni1 - pusher.i1_prev(p); + const int delta_i1_to_wall = pusher.ni1 - 1 - pusher.i1_prev(p); const prtldx_t delta_dx1_to_wall = ONE - pusher.dx1_prev(p); const real_t dx_to_wall = i_di_to_Xi(delta_i1_to_wall, delta_dx1_to_wall); const real_t dt_to_wall = dx_to_wall / pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, beta_x_p); + // update the velocity to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) pusher.ux1(p) = -math::abs(v[0]); pusher.ux2(p) = v[1]; pusher.ux3(p) = v[2]; @@ -198,31 +191,23 @@ namespace user { const real_t remaining_dt = pusher.dt - dt_to_wall; const real_t remaining_dt_inv_energy = remaining_dt / math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); - - pusher.i1(p) = pusher.ni1 - 1; + // update the position to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) + pusher.i1(p) = pusher.ni1 - 2; pusher.dx1(p) = ONE - pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, math::abs(pusher.ux1(p)))*remaining_dt_inv_energy; - // // alternative: simply kill the particle at the boundary without reflection - // pusher.ux1(p) = ZERO; - // pusher.ux2(p) = ZERO; - // pusher.ux3(p) = ZERO; - - // pusher.i1(p) = pusher.i1_prev(p); - // pusher.dx1(p) = pusher.dx1_prev(p); - // Re-sync coordinate variable xp for subsequent boundary procedures xp[0] = static_cast(pusher.i1(p)) + static_cast(pusher.dx1(p)); - } } }; template auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const { + return CustomPrtlUpdate { domain.random_pool(), - temperature / domain.species[sp].mass(), - temperature_gradient * temperature / domain.species[sp].mass(), + temperature, // / domain.species[sp].mass(), + temperature_gradient * temperature, // / domain.species[sp].mass(), global_domain.mesh().extent(in::x1).first, // xmin global_domain.mesh().extent(in::x1).second // xmax }; From cba3e358896ccc43975c91241ba7d5352e78f513 Mon Sep 17 00:00:00 2001 From: A Sullivan Date: Wed, 1 Apr 2026 10:50:07 -0700 Subject: [PATCH 07/14] piston archetypes and pgen added --- pgens/piston_shock/pgen.hpp | 142 ++++++++++++++++++++++++++++++ src/archetypes/piston.h | 171 ++++++++++++++++++++++++++++++++++++ 2 files changed, 313 insertions(+) create mode 100644 pgens/piston_shock/pgen.hpp create mode 100644 src/archetypes/piston.h diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp new file mode 100644 index 000000000..86839f0c9 --- /dev/null +++ b/pgens/piston_shock/pgen.hpp @@ -0,0 +1,142 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "arch/traits.h" +#include "archetypes/piston.h" + +#include "archetypes/particle_injector.h" +#include "archetypes/problem_generator.h" +#include "archetypes/energy_dist.h" +#include "framework/domain/domain.h" +#include "archetypes/spatial_dist.h" +#include "framework/domain/metadomain.h" +#include "archetypes/utils.h" + +#include + + + +namespace user { + using namespace ntt; + + template + struct PGen : public arch::ProblemGenerator { + + // compatibility traits for the problem generator + static constexpr auto engines { + arch::traits::pgen::compatible_with::value + }; + static constexpr auto metrics { + arch::traits::pgen::compatible_with::value + }; + static constexpr auto dimensions { + arch::traits::pgen::compatible_with::value + }; + + // for easy access to variables in the child class + using arch::ProblemGenerator::D; + using arch::ProblemGenerator::C; + using arch::ProblemGenerator::params; + + const Metadomain& global_domain; + + // domain properties + const real_t global_xmin, global_xmax; + + // plasma + const real_t temperature, temperature_ratio; + // piston properties + const real_t piston_velocity, piston_initial_position; + + inline PGen(const SimulationParams& p, const Metadomain& global_domain) + : arch::ProblemGenerator { p } + , global_domain { global_domain } + , global_xmin { global_domain.mesh().extent(in::x1).first } + , global_xmax { global_domain.mesh().extent(in::x1).second } + , piston_velocity {p.template get("setup.piston_velocity", 0.0)} + , piston_initial_position {p.template get("setup.piston_initial_position", 0.0)} + , temperature { p.template get("setup.temperature", 0.0) } + , temperature_ratio { p.template get("setup.temperature_ratio") }{} + + inline void InitPrtls(Domain& local_domain) { + real_t xg_min = global_xmin+piston_initial_position; + real_t xg_max = global_xmax ; + + // define box to inject into + boundaries_t box; + // loop over all dimensions + for (auto d { 0u }; d < (unsigned int)M::Dim; ++d) { + // compute the range for the x-direction + if (d == static_cast(in::x1)) { + box.push_back({ xg_min, xg_max }); + } else { + // inject into full range in other directions + box.push_back(Range::All); + } + } + + // define temperatures of species + const auto temperatures = std::make_pair(temperature, + temperature_ratio * temperature); + // define drift speed of species + const auto drifts = std::make_pair( + std::vector { ZERO, ZERO, ZERO }, + std::vector { ZERO, ZERO, ZERO }); + + // inject particles + arch::InjectUniformMaxwellians(params, + local_domain, + ONE, + temperatures, + { 1, 2 }, + drifts, + false, + box); + } + + + + struct CustomPrtlUpdate { + real_t piston_position_current; // current position of piston + real_t piston_velocity_current; + real_t xg_max; + + bool is_left; + + bool massive; + + + + template + Inline void operator()(index_t p, Coord& xp, PusherKernel& pusher) const { + + real_t piston_position_use; + + if (piston_position_current(p, pusher, piston_position_use, piston_velocity_current, is_left)){ + arch::PistonUpdate(p, pusher, piston_position_use, piston_velocity_current, massive); + } + + + } + }; + + template + auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const { + return CustomPrtlUpdate { piston_initial_position+static_cast(time)*piston_velocity, piston_velocity, global_domain.mesh().extent(in::x1).second, true, true}; + + }; + + }; + +} // namespace user + +#endif // PROBLEM_GENERATOR_H \ No newline at end of file diff --git a/src/archetypes/piston.h b/src/archetypes/piston.h new file mode 100644 index 000000000..a33aa49d8 --- /dev/null +++ b/src/archetypes/piston.h @@ -0,0 +1,171 @@ +/** + * @file archetypes/piston.h + * @brief Piston functions for implementing the piston in the CustomPrtlUpdate in pgran + * @implements + * - arch::PistonUpdate<> -> void + * - arch::CrossesPiston<> -> Bool + * @namespaces: + * - arch:: + */ + +#ifndef ARCHETYPES_PISTON_H +#define ARCHETYPES_PISTON_H + +#include "enums.h" +#include "global.h" + +#include "archetypes/energy_dist.h" +#include "framework/domain/domain.h" +#include "framework/parameters/parameters.h" +#include "framework/domain/metadomain.h" + +#include + +/* -------------------------------------------------------------------------- */ +/* Local macros (same as in particle_pusher_sr.hpp) */ +/* -------------------------------------------------------------------------- */ +#define from_Xi_to_i(XI, I) \ + { \ + I = static_cast((XI + 1)) - 1; \ + } + +#define from_Xi_to_i_di(XI, I, DI) \ + { \ + from_Xi_to_i((XI), (I)); \ + DI = static_cast((XI)) - static_cast(I); \ + } + +#define i_di_to_Xi(I, DI) static_cast((I)) + static_cast((DI)) + +/* -------------------------------------------------------------------------- */ + +namespace arch { + + /** + * @brief Updates particle position and velocity if it reflects off a moiving piston, called to correct patticle position + * + * @param p Index of particle + * @param pusher Particle pusher engine for particle update + * @param piston_position Position of the piston at the start of timestep in global coordinates + * @param piston_v Velocity of piston at current timestep + */ + template + Inline void PistonUpdate(const index_t p, + const PusherKernel& pusher, + const real_t piston_position, + const real_t piston_v, + const bool massive){ + + coord_t xp_Cd { ZERO }; + pusher.getParticlePosition(p, xp_Cd); + + // step 1: calculate the particle 3 velocity + const real_t gamma_p { + massive + ? (math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) + : (math::sqrt(SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) + }; + + const real_t beta_x_p = pusher.ux1(p)/gamma_p; + const real_t xp_prev = i_di_to_Xi(pusher.i1_prev(p), pusher.dx1_prev(p)); + + const real_t piston_position_local = pusher.metric.template convert<1, Crd::Ph, Crd::Cd>(piston_position); + int i_w ; + prtldx_t dx_w; + + from_Xi_to_i_di(piston_position_local, i_w, dx_w); + + const real_t piston_gamma = ONE/math::sqrt(ONE-SQR(piston_v)); + + // step 2: calculate the time for the particle to reach the piston + const int delta_i1_to_piston = (i_w - pusher.i1_prev(p)); + const prtldx_t delta_dx1_to_piston = (dx_w-pusher.dx1_prev(p)); + const real_t dx_to_piston = i_di_to_Xi(delta_i1_to_piston, delta_dx1_to_piston); + + const real_t dt_to_piston = dx_to_piston / pusher.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, beta_x_p-piston_v); + // step 3: calculate remaining time after the collision + const real_t remaining_dt = pusher.dt - dt_to_piston; + // step 4: update the particle's velocity and position after the collision + + pusher.ux1(p) = -SQR(piston_gamma)*( (ONE+SQR(piston_v))*pusher.ux1(p)-TWO * piston_v*gamma_p ); + + + const real_t remaining_dt_inv_energy { + massive + ? (remaining_dt / math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) + : (remaining_dt / math::sqrt(SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) + }; + // define piston integer and fractional coordinate + + int i_w_coll = i_w; + prtldx_t dx_w_coll = dx_w + pusher.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, piston_v)*dt_to_piston; + + i_w_coll += static_cast( dx_w_coll >= ONE) - + static_cast( dx_w_coll< ZERO); + dx_w_coll -= (dx_w_coll >= ONE); + dx_w_coll+= (dx_w_coll< ZERO); + + pusher.i1(p) = i_w_coll; + pusher.dx1(p) = pusher.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, pusher.ux1(p))*remaining_dt_inv_energy+dx_w_coll ; + + pusher.i1(p) += static_cast(pusher.dx1(p) >= ONE) - + static_cast(pusher.dx1(p) < ZERO); + pusher.dx1(p) -= (pusher.dx1(p) >= ONE); + pusher.dx1(p) += (pusher.dx1(p) < ZERO); + + + } + /** + * @brief Checks whether a particle reflects off a moving piston, called after particle has been moved by regular pusher + * + * @param p Index of particle + * @param pusher Particle pusher engine for particle update + * @param piston_position Position of the piston at the start of timestep in global coordinates + * @param piston_v Velocity of piston at current timestep + * @param is_left Is piston on the left side of the box or right side of the box + */ + template + Inline bool CrossesPiston(const index_t p, + const PusherKernel& pusher, + const real_t piston_position, + const real_t piston_v, + const bool is_left ){ + const real_t x1_Cd = i_di_to_Xi(pusher.i1(p), pusher.dx1(p)); + coord_t xp_Cd { ZERO }; + pusher.getParticlePosition(p, xp_Cd); + // x1_Cd_wallmove is not the actual particle coordinate + // it is particle position minus how much the wall has moved in this timestep + // This is a computational trick + const real_t x1_Cd_wallmove = x1_Cd-pusher.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, piston_v)*pusher.dt; + const real_t x1_Ph_wallmove = pusher.metric.template convert<1, Crd::Cd, Crd::Ph>(x1_Cd_wallmove); + + + if (is_left){ // if piston is moving from left, ask if particle is to the left of piston + if(piston_position>x1_Ph_wallmove){ + return true; + } else { + return false; + } + } else { // if piston is moving from the right, so ask is particle to right of piston + if(piston_position Date: Wed, 1 Apr 2026 15:55:33 -0700 Subject: [PATCH 08/14] piston_shock.toml added --- pgens/piston_shock/piston_shock.toml | 75 ++++++++++++++++++++++++++++ pgens/piston_shock/shock.toml | 74 +++++++++++++++++++++++++++ 2 files changed, 149 insertions(+) create mode 100644 pgens/piston_shock/piston_shock.toml create mode 100644 pgens/piston_shock/shock.toml diff --git a/pgens/piston_shock/piston_shock.toml b/pgens/piston_shock/piston_shock.toml new file mode 100644 index 000000000..6a81d1ab5 --- /dev/null +++ b/pgens/piston_shock/piston_shock.toml @@ -0,0 +1,75 @@ +[simulation] + name = "piston_shock" + engine = "srpic" + runtime = 100.0 + +[grid] + resolution = [300, 24] + extent = [[0.0, 100.0], [0.0, 8.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["CONDUCTOR", "CONDUCTOR"], ["PERIODIC"]] + particles = [["REFLECT", "REFLECT"], ["PERIODIC"]] + +[scales] + larmor0 = 5.7735 + skindepth0 = 1.0 + +[algorithms] + current_filters = 8 + + [algorithms.timestep] + CFL = 0.7 + +[particles] + ppc0 = 8.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 2e7 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = 1.0 + maxnpart = 2e7 + +[setup] + piston_velocity = 0.99 # speed of leftmoving piston + piston_initial_position = 0.0 # speed towards the wall [c] + temperature = 0.168 # temperature of maxwell distribution [kB T / (m_i c^2)] + temperature_ratio = 1.0 # temperature ratio of electrons to protons + #Bmag = 1.0 # magnetic field strength as fraction of magnetisation + #Btheta = 63.0 # magnetic field angle in the plane + #Bphi = 0.0 # magnetic field angle out of plane + #filling_fraction = 0.99 # fraction of the shock piston filled with plasma + #injector_velocity = 0.0 # speed of injector [c] + #injection_start = 0.0 # start time of moving injector + #injection_frequency = 100 + +[output] + interval_time = 1000.0 + format = "BPFile" + + [output.fields] + quantities = ["N", "B", "E"] + + [output.particles] + enable = true + stride = 10 + + [output.spectra] + enable = false + +[diagnostics] + log_level = "WARNING" + blocking_timers = true + +[checkpoint] + interval = 10000 + keep = 1 \ No newline at end of file diff --git a/pgens/piston_shock/shock.toml b/pgens/piston_shock/shock.toml new file mode 100644 index 000000000..8cde86990 --- /dev/null +++ b/pgens/piston_shock/shock.toml @@ -0,0 +1,74 @@ +[simulation] + name = "shock" + engine = "srpic" + runtime = 15200.0 + +[grid] + resolution = [32768, 256] + extent = [[0.0, 8192.0], [-32.0, 32.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["CONDUCTOR", "MATCH"], ["PERIODIC"]] + particles = [["REFLECT", "ABSORB"], ["PERIODIC"]] + +[scales] + larmor0 = 5.7735 + skindepth0 = 1.0 + +[algorithms] + current_filters = 8 + + [algorithms.timestep] + CFL = 0.7 + +[particles] + ppc0 = 8.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 2e7 + + [[particles.species]] + label = "p+" + mass = 100.0 + charge = 1.0 + maxnpart = 2e7 + +[setup] + drift_ux = 0.15 # speed towards the wall [c] + temperature = 0.168 # temperature of maxwell distribution [kB T / (m_i c^2)] + temperature_ratio = 1.0 # temperature ratio of electrons to protons + Bmag = 1.0 # magnetic field strength as fraction of magnetisation + Btheta = 63.0 # magnetic field angle in the plane + Bphi = 0.0 # magnetic field angle out of plane + filling_fraction = 0.99 # fraction of the shock piston filled with plasma + injector_velocity = 0.0 # speed of injector [c] + injection_start = 0.0 # start time of moving injector + injection_frequency = 100 + +[output] + interval_time = 1000.0 + format = "BPFile" + + [output.fields] + quantities = ["N", "B", "E"] + + [output.particles] + enable = true + stride = 10 + + [output.spectra] + enable = false + +[diagnostics] + log_level = "WARNING" + blocking_timers = true + +[checkpoint] + interval = 10000 + keep = 1 \ No newline at end of file From 35d261f8a77021169669f881521e38f66656e0cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Mon, 6 Apr 2026 11:47:25 -0500 Subject: [PATCH 09/14] revert to previous version --- src/kernels/digital_filter.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/kernels/digital_filter.hpp b/src/kernels/digital_filter.hpp index 7f3e52f1f..5ac60327d 100644 --- a/src/kernels/digital_filter.hpp +++ b/src/kernels/digital_filter.hpp @@ -93,8 +93,8 @@ namespace kernel { ? (boundaries[2].second == FldsBC::CONDUCTOR) : false } , i1_max { size_[0] + N_GHOSTS } - , i2_max { (short)D > 1 ? (size_[(short)D > 1 ? 1 : 0] + N_GHOSTS) : 0 } - , i3_max { (short)D > 2 ? (size_[(short)D > 2 ? 2 : 0] + N_GHOSTS) : 0 } {} + , i2_max { (short)D > 1 ? (size_[1] + N_GHOSTS) : 0 } + , i3_max { (short)D > 2 ? (size_[2] + N_GHOSTS) : 0 } {} Inline void operator()(index_t i1) const { if constexpr ((D == Dim::_1D) && (C == Coord::Cart)) { From a0dcff3a9a4dd8fc26e997a1bdab5f89bdcbcad7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Mon, 6 Apr 2026 11:54:49 -0500 Subject: [PATCH 10/14] applied clang-format --- .../examples/custom_particle_update/pgen.hpp | 91 ++++++---- pgens/piston_shock/pgen.hpp | 70 ++++---- src/archetypes/piston.h | 160 ++++++++++-------- src/archetypes/traits.h | 6 +- src/engines/srpic/particle_pusher.h | 29 ++-- src/kernels/particle_pusher_sr.hpp | 75 ++++---- 6 files changed, 238 insertions(+), 193 deletions(-) diff --git a/pgens/examples/custom_particle_update/pgen.hpp b/pgens/examples/custom_particle_update/pgen.hpp index 34310e037..76a63b1e1 100644 --- a/pgens/examples/custom_particle_update/pgen.hpp +++ b/pgens/examples/custom_particle_update/pgen.hpp @@ -6,9 +6,9 @@ #include "arch/traits.h" +#include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" #include "archetypes/problem_generator.h" -#include "archetypes/energy_dist.h" #include "framework/domain/domain.h" #include "framework/domain/metadomain.h" @@ -30,7 +30,6 @@ #define i_di_to_Xi(I, DI) static_cast((I)) + static_cast((DI)) - namespace user { using namespace ntt; @@ -59,19 +58,20 @@ namespace user { inline PGen(const SimulationParams& p, const Metadomain& global_domain) : arch::ProblemGenerator { p } - , global_domain { global_domain } + , global_domain { global_domain } , temperature { params.template get("setup.temperature", 0.0) } - , temperature_gradient { params.template get("setup.temperature_gradient", 0.0) } - {} + , temperature_gradient { + params.template get("setup.temperature_gradient", 0.0) + } {} inline void InitPrtls(Domain& local_domain) { const auto empty = std::vector {}; const auto x1_e = params.template get>("setup.x1_e", - empty); + empty); const auto x2_e = params.template get>("setup.x2_e", - empty); + empty); const auto x3_e = params.template get>("setup.x3_e", - empty); + empty); const auto phi_e = params.template get>("setup.phi_e", empty); const auto ux1_e = params.template get>("setup.ux1_e", @@ -82,11 +82,11 @@ namespace user { empty); const auto x1_i = params.template get>("setup.x1_i", - empty); + empty); const auto x2_i = params.template get>("setup.x2_i", - empty); + empty); const auto x3_i = params.template get>("setup.x3_i", - empty); + empty); const auto phi_i = params.template get>("setup.phi_i", empty); const auto ux1_i = params.template get>("setup.ux1_i", @@ -127,22 +127,23 @@ namespace user { struct CustomPrtlUpdate { random_number_pool_t pool; - real_t temp_cold, temp_hot; - real_t xmin, xmax; + real_t temp_cold, temp_hot; + real_t xmin, xmax; template Inline void operator()(index_t p, Coord& xp, PusherKernel& pusher) const { const auto x_Cd = static_cast(pusher.i1(p)) + - static_cast(pusher.dx1(p)); + static_cast(pusher.dx1(p)); const auto x_Ph = pusher.metric.template convert<1, Crd::Cd, Crd::XYZ>(x_Cd); - vec_t v {ZERO}; + vec_t v { ZERO }; const coord_t x_dummy { ZERO }; // step 1: calculate the particle 3 velocity - const real_t gamma_p = math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); + const real_t gamma_p = math::sqrt( + ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); - const real_t beta_x_p = math::abs(pusher.ux1(p))/gamma_p; + const real_t beta_x_p = math::abs(pusher.ux1(p)) / gamma_p; const real_t xp_prev = i_di_to_Xi(pusher.i1_prev(p), pusher.dx1_prev(p)); // Reflecting boundary that resamples velocity @@ -150,10 +151,13 @@ namespace user { arch::JuttnerSinge(v, temp_cold, pool); // calculate the time for the particle to reach the wall - const int delta_i1_to_wall = pusher.i1_prev(p); + const int delta_i1_to_wall = pusher.i1_prev(p); const prtldx_t delta_dx1_to_wall = pusher.dx1_prev(p); const real_t dx_to_wall = i_di_to_Xi(delta_i1_to_wall, delta_dx1_to_wall); - const real_t dt_to_wall = dx_to_wall / pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, beta_x_p); + const real_t dt_to_wall = + dx_to_wall / + pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, + beta_x_p); // update the velocity to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) pusher.ux1(p) = math::abs(v[0]); @@ -161,12 +165,19 @@ namespace user { pusher.ux3(p) = v[2]; // calculate remaining time after the collision - const real_t remaining_dt = pusher.dt - dt_to_wall; - const real_t remaining_dt_inv_energy = remaining_dt / math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); - + const real_t remaining_dt = pusher.dt - dt_to_wall; + const real_t remaining_dt_inv_energy = remaining_dt / + math::sqrt(ONE + + SQR(pusher.ux1(p)) + + SQR(pusher.ux2(p)) + + SQR(pusher.ux3(p))); + // update the position to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) pusher.i1(p) = 0; - pusher.dx1(p) = pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, pusher.ux1(p))*remaining_dt_inv_energy; + pusher.dx1(p) = pusher.metric.template transform<1, Idx::XYZ, Idx::U>( + x_dummy, + pusher.ux1(p)) * + remaining_dt_inv_energy; // Re-sync coordinate variable xp for subsequent boundary procedures xp[0] = static_cast(pusher.i1(p)) + @@ -176,11 +187,14 @@ namespace user { arch::JuttnerSinge(v, temp_hot, pool); // step 2: calculate the time for the particle to reach the piston - const int delta_i1_to_wall = pusher.ni1 - 1 - pusher.i1_prev(p); + const int delta_i1_to_wall = pusher.ni1 - 1 - pusher.i1_prev(p); const prtldx_t delta_dx1_to_wall = ONE - pusher.dx1_prev(p); const real_t dx_to_wall = i_di_to_Xi(delta_i1_to_wall, delta_dx1_to_wall); - - const real_t dt_to_wall = dx_to_wall / pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, beta_x_p); + + const real_t dt_to_wall = + dx_to_wall / + pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, + beta_x_p); // update the velocity to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) pusher.ux1(p) = -math::abs(v[0]); @@ -188,12 +202,20 @@ namespace user { pusher.ux3(p) = v[2]; // step 3: calculate remaining time after the collision - const real_t remaining_dt = pusher.dt - dt_to_wall; - const real_t remaining_dt_inv_energy = remaining_dt / math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); + const real_t remaining_dt = pusher.dt - dt_to_wall; + const real_t remaining_dt_inv_energy = remaining_dt / + math::sqrt(ONE + + SQR(pusher.ux1(p)) + + SQR(pusher.ux2(p)) + + SQR(pusher.ux3(p))); // update the position to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) pusher.i1(p) = pusher.ni1 - 2; - pusher.dx1(p) = ONE - pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, math::abs(pusher.ux1(p)))*remaining_dt_inv_energy; + pusher.dx1(p) = ONE - + pusher.metric.template transform<1, Idx::XYZ, Idx::U>( + x_dummy, + math::abs(pusher.ux1(p))) * + remaining_dt_inv_energy; // Re-sync coordinate variable xp for subsequent boundary procedures xp[0] = static_cast(pusher.i1(p)) + @@ -205,11 +227,12 @@ namespace user { template auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const { - return CustomPrtlUpdate { domain.random_pool(), - temperature, // / domain.species[sp].mass(), - temperature_gradient * temperature, // / domain.species[sp].mass(), - global_domain.mesh().extent(in::x1).first, // xmin - global_domain.mesh().extent(in::x1).second // xmax + return CustomPrtlUpdate { + domain.random_pool(), + temperature, // / domain.species[sp].mass(), + temperature_gradient * temperature, // / domain.species[sp].mass(), + global_domain.mesh().extent(in::x1).first, // xmin + global_domain.mesh().extent(in::x1).second // xmax }; } }; diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index 86839f0c9..745f751d2 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -5,20 +5,18 @@ #include "global.h" #include "arch/traits.h" -#include "archetypes/piston.h" +#include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" +#include "archetypes/piston.h" #include "archetypes/problem_generator.h" -#include "archetypes/energy_dist.h" -#include "framework/domain/domain.h" #include "archetypes/spatial_dist.h" -#include "framework/domain/metadomain.h" #include "archetypes/utils.h" +#include "framework/domain/domain.h" +#include "framework/domain/metadomain.h" #include - - namespace user { using namespace ntt; @@ -44,7 +42,7 @@ namespace user { const Metadomain& global_domain; // domain properties - const real_t global_xmin, global_xmax; + const real_t global_xmin, global_xmax; // plasma const real_t temperature, temperature_ratio; @@ -53,17 +51,20 @@ namespace user { inline PGen(const SimulationParams& p, const Metadomain& global_domain) : arch::ProblemGenerator { p } - , global_domain { global_domain } + , global_domain { global_domain } , global_xmin { global_domain.mesh().extent(in::x1).first } , global_xmax { global_domain.mesh().extent(in::x1).second } - , piston_velocity {p.template get("setup.piston_velocity", 0.0)} - , piston_initial_position {p.template get("setup.piston_initial_position", 0.0)} + , piston_velocity { p.template get("setup.piston_velocity", 0.0) } + , piston_initial_position { p.template get( + "setup.piston_initial_position", + 0.0) } , temperature { p.template get("setup.temperature", 0.0) } - , temperature_ratio { p.template get("setup.temperature_ratio") }{} + , temperature_ratio { p.template get( + "setup.temperature_ratio") } {} inline void InitPrtls(Domain& local_domain) { - real_t xg_min = global_xmin+piston_initial_position; - real_t xg_max = global_xmax ; + real_t xg_min = global_xmin + piston_initial_position; + real_t xg_max = global_xmax; // define box to inject into boundaries_t box; @@ -82,9 +83,8 @@ namespace user { const auto temperatures = std::make_pair(temperature, temperature_ratio * temperature); // define drift speed of species - const auto drifts = std::make_pair( - std::vector { ZERO, ZERO, ZERO }, - std::vector { ZERO, ZERO, ZERO }); + const auto drifts = std::make_pair(std::vector { ZERO, ZERO, ZERO }, + std::vector { ZERO, ZERO, ZERO }); // inject particles arch::InjectUniformMaxwellians(params, @@ -97,8 +97,6 @@ namespace user { box); } - - struct CustomPrtlUpdate { real_t piston_position_current; // current position of piston real_t piston_velocity_current; @@ -107,34 +105,42 @@ namespace user { bool is_left; bool massive; - - template Inline void operator()(index_t p, Coord& xp, PusherKernel& pusher) const { - + real_t piston_position_use; - if (piston_position_current(p, pusher, piston_position_use, piston_velocity_current, is_left)){ - arch::PistonUpdate(p, pusher, piston_position_use, piston_velocity_current, massive); + if (arch::CrossesPiston(p, + pusher, + piston_position_use, + piston_velocity_current, + is_left)) { + arch::PistonUpdate(p, + pusher, + piston_position_use, + piston_velocity_current, + massive); } - - } }; template auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const { - return CustomPrtlUpdate { piston_initial_position+static_cast(time)*piston_velocity, piston_velocity, global_domain.mesh().extent(in::x1).second, true, true}; - - }; - + return CustomPrtlUpdate { piston_initial_position + + static_cast(time) * piston_velocity, + piston_velocity, + global_domain.mesh().extent(in::x1).second, + true, + true }; + }; }; } // namespace user diff --git a/src/archetypes/piston.h b/src/archetypes/piston.h index a33aa49d8..171d1728d 100644 --- a/src/archetypes/piston.h +++ b/src/archetypes/piston.h @@ -16,8 +16,8 @@ #include "archetypes/energy_dist.h" #include "framework/domain/domain.h" -#include "framework/parameters/parameters.h" #include "framework/domain/metadomain.h" +#include "framework/parameters/parameters.h" #include @@ -42,82 +42,98 @@ namespace arch { /** - * @brief Updates particle position and velocity if it reflects off a moiving piston, called to correct patticle position - * + * @brief Updates particle position and velocity if it reflects off a moiving + * piston, called to correct patticle position + * * @param p Index of particle * @param pusher Particle pusher engine for particle update - * @param piston_position Position of the piston at the start of timestep in global coordinates + * @param piston_position Position of the piston at the start of timestep in + * global coordinates * @param piston_v Velocity of piston at current timestep */ template - Inline void PistonUpdate(const index_t p, - const PusherKernel& pusher, - const real_t piston_position, - const real_t piston_v, - const bool massive){ - + Inline void PistonUpdate(const index_t p, + const PusherKernel& pusher, + const real_t piston_position, + const real_t piston_v, + const bool massive) { + coord_t xp_Cd { ZERO }; pusher.getParticlePosition(p, xp_Cd); - + // step 1: calculate the particle 3 velocity const real_t gamma_p { - massive - ? (math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) + massive + ? (math::sqrt( + ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) : (math::sqrt(SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) - }; - - const real_t beta_x_p = pusher.ux1(p)/gamma_p; - const real_t xp_prev = i_di_to_Xi(pusher.i1_prev(p), pusher.dx1_prev(p)); - - const real_t piston_position_local = pusher.metric.template convert<1, Crd::Ph, Crd::Cd>(piston_position); - int i_w ; - prtldx_t dx_w; + }; + + const real_t beta_x_p = pusher.ux1(p) / gamma_p; + const real_t xp_prev = i_di_to_Xi(pusher.i1_prev(p), pusher.dx1_prev(p)); + + const real_t piston_position_local = + pusher.metric.template convert<1, Crd::Ph, Crd::Cd>(piston_position); + int i_w; + prtldx_t dx_w; from_Xi_to_i_di(piston_position_local, i_w, dx_w); - const real_t piston_gamma = ONE/math::sqrt(ONE-SQR(piston_v)); + const real_t piston_gamma = ONE / math::sqrt(ONE - SQR(piston_v)); // step 2: calculate the time for the particle to reach the piston - const int delta_i1_to_piston = (i_w - pusher.i1_prev(p)); - const prtldx_t delta_dx1_to_piston = (dx_w-pusher.dx1_prev(p)); + const int delta_i1_to_piston = (i_w - pusher.i1_prev(p)); + const prtldx_t delta_dx1_to_piston = (dx_w - pusher.dx1_prev(p)); const real_t dx_to_piston = i_di_to_Xi(delta_i1_to_piston, delta_dx1_to_piston); - - const real_t dt_to_piston = dx_to_piston / pusher.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, beta_x_p-piston_v); + + const real_t dt_to_piston = dx_to_piston / + pusher.metric.template transform<1, Idx::XYZ, Idx::U>( + xp_Cd, + beta_x_p - piston_v); // step 3: calculate remaining time after the collision const real_t remaining_dt = pusher.dt - dt_to_piston; // step 4: update the particle's velocity and position after the collision - pusher.ux1(p) = -SQR(piston_gamma)*( (ONE+SQR(piston_v))*pusher.ux1(p)-TWO * piston_v*gamma_p ); + pusher.ux1(p) = -SQR(piston_gamma) * ((ONE + SQR(piston_v)) * pusher.ux1(p) - + TWO * piston_v * gamma_p); - const real_t remaining_dt_inv_energy { - massive - ? (remaining_dt / math::sqrt(ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) - : (remaining_dt / math::sqrt(SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) - }; + massive + ? (remaining_dt / math::sqrt(ONE + SQR(pusher.ux1(p)) + + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) + : (remaining_dt / math::sqrt(SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + + SQR(pusher.ux3(p)))) + }; // define piston integer and fractional coordinate - int i_w_coll = i_w; - prtldx_t dx_w_coll = dx_w + pusher.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, piston_v)*dt_to_piston; + int i_w_coll = i_w; + prtldx_t dx_w_coll = dx_w + + pusher.metric.template transform<1, Idx::XYZ, Idx::U>( + xp_Cd, + piston_v) * + dt_to_piston; - i_w_coll += static_cast( dx_w_coll >= ONE) - - static_cast( dx_w_coll< ZERO); + i_w_coll += static_cast(dx_w_coll >= ONE) - + static_cast(dx_w_coll < ZERO); dx_w_coll -= (dx_w_coll >= ONE); - dx_w_coll+= (dx_w_coll< ZERO); + dx_w_coll += (dx_w_coll < ZERO); - pusher.i1(p) = i_w_coll; - pusher.dx1(p) = pusher.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, pusher.ux1(p))*remaining_dt_inv_energy+dx_w_coll ; + pusher.i1(p) = i_w_coll; + pusher.dx1(p) = pusher.metric.template transform<1, Idx::XYZ, Idx::U>( + xp_Cd, + pusher.ux1(p)) * + remaining_dt_inv_energy + + dx_w_coll; - pusher.i1(p) += static_cast(pusher.dx1(p) >= ONE) - - static_cast(pusher.dx1(p) < ZERO); + pusher.i1(p) += static_cast(pusher.dx1(p) >= ONE) - + static_cast(pusher.dx1(p) < ZERO); pusher.dx1(p) -= (pusher.dx1(p) >= ONE); pusher.dx1(p) += (pusher.dx1(p) < ZERO); + } - - } /** * @brief Checks whether a particle reflects off a moving piston, called after particle has been moved by regular pusher - * + * * @param p Index of particle * @param pusher Particle pusher engine for particle update * @param piston_position Position of the piston at the start of timestep in global coordinates @@ -125,45 +141,41 @@ namespace arch { * @param is_left Is piston on the left side of the box or right side of the box */ template - Inline bool CrossesPiston(const index_t p, - const PusherKernel& pusher, - const real_t piston_position, - const real_t piston_v, - const bool is_left ){ - const real_t x1_Cd = i_di_to_Xi(pusher.i1(p), pusher.dx1(p)); + Inline bool CrossesPiston(const index_t p, + const PusherKernel& pusher, + const real_t piston_position, + const real_t piston_v, + const bool is_left) { + const real_t x1_Cd = i_di_to_Xi(pusher.i1(p), pusher.dx1(p)); coord_t xp_Cd { ZERO }; pusher.getParticlePosition(p, xp_Cd); // x1_Cd_wallmove is not the actual particle coordinate - // it is particle position minus how much the wall has moved in this timestep - // This is a computational trick - const real_t x1_Cd_wallmove = x1_Cd-pusher.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, piston_v)*pusher.dt; - const real_t x1_Ph_wallmove = pusher.metric.template convert<1, Crd::Cd, Crd::Ph>(x1_Cd_wallmove); - - - if (is_left){ // if piston is moving from left, ask if particle is to the left of piston - if(piston_position>x1_Ph_wallmove){ - return true; - } else { - return false; - } + // it is particle position minus how much the wall has moved in this + // timestep This is a computational trick + const real_t x1_Cd_wallmove = + x1_Cd - + pusher.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, piston_v) * + pusher.dt; + const real_t x1_Ph_wallmove = pusher.metric.template convert<1, Crd::Cd, Crd::Ph>( + x1_Cd_wallmove); + + if (is_left) { // if piston is moving from left, ask if particle is to the left of piston + if (piston_position > x1_Ph_wallmove) { + return true; + } else { + return false; + } } else { // if piston is moving from the right, so ask is particle to right of piston - if(piston_position concept HasCustomPrtlUpdate = requires(const PG& pgen, - simtime_t time, - spidx_t sp, - D& domain) { + simtime_t time, + spidx_t sp, + D& domain) { pgen.CustomParticleUpdate(time, sp, domain); }; diff --git a/src/engines/srpic/particle_pusher.h b/src/engines/srpic/particle_pusher.h index abd5f3a59..a63cd4134 100644 --- a/src/engines/srpic/particle_pusher.h +++ b/src/engines/srpic/particle_pusher.h @@ -26,20 +26,23 @@ namespace ntt { namespace srpic { template - void CallPusher_WithExternalFieldFlag(Domain& domain, - const SimulationParams& params, - const kernel::sr::PusherParams& pusher_params, - kernel::sr::PusherArrays& pusher_arrays, - EmissionTypeFlag emission_policy_flag, - const range_t& range, - const ndfield_t& EB, - const M& metric, - const PG& pgen, - const F& external_fields) { + void CallPusher_WithExternalFieldFlag( + Domain& domain, + const SimulationParams& params, + const kernel::sr::PusherParams& pusher_params, + kernel::sr::PusherArrays& pusher_arrays, + EmissionTypeFlag emission_policy_flag, + const range_t& range, + const ndfield_t& EB, + const M& metric, + const PG& pgen, + const F& external_fields) { auto get_custom_prtl_update = [&]() { - if constexpr (arch::traits::pgen::HasCustomPrtlUpdate) { - return pgen.CustomParticleUpdate( - pusher_params.time, pusher_params.species_index, domain); + if constexpr ( + arch::traits::pgen::HasCustomPrtlUpdate) { + return pgen.CustomParticleUpdate(pusher_params.time, + pusher_params.species_index, + domain); } else { return kernel::sr::NoCustomPrtlUpdate_t {}; } diff --git a/src/kernels/particle_pusher_sr.hpp b/src/kernels/particle_pusher_sr.hpp index c15a89e0d..31055eff0 100644 --- a/src/kernels/particle_pusher_sr.hpp +++ b/src/kernels/particle_pusher_sr.hpp @@ -104,7 +104,8 @@ namespace kernel::sr { template struct NoCustomPrtlUpdate_t { template - Inline void operator()(index_t, coord_t&, const PusherKernel&) const {} + Inline void operator()(index_t, coord_t&, const PusherKernel&) const { + } }; /** @@ -133,8 +134,8 @@ namespace kernel::sr { public: const spidx_t sp; - simtime_t time; - const real_t dt; + simtime_t time; + const real_t dt; array_t i1, i2, i3; array_t i1_prev, i2_prev, i3_prev; array_t dx1, dx2, dx3; @@ -144,26 +145,26 @@ namespace kernel::sr { array_t weight; array_t tag; const M metric; - const int ni1, ni2, ni3; + const int ni1, ni2, ni3; private: - const F external_fields; + const F external_fields; - const E emission_policy; + const E emission_policy; const PUPD custom_prtl_update; const real_t coeff; - bool is_absorb_i1min { false }, is_absorb_i1max { false }; - bool is_absorb_i2min { false }, is_absorb_i2max { false }; - bool is_absorb_i3min { false }, is_absorb_i3max { false }; - bool is_periodic_i1min { false }, is_periodic_i1max { false }; - bool is_periodic_i2min { false }, is_periodic_i2max { false }; - bool is_periodic_i3min { false }, is_periodic_i3max { false }; - bool is_reflect_i1min { false }, is_reflect_i1max { false }; - bool is_reflect_i2min { false }, is_reflect_i2max { false }; - bool is_reflect_i3min { false }, is_reflect_i3max { false }; - bool is_axis_i2min { false }, is_axis_i2max { false }; + bool is_absorb_i1min { false }, is_absorb_i1max { false }; + bool is_absorb_i2min { false }, is_absorb_i2max { false }; + bool is_absorb_i3min { false }, is_absorb_i3max { false }; + bool is_periodic_i1min { false }, is_periodic_i1max { false }; + bool is_periodic_i2min { false }, is_periodic_i2max { false }; + bool is_periodic_i3min { false }, is_periodic_i3max { false }; + bool is_reflect_i1min { false }, is_reflect_i1max { false }; + bool is_reflect_i2min { false }, is_reflect_i2max { false }; + bool is_reflect_i3min { false }, is_reflect_i3max { false }; + bool is_axis_i2min { false }, is_axis_i2max { false }; // gca parameters const real_t gca_larmor, gca_EovrB_sqr; @@ -579,8 +580,8 @@ namespace kernel::sr { dx1_prev(p) = dx1(p); dx1(p) += metric.template transform<1, Idx::XYZ, Idx::U>(xp, ux1(p)) * dt_inv_energy; - i1(p) += static_cast(dx1(p) >= ONE) - - static_cast(dx1(p) < ZERO); + i1(p) += static_cast(dx1(p) >= ONE) - + static_cast(dx1(p) < ZERO); dx1(p) -= (dx1(p) >= ONE); dx1(p) += (dx1(p) < ZERO); } @@ -589,8 +590,8 @@ namespace kernel::sr { dx2_prev(p) = dx2(p); dx2(p) += metric.template transform<2, Idx::XYZ, Idx::U>(xp, ux2(p)) * dt_inv_energy; - i2(p) += static_cast(dx2(p) >= ONE) - - static_cast(dx2(p) < ZERO); + i2(p) += static_cast(dx2(p) >= ONE) - + static_cast(dx2(p) < ZERO); dx2(p) -= (dx2(p) >= ONE); dx2(p) += (dx2(p) < ZERO); } @@ -599,8 +600,8 @@ namespace kernel::sr { dx3_prev(p) = dx3(p); dx3(p) += metric.template transform<3, Idx::XYZ, Idx::U>(xp, ux3(p)) * dt_inv_energy; - i3(p) += static_cast(dx3(p) >= ONE) - - static_cast(dx3(p) < ZERO); + i3(p) += static_cast(dx3(p) >= ONE) - + static_cast(dx3(p) < ZERO); dx3(p) -= (dx3(p) >= ONE); dx3(p) += (dx3(p) < ZERO); } @@ -611,8 +612,8 @@ namespace kernel::sr { : ONE / math::sqrt(SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p))) }; vec_t vp_Cart { ux1(p) * inv_energy, - ux2(p) * inv_energy, - ux3(p) * inv_energy }; + ux2(p) * inv_energy, + ux3(p) * inv_energy }; // get cartesian position coord_t xp_Cart { ZERO }; metric.template convert_xyz(xp, xp_Cart); @@ -650,7 +651,7 @@ namespace kernel::sr { // call custom particle position update custom_prtl_update(p, xp, *this); - + // apply boundary conditions boundaryConditions(p, xp); } @@ -716,11 +717,11 @@ namespace kernel::sr { auto COEFF2 { ONE + NORM_SQR(u1[0], u1[1], u1[2]) - NORM_SQR(b0[0], b0[1], b0[2]) }; - COEFF = ONE / - math::sqrt( - INV_2 * (COEFF2 + math::sqrt(SQR(COEFF2) + - FOUR * (SQR(b0[0]) + SQR(b0[1]) + - SQR(b0[2]) + SQR(COEFF))))); + COEFF = ONE / + math::sqrt( + INV_2 * (COEFF2 + math::sqrt(SQR(COEFF2) + + FOUR * (SQR(b0[0]) + SQR(b0[1]) + + SQR(b0[2]) + SQR(COEFF))))); COEFF2 = ONE / (ONE + SQR(b0[0] * COEFF) + SQR(b0[1] * COEFF) + SQR(b0[2] * COEFF)); @@ -1515,8 +1516,8 @@ namespace kernel::sr { const vec_t& e0, const vec_t& b0) const { real_t gamma_prime_sqr = ONE / math::sqrt(ONE + NORM_SQR(u_prime[0], - u_prime[1], - u_prime[2])); + u_prime[1], + u_prime[2])); u_prime[0] *= gamma_prime_sqr; u_prime[1] *= gamma_prime_sqr; u_prime[2] *= gamma_prime_sqr; @@ -1566,8 +1567,8 @@ namespace kernel::sr { Inline void inverseComptonDrag(index_t p, vec_t& u_prime) const { real_t gamma_prime_sqr = ONE / math::sqrt(ONE + NORM_SQR(u_prime[0], - u_prime[1], - u_prime[2])); + u_prime[1], + u_prime[2])); u_prime[0] *= gamma_prime_sqr; u_prime[1] *= gamma_prime_sqr; u_prime[2] *= gamma_prime_sqr; @@ -1588,7 +1589,7 @@ namespace kernel::sr { { typename E::Payload payload; vec_t delta_u_Ph { ZERO }; - const auto emission_response = emission_policy.shouldEmit(xp_Cd, + const auto emission_response = emission_policy.shouldEmit(xp_Cd, xp_Ph, u_prime, ep_Cart, @@ -1605,8 +1606,8 @@ namespace kernel::sr { if (emission_response.first) { vec_t direction { ZERO }; const auto delta_u_Ph_mag = NORM(delta_u_Ph[0], - delta_u_Ph[1], - delta_u_Ph[2]); + delta_u_Ph[1], + delta_u_Ph[2]); direction[0] = -delta_u_Ph[0] / delta_u_Ph_mag; direction[1] = -delta_u_Ph[1] / delta_u_Ph_mag; direction[2] = -delta_u_Ph[2] / delta_u_Ph_mag; From 1013069e173f8cb1e8b68e9cc23f9a043e103d3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Mon, 6 Apr 2026 16:37:37 -0500 Subject: [PATCH 11/14] correct mass assignment --- .../examples/custom_particle_update/pgen.hpp | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/pgens/examples/custom_particle_update/pgen.hpp b/pgens/examples/custom_particle_update/pgen.hpp index 76a63b1e1..9e5a8d756 100644 --- a/pgens/examples/custom_particle_update/pgen.hpp +++ b/pgens/examples/custom_particle_update/pgen.hpp @@ -154,12 +154,11 @@ namespace user { const int delta_i1_to_wall = pusher.i1_prev(p); const prtldx_t delta_dx1_to_wall = pusher.dx1_prev(p); const real_t dx_to_wall = i_di_to_Xi(delta_i1_to_wall, delta_dx1_to_wall); - const real_t dt_to_wall = - dx_to_wall / + const real_t dt_to_wall = dx_to_wall / pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, beta_x_p); - // update the velocity to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) + // update the velocity to the post-collision value (reflection with new speed sampled from Juttner distribution) pusher.ux1(p) = math::abs(v[0]); pusher.ux2(p) = v[1]; pusher.ux3(p) = v[2]; @@ -172,7 +171,7 @@ namespace user { SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); - // update the position to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) + // update the position to the post-collision value (reflection with new speed sampled from Juttner distribution) pusher.i1(p) = 0; pusher.dx1(p) = pusher.metric.template transform<1, Idx::XYZ, Idx::U>( x_dummy, @@ -190,13 +189,11 @@ namespace user { const int delta_i1_to_wall = pusher.ni1 - 1 - pusher.i1_prev(p); const prtldx_t delta_dx1_to_wall = ONE - pusher.dx1_prev(p); const real_t dx_to_wall = i_di_to_Xi(delta_i1_to_wall, delta_dx1_to_wall); - - const real_t dt_to_wall = - dx_to_wall / + const real_t dt_to_wall = dx_to_wall / pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, beta_x_p); - // update the velocity to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) + // update the velocity to the post-collision value (reflection with new speed sampled from Juttner distribution) pusher.ux1(p) = -math::abs(v[0]); pusher.ux2(p) = v[1]; pusher.ux3(p) = v[2]; @@ -209,7 +206,7 @@ namespace user { SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); - // update the position to the post-collision value (perfect reflection with new speed sampled from Juttner distribution) + // update the position to the post-collision value (reflection with new speed sampled from Juttner distribution) pusher.i1(p) = pusher.ni1 - 2; pusher.dx1(p) = ONE - pusher.metric.template transform<1, Idx::XYZ, Idx::U>( @@ -229,8 +226,8 @@ namespace user { return CustomPrtlUpdate { domain.random_pool(), - temperature, // / domain.species[sp].mass(), - temperature_gradient * temperature, // / domain.species[sp].mass(), + temperature / domain.species[sp - 1].mass(), // sp is 1-indexed + temperature_gradient * temperature / domain.species[sp - 1].mass(), global_domain.mesh().extent(in::x1).first, // xmin global_domain.mesh().extent(in::x1).second // xmax }; From 0078b502af3d4ebc44314180d3cd473485e2631c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Mon, 6 Apr 2026 16:38:53 -0500 Subject: [PATCH 12/14] added missing undefs of macros --- pgens/examples/custom_particle_update/pgen.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pgens/examples/custom_particle_update/pgen.hpp b/pgens/examples/custom_particle_update/pgen.hpp index 9e5a8d756..9d7113d42 100644 --- a/pgens/examples/custom_particle_update/pgen.hpp +++ b/pgens/examples/custom_particle_update/pgen.hpp @@ -236,4 +236,8 @@ namespace user { } // namespace user +#undef from_Xi_to_i_di +#undef from_Xi_to_i +#undef i_di_to_Xi + #endif // PROBLEM_GENERATOR_H \ No newline at end of file From 1abec0658a4228ef04672d25c40975b708812a3e Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 7 Apr 2026 19:29:56 +0000 Subject: [PATCH 13/14] update tests to match new pusher definition --- src/kernels/tests/custom_emission.cpp | 6 ++++-- src/kernels/tests/ext_force.cpp | 6 ++++-- src/kernels/tests/gca_pusher.cpp | 6 ++++-- src/kernels/tests/prtl_bc.cpp | 6 ++++-- src/kernels/tests/pusher.cpp | 11 +++++++---- 5 files changed, 23 insertions(+), 12 deletions(-) diff --git a/src/kernels/tests/custom_emission.cpp b/src/kernels/tests/custom_emission.cpp index 94d631122..dfc88deda 100644 --- a/src/kernels/tests/custom_emission.cpp +++ b/src/kernels/tests/custom_emission.cpp @@ -248,6 +248,7 @@ auto main(int argc, char* argv[]) -> int { }; ndfield_t EB { "EB", 128u + 2u * N_GHOSTS }; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t> {}; for (auto step = 0u; step < 7u; ++step) { pusher_params.time = static_cast(step) * delta_t; @@ -286,13 +287,14 @@ auto main(int argc, char* argv[]) -> int { Kokkos::parallel_for( "ParticlePusher", 2u, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, kernel::sr::NoField_t {}, - emission_policy)); + emission_policy, + no_custom_update)); const auto n_injected = emission_policy.numbers_injected(); emitted_species_1.set_counter(emitted_species_1.counter() + n_injected[0]); emitted_species_1.set_npart(emitted_species_1.npart() + n_injected[0]); diff --git a/src/kernels/tests/ext_force.cpp b/src/kernels/tests/ext_force.cpp index 4e1b4ac32..29b22d528 100644 --- a/src/kernels/tests/ext_force.cpp +++ b/src/kernels/tests/ext_force.cpp @@ -198,6 +198,7 @@ void testPusher(const std::vector& res) { const auto no_emission = kernel::NoEmissionPolicy_t> {}; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t> {}; for (auto t { 0u }; t < 100; ++t) { const real_t time = t * dt; @@ -206,13 +207,14 @@ void testPusher(const std::vector& res) { Kokkos::parallel_for( "pusher", CreateRangePolicy({ 0 }, { 2 }), - kernel::sr::Pusher_kernel, decltype(ext_force), false, decltype(no_emission)>( + kernel::sr::Pusher_kernel, decltype(ext_force), false, decltype(no_emission), decltype(no_custom_update)>( pusher_params, pusher_arrays, emfield, metric, ext_force, - no_emission)); + no_emission, + no_custom_update)); auto i1_prev_ = Kokkos::create_mirror_view(i1_prev); auto i2_prev_ = Kokkos::create_mirror_view(i2_prev); diff --git a/src/kernels/tests/gca_pusher.cpp b/src/kernels/tests/gca_pusher.cpp index a5bd34799..1c198e028 100644 --- a/src/kernels/tests/gca_pusher.cpp +++ b/src/kernels/tests/gca_pusher.cpp @@ -181,6 +181,7 @@ void testPusher(const std::vector& res) { pusher_arrays.tag = tag; const auto no_emission = kernel::NoEmissionPolicy_t> {}; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t> {}; for (auto t { 0u }; t < 2000; ++t) { pusher_params.time = t * dt; @@ -188,13 +189,14 @@ void testPusher(const std::vector& res) { Kokkos::parallel_for( "pusher", CreateRangePolicy({ 0 }, { 2 }), - kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission)>( + kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission), decltype(no_custom_update)>( pusher_params, pusher_arrays, emfield, metric, kernel::sr::NoField_t {}, - no_emission)); + no_emission, + no_custom_update)); auto ux1_ = Kokkos::create_mirror_view(ux1); auto ux2_ = Kokkos::create_mirror_view(ux2); diff --git a/src/kernels/tests/prtl_bc.cpp b/src/kernels/tests/prtl_bc.cpp index bbdfa57b8..c459baa51 100644 --- a/src/kernels/tests/prtl_bc.cpp +++ b/src/kernels/tests/prtl_bc.cpp @@ -266,18 +266,20 @@ void testPeriodicBC(const std::vector& res, pusher_arrays.phi = phi; pusher_arrays.tag = tag; const auto no_emission = kernel::NoEmissionPolicy_t {}; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t {}; for (auto n { 0 }; n < n_iter; ++n) { Kokkos::parallel_for( "pusher", CreateRangePolicy({ 0 }, { 2 }), - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, emfield, metric, kernel::sr::NoField_t {}, - no_emission)); + no_emission, + no_custom_update)); auto i1_ = Kokkos::create_mirror_view(i1); auto i2_ = Kokkos::create_mirror_view(i2); auto i3_ = Kokkos::create_mirror_view(i3); diff --git a/src/kernels/tests/pusher.cpp b/src/kernels/tests/pusher.cpp index 3bd951de2..ef5f69b8f 100644 --- a/src/kernels/tests/pusher.cpp +++ b/src/kernels/tests/pusher.cpp @@ -179,6 +179,7 @@ void testPusher(const std::vector& res) { pusher_arrays.tag = tag; const auto no_emission = kernel::NoEmissionPolicy_t> {}; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t> {}; for (auto t { 0u }; t < 2000; ++t) { const real_t time = t * dt; @@ -188,26 +189,28 @@ void testPusher(const std::vector& res) { Kokkos::parallel_for( "pusher", CreateRangePolicy({ 0 }, { 1 }), - kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission)>( + kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission), decltype(no_custom_update)>( pusher_params, pusher_arrays, emfield, metric, kernel::sr::NoField_t {}, - no_emission)); + no_emission, + no_custom_update)); pusher_params.pusher_flags = ParticlePusher::VAY; Kokkos::parallel_for( "pusher", CreateRangePolicy({ 1 }, { 2 }), - kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission)>( + kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission), decltype(no_custom_update)>( pusher_params, pusher_arrays, emfield, metric, kernel::sr::NoField_t {}, - no_emission)); + no_emission, + no_custom_update)); auto i1_prev_ = Kokkos::create_mirror_view(i1_prev); auto i2_prev_ = Kokkos::create_mirror_view(i2_prev); From 394b7083c165a0cbe2330bfa44a116e2dd411f76 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 7 Apr 2026 19:30:25 +0000 Subject: [PATCH 14/14] added test for CustomParticleUpdate --- src/kernels/tests/CMakeLists.txt | 1 + src/kernels/tests/custom_prtl_update.cpp | 291 +++++++++++++++++++++++ 2 files changed, 292 insertions(+) create mode 100644 src/kernels/tests/custom_prtl_update.cpp diff --git a/src/kernels/tests/CMakeLists.txt b/src/kernels/tests/CMakeLists.txt index 7b73642ac..f1438e108 100644 --- a/src/kernels/tests/CMakeLists.txt +++ b/src/kernels/tests/CMakeLists.txt @@ -37,3 +37,4 @@ gen_test(pusher) gen_test(ext_force) gen_test(reduced_stats) gen_test(custom_emission) +gen_test(custom_prtl_update) diff --git a/src/kernels/tests/custom_prtl_update.cpp b/src/kernels/tests/custom_prtl_update.cpp new file mode 100644 index 000000000..c26ad76bf --- /dev/null +++ b/src/kernels/tests/custom_prtl_update.cpp @@ -0,0 +1,291 @@ +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/formatting.h" +#include "utils/numeric.h" + +#include "metrics/minkowski.h" + +#include "kernels/emission/emission.hpp" +#include "kernels/particle_pusher_sr.hpp" + +#include +#include + +#include +#include +#include +#include +#include + +using namespace ntt; +using namespace metric; + +void errorIf(bool condition, const std::string& message = "") { + if (condition) { + throw std::runtime_error(message); + } +} + +Inline auto equal(real_t a, real_t b, const std::string& msg) -> bool { + if (not(math::abs(a - b) < 1e-4)) { + Kokkos::printf("%.12e != %.12e %s\n", a, b, msg.c_str()); + return false; + } + return true; +} + +template +void put_value(array_t& arr, T v, index_t p) { + auto h = Kokkos::create_mirror_view(arr); + Kokkos::deep_copy(h, arr); + h(p) = v; + Kokkos::deep_copy(arr, h); +} + +struct TestCustomPrtlUpdate { + template + Inline void operator()(index_t p, Coord& xp, const PusherKernel& pusher) const { + if constexpr (PusherKernel::D == Dim::_1D || PusherKernel::D == Dim::_2D || PusherKernel::D == Dim::_3D) { + auto invert_vel = false; + if (pusher.i1(p) < 0) { + pusher.i1(p) = 0; + pusher.dx1(p) = ONE - pusher.dx1(p); + invert_vel = true; + } else if (pusher.i1(p) >= pusher.ni1) { + pusher.i1(p) = pusher.ni1 - 1; + pusher.dx1(p) = ONE - pusher.dx1(p); + invert_vel = true; + } + if (invert_vel) { + pusher.ux1(p) = -pusher.ux1(p); + } + } + if constexpr (PusherKernel::D == Dim::_2D || PusherKernel::D == Dim::_3D) { + auto invert_vel = false; + if (pusher.i2(p) < 0) { + pusher.i2(p) = 0; + pusher.dx2(p) = ONE - pusher.dx2(p); + invert_vel = true; + } else if (pusher.i2(p) >= pusher.ni2) { + pusher.i2(p) = pusher.ni2 - 1; + pusher.dx2(p) = ONE - pusher.dx2(p); + invert_vel = true; + } + if (invert_vel) { + pusher.ux2(p) = -pusher.ux2(p); + } + } + if constexpr (PusherKernel::D == Dim::_3D) { + auto invert_vel = false; + if (pusher.i3(p) < 0) { + pusher.i3(p) = 0; + pusher.dx3(p) = ONE - pusher.dx3(p); + invert_vel = true; + } else if (pusher.i3(p) >= pusher.ni3) { + pusher.i3(p) = pusher.ni3 - 1; + pusher.dx3(p) = ONE - pusher.dx3(p); + invert_vel = true; + } + if (invert_vel) { + pusher.ux3(p) = -pusher.ux3(p); + } + } + } +}; + +template +void testCustomPrtlUpdate(const std::vector& res, + const boundaries_t& ext, + const std::map& params = {}) { + + errorIf(res.size() != M::Dim, "res.size() != M::Dim"); + errorIf(M::CoordType != Coord::Cart, "M::CoordType != Coord::Cart"); + + M metric { res, ext, params }; + + int nx1 = 0, nx2 = 0, nx3 = 0; + if (res.size() > 0) nx1 = static_cast(res.at(0)); + if (res.size() > 1) nx2 = static_cast(res.at(1)); + if (res.size() > 2) nx3 = static_cast(res.at(2)); + + const real_t dt = 0.1 * (ext.at(0).second - ext.at(0).first) / static_cast(nx1); + + ndfield_t emfield; + if constexpr (M::Dim == Dim::_1D) { + emfield = ndfield_t { "emfield", res.at(0) + 2 * N_GHOSTS }; + } else if constexpr (M::Dim == Dim::_2D) { + emfield = ndfield_t { "emfield", res.at(0) + 2 * N_GHOSTS, res.at(1) + 2 * N_GHOSTS }; + } else { + emfield = ndfield_t { "emfield", res.at(0) + 2 * N_GHOSTS, res.at(1) + 2 * N_GHOSTS, res.at(2) + 2 * N_GHOSTS }; + } + + // Arrays for REFLECT boundaries + array_t r_i1 { "r_i1", 1 }, r_i2 { "r_i2", 1 }, r_i3 { "r_i3", 1 }; + array_t r_i1_prev { "r_i1_prev", 1 }, r_i2_prev { "r_i2_prev", 1 }, r_i3_prev { "r_i3_prev", 1 }; + array_t r_dx1 { "r_dx1", 1 }, r_dx2 { "r_dx2", 1 }, r_dx3 { "r_dx3", 1 }; + array_t r_dx1_prev { "r_dx1_prev", 1 }, r_dx2_prev { "r_dx2_prev", 1 }, r_dx3_prev { "r_dx3_prev", 1 }; + array_t r_ux1 { "r_ux1", 1 }, r_ux2 { "r_ux2", 1 }, r_ux3 { "r_ux3", 1 }; + array_t r_tag { "r_tag", 1 }; + + // Arrays for CUSTOM boundaries + array_t c_i1 { "c_i1", 1 }, c_i2 { "c_i2", 1 }, c_i3 { "c_i3", 1 }; + array_t c_i1_prev { "c_i1_prev", 1 }, c_i2_prev { "c_i2_prev", 1 }, c_i3_prev { "c_i3_prev", 1 }; + array_t c_dx1 { "c_dx1", 1 }, c_dx2 { "c_dx2", 1 }, c_dx3 { "c_dx3", 1 }; + array_t c_dx1_prev { "c_dx1_prev", 1 }, c_dx2_prev { "c_dx2_prev", 1 }, c_dx3_prev { "c_dx3_prev", 1 }; + array_t c_ux1 { "c_ux1", 1 }, c_ux2 { "c_ux2", 1 }, c_ux3 { "c_ux3", 1 }; + array_t c_tag { "c_tag", 1 }; + + array_t phi; + + real_t ux_1 = 0.569197; + real_t uy_1 = 0.716085; + real_t uz_1 = -0.760101; + put_value(r_ux1, ux_1, 0); put_value(c_ux1, ux_1, 0); + put_value(r_ux2, uy_1, 0); put_value(c_ux2, uy_1, 0); + put_value(r_ux3, uz_1, 0); put_value(c_ux3, uz_1, 0); + + put_value(r_tag, ParticleTag::alive, 0); + put_value(c_tag, ParticleTag::alive, 0); + + if constexpr (M::PrtlDim >= Dim::_1D) { + put_value(r_i1, nx1 - 1, 0); put_value(c_i1, nx1 - 1, 0); + put_value(r_dx1, 0.9, 0); put_value(c_dx1, 0.9, 0); + } + if constexpr (M::PrtlDim >= Dim::_2D) { + put_value(r_i2, nx2 - 1, 0); put_value(c_i2, nx2 - 1, 0); + put_value(r_dx2, 0.9, 0); put_value(c_dx2, 0.9, 0); + } + if constexpr (M::PrtlDim == Dim::_3D) { + put_value(r_i3, 0, 0); put_value(c_i3, 0, 0); + put_value(r_dx3, 0.1, 0); put_value(c_dx3, 0.1, 0); + } + + kernel::sr::PusherParams r_params {}; + r_params.pusher_flags = ParticlePusher::BORIS; + r_params.mass = ONE; + r_params.charge = ONE; + r_params.dt = dt; + r_params.omegaB0 = ONE; + r_params.ni1 = nx1; + r_params.ni2 = nx2; + r_params.ni3 = nx3; + r_params.boundaries = { + { PrtlBC::REFLECT, PrtlBC::REFLECT }, + { PrtlBC::REFLECT, PrtlBC::REFLECT }, + { PrtlBC::REFLECT, PrtlBC::REFLECT } + }; + + kernel::sr::PusherParams c_params = r_params; + // initialize with periodic boundaries so that reflection is only handled by the custom update + c_params.boundaries = { + { PrtlBC::PERIODIC, PrtlBC::PERIODIC }, + { PrtlBC::PERIODIC, PrtlBC::PERIODIC }, + { PrtlBC::PERIODIC, PrtlBC::PERIODIC } + }; + + kernel::sr::PusherArrays r_arrays {}; + r_arrays.sp = 1u; + r_arrays.i1 = r_i1; r_arrays.i2 = r_i2; r_arrays.i3 = r_i3; + r_arrays.i1_prev = r_i1_prev; r_arrays.i2_prev = r_i2_prev; r_arrays.i3_prev = r_i3_prev; + r_arrays.dx1 = r_dx1; r_arrays.dx2 = r_dx2; r_arrays.dx3 = r_dx3; + r_arrays.dx1_prev = r_dx1_prev; r_arrays.dx2_prev = r_dx2_prev; r_arrays.dx3_prev = r_dx3_prev; + r_arrays.ux1 = r_ux1; r_arrays.ux2 = r_ux2; r_arrays.ux3 = r_ux3; + r_arrays.phi = phi; + r_arrays.tag = r_tag; + + kernel::sr::PusherArrays c_arrays {}; + c_arrays.sp = 1u; + c_arrays.i1 = c_i1; c_arrays.i2 = c_i2; c_arrays.i3 = c_i3; + c_arrays.i1_prev = c_i1_prev; c_arrays.i2_prev = c_i2_prev; c_arrays.i3_prev = c_i3_prev; + c_arrays.dx1 = c_dx1; c_arrays.dx2 = c_dx2; c_arrays.dx3 = c_dx3; + c_arrays.dx1_prev = c_dx1_prev; c_arrays.dx2_prev = c_dx2_prev; c_arrays.dx3_prev = c_dx3_prev; + c_arrays.ux1 = c_ux1; c_arrays.ux2 = c_ux2; c_arrays.ux3 = c_ux3; + c_arrays.phi = phi; + c_arrays.tag = c_tag; + + const auto no_emission = kernel::NoEmissionPolicy_t {}; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t {}; + const auto custom_update = TestCustomPrtlUpdate {}; + + const auto n_iter = 100; + + for (auto n { 0 }; n < n_iter; ++n) { + Kokkos::parallel_for( + "pusher_reflect", + CreateRangePolicy({ 0 }, { 1 }), + kernel::sr::Pusher_kernel( + r_params, r_arrays, emfield, metric, kernel::sr::NoField_t {}, no_emission, no_custom_update) + ); + Kokkos::parallel_for( + "pusher_custom", + CreateRangePolicy({ 0 }, { 1 }), + kernel::sr::Pusher_kernel( + c_params, c_arrays, emfield, metric, kernel::sr::NoField_t {}, no_emission, custom_update) + ); + + auto hr_i1 = Kokkos::create_mirror_view(r_i1); Kokkos::deep_copy(hr_i1, r_i1); + auto hc_i1 = Kokkos::create_mirror_view(c_i1); Kokkos::deep_copy(hc_i1, c_i1); + auto hr_dx1 = Kokkos::create_mirror_view(r_dx1); Kokkos::deep_copy(hr_dx1, r_dx1); + auto hc_dx1 = Kokkos::create_mirror_view(c_dx1); Kokkos::deep_copy(hc_dx1, c_dx1); + auto hr_ux1 = Kokkos::create_mirror_view(r_ux1); Kokkos::deep_copy(hr_ux1, r_ux1); + auto hc_ux1 = Kokkos::create_mirror_view(c_ux1); Kokkos::deep_copy(hc_ux1, c_ux1); + + if constexpr (M::PrtlDim >= Dim::_1D) { + errorIf(hr_i1(0) != hc_i1(0), "i1 mismatch"); + errorIf(not equal(hr_dx1(0), hc_dx1(0), "dx1 mismatch")); + errorIf(not equal(hr_ux1(0), hc_ux1(0), "ux1 mismatch")); + } + + if constexpr (M::PrtlDim >= Dim::_2D) { + auto hr_i2 = Kokkos::create_mirror_view(r_i2); Kokkos::deep_copy(hr_i2, r_i2); + auto hc_i2 = Kokkos::create_mirror_view(c_i2); Kokkos::deep_copy(hc_i2, c_i2); + auto hr_dx2 = Kokkos::create_mirror_view(r_dx2); Kokkos::deep_copy(hr_dx2, r_dx2); + auto hc_dx2 = Kokkos::create_mirror_view(c_dx2); Kokkos::deep_copy(hc_dx2, c_dx2); + auto hr_ux2 = Kokkos::create_mirror_view(r_ux2); Kokkos::deep_copy(hr_ux2, r_ux2); + auto hc_ux2 = Kokkos::create_mirror_view(c_ux2); Kokkos::deep_copy(hc_ux2, c_ux2); + errorIf(hr_i2(0) != hc_i2(0), "i2 mismatch"); + errorIf(not equal(hr_dx2(0), hc_dx2(0), "dx2 mismatch")); + errorIf(not equal(hr_ux2(0), hc_ux2(0), "ux2 mismatch")); + } + + if constexpr (M::PrtlDim == Dim::_3D) { + auto hr_i3 = Kokkos::create_mirror_view(r_i3); Kokkos::deep_copy(hr_i3, r_i3); + auto hc_i3 = Kokkos::create_mirror_view(c_i3); Kokkos::deep_copy(hc_i3, c_i3); + auto hr_dx3 = Kokkos::create_mirror_view(r_dx3); Kokkos::deep_copy(hr_dx3, r_dx3); + auto hc_dx3 = Kokkos::create_mirror_view(c_dx3); Kokkos::deep_copy(hc_dx3, c_dx3); + auto hr_ux3 = Kokkos::create_mirror_view(r_ux3); Kokkos::deep_copy(hr_ux3, r_ux3); + auto hc_ux3 = Kokkos::create_mirror_view(c_ux3); Kokkos::deep_copy(hc_ux3, c_ux3); + errorIf(hr_i3(0) != hc_i3(0), "i3 mismatch"); + errorIf(not equal(hr_dx3(0), hc_dx3(0), "dx3 mismatch")); + errorIf(not equal(hr_ux3(0), hc_ux3(0), "ux3 mismatch")); + } + } +} + +auto main(int argc, char* argv[]) -> int { + Kokkos::initialize(argc, argv); + + try { + using namespace ntt; + + const std::vector res1d { 50 }; + const boundaries_t ext1d { { 0.0, 1000.0 } }; + const std::vector res2d { 30, 20 }; + const boundaries_t ext2d { { -15.0, 15.0 }, { -10.0, 10.0 } }; + const std::vector res3d { 10, 10, 10 }; + const boundaries_t ext3d { { 0.0, 1.0 }, { 0.0, 1.0 }, { 0.0, 1.0 } }; + + testCustomPrtlUpdate>(res1d, ext1d, {}); + testCustomPrtlUpdate>(res2d, ext2d, {}); + testCustomPrtlUpdate>(res3d, ext3d, {}); + + } catch (std::exception& e) { + std::cerr << e.what() << std::endl; + Kokkos::finalize(); + return 1; + } + Kokkos::finalize(); + return 0; +}