Skip to content

Custom particle BCs #190

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
LudwigBoess wants to merge 2 commits into
base: 1.4.0rc
Choose a base branch
from
Draft

Custom particle BCs #190

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
57f706a
added option to define custom particle BCs in pgen
LudwigBoess Mar 21, 2026
4a5bca1
Merge branch '1.4.0rc' into dev/custom_particle_bcs
LudwigBoess Mar 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions src/archetypes/traits.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,14 @@ namespace arch {
pgen.EmissionPolicy(time, sp, domain);
};

template <class PG, class M, class D>
concept HasCustomPrtlBC = requires(const PG& pgen,
simtime_t time,
spidx_t sp,
D& domain) {
pgen.CustomParticleBoundary(time, sp, domain);
};

template <class PG, class D>
concept HasInitPrtls = requires(PG& pgen, D& domain) {
{ pgen.InitPrtls(domain) } -> std::same_as<void>;
Expand Down
51 changes: 32 additions & 19 deletions src/engines/srpic/particle_pusher.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,29 +26,39 @@ namespace ntt {
namespace srpic {

template <class M, class F, class PG, bool Atm>
void CallPusher_WithExternalFieldFlag(
Domain<SimEngine::SRPIC, M>& domain,
const SimulationParams& params,
const kernel::sr::PusherParams& pusher_params,
kernel::sr::PusherArrays& pusher_arrays,
EmissionTypeFlag emission_policy_flag,
const range_t<Dim::_1D>& range,
const ndfield_t<M::Dim, 6>& EB,
const M& metric,
const PG& pgen,
const F& external_fields) {
void CallPusher(Domain<SimEngine::SRPIC, M>& domain,
const SimulationParams& params,
const kernel::sr::PusherParams& pusher_params,
kernel::sr::PusherArrays& pusher_arrays,
EmissionTypeFlag emission_policy_flag,
const range_t<Dim::_1D>& range,
const ndfield_t<M::Dim, 6>& EB,
const M& metric,
const PG& pgen,
const F& external_fields) {
auto get_custom_prtl_bc = [&]() {
if constexpr (arch::traits::pgen::HasCustomPrtlBC<PG, M, decltype(domain)>) {
return pgen.CustomParticleBoundary(
pusher_params.time, pusher_params.species_index, domain);
} else {
return kernel::sr::NoCustomPrtlBC_t<SimEngine::SRPIC, M> {};
}
};
const auto custom_prtl_bc = get_custom_prtl_bc();

if (emission_policy_flag == EmissionType::NONE) {
const auto no_emission = kernel::NoEmissionPolicy_t<SimEngine::SRPIC, M> {};
Kokkos::parallel_for(
"ParticlePusher",
range,
kernel::sr::Pusher_kernel<M, F, Atm, decltype(no_emission)>(
kernel::sr::Pusher_kernel<M, F, Atm, decltype(no_emission), decltype(custom_prtl_bc)>(
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<spidx_t>(
"radiation.emission.synchrotron.photon_species");
Expand Down Expand Up @@ -78,13 +88,14 @@ namespace ntt {
Kokkos::parallel_for(
"ParticlePusher",
range,
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy)>(
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy), decltype(custom_prtl_bc)>(
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",
Expand Down Expand Up @@ -121,13 +132,14 @@ namespace ntt {
Kokkos::parallel_for(
"ParticlePusher",
range,
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy)>(
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy), decltype(custom_prtl_bc)>(
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",
Expand All @@ -148,13 +160,14 @@ namespace ntt {
Kokkos::parallel_for(
"ParticlePusher",
range,
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy)>(
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy), decltype(custom_prtl_bc)>(
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(),
Expand Down
4 changes: 2 additions & 2 deletions src/kernels/digital_filter.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)) {
Expand Down
42 changes: 38 additions & 4 deletions src/kernels/particle_pusher_sr.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,11 +101,17 @@ namespace kernel::sr {
array_t<short*> tag;
};

template <SimEngine::type S, class M>
struct NoCustomPrtlBC_t {
template <class PusherKernel>
Inline void operator()(index_t, int, bool, const PusherKernel&) const {}
};

/**
* @tparam M Metric
* @tparam F Additional force
*/
template <class M, class F, bool Atm, class E>
template <class M, class F, bool Atm, class E, class PBC = NoCustomPrtlBC_t<SimEngine::SRPIC, M>>
requires metric::traits::HasD<M> && metric::traits::HasTransformXYZ<M> &&
metric::traits::HasConvertXYZ<M> &&
metric::traits::HasTransform_i<M> && metric::traits::HasConvert_i<M>
Expand All @@ -120,12 +126,13 @@ namespace kernel::sr {
std::is_same<E, const NoEmissionPolicy_t<SimEngine::SRPIC, M>>::value,
"Invalid emission policy E for Pusher_kernel");

private:
const ParticlePusherFlags pusher_flags;
const RadiativeDragFlags radiative_drag_flags;

const randacc_ndfield_t<D, 6> EB;

public:
const spidx_t sp;
array_t<int*> i1, i2, i3;
array_t<int*> i1_prev, i2_prev, i3_prev;
array_t<prtldx_t*> dx1, dx2, dx3;
Expand All @@ -135,13 +142,16 @@ namespace kernel::sr {
array_t<real_t*> weight;
array_t<short*> tag;
const M metric;
const int ni1, ni2, ni3;

private:
const F external_fields;

const E emission_policy;
const PBC custom_prtl_bc;

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 };
Expand All @@ -151,6 +161,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
Expand All @@ -167,7 +180,8 @@ namespace kernel::sr {
const randacc_ndfield_t<D, 6>& 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 }
Expand All @@ -192,6 +206,8 @@ 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 }
, dt { pusher_params.dt }
Expand Down Expand Up @@ -252,6 +268,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",
Expand All @@ -265,6 +283,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);
}
Expand All @@ -281,6 +301,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);
}
}

Expand Down Expand Up @@ -1352,6 +1374,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) {
Expand All @@ -1363,6 +1387,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) {
Expand Down Expand Up @@ -1394,6 +1420,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);
Expand All @@ -1408,6 +1436,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);
Expand Down Expand Up @@ -1442,6 +1472,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) {
Expand All @@ -1453,6 +1485,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) {
Expand Down
Loading