Skip to content

Commit 57f706a

Browse files
LudwigBoessLudwigBoess
committed
added option to define custom particle BCs in pgen
1 parent 8bfa9d6 commit 57f706a

4 files changed

Lines changed: 68 additions & 14 deletions

File tree

‎src/archetypes/traits.h‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,14 @@ namespace arch {
8585
pgen.EmissionPolicy(time, sp, domain);
8686
};
8787

88+
template <class PG, class M, class D>
89+
concept HasCustomPrtlBC = requires(const PG& pgen,
90+
simtime_t time,
91+
spidx_t sp,
92+
D& domain) {
93+
pgen.CustomParticleBoundary(time, sp, domain);
94+
};
95+
8896
template <class PG, class D>
8997
concept HasInitPrtls = requires(PG& pgen, D& domain) {
9098
{ pgen.InitPrtls(domain) } -> std::same_as<void>;

‎src/engines/srpic/particle_pusher.h‎

Lines changed: 22 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -36,18 +36,29 @@ namespace ntt {
3636
const M& metric,
3737
const PG& pgen,
3838
const F& external_fields) {
39+
auto get_custom_prtl_bc = [&]() {
40+
if constexpr (arch::traits::pgen::HasCustomPrtlBC<PG, M, decltype(domain)>) {
41+
return pgen.CustomParticleBoundary(
42+
pusher_params.time, pusher_params.species_index, domain);
43+
} else {
44+
return kernel::sr::NoCustomPrtlBC_t<SimEngine::SRPIC, M> {};
45+
}
46+
};
47+
const auto custom_prtl_bc = get_custom_prtl_bc();
48+
3949
if (emission_policy_flag == EmissionType::NONE) {
4050
const auto no_emission = kernel::NoEmissionPolicy_t<SimEngine::SRPIC, M> {};
4151
Kokkos::parallel_for(
4252
"ParticlePusher",
4353
range,
44-
kernel::sr::Pusher_kernel<M, F, Atm, decltype(no_emission)>(
54+
kernel::sr::Pusher_kernel<M, F, Atm, decltype(no_emission), decltype(custom_prtl_bc)>(
4555
pusher_params,
4656
pusher_arrays,
4757
EB,
4858
metric,
4959
external_fields,
50-
no_emission));
60+
no_emission,
61+
custom_prtl_bc));
5162
} else if (emission_policy_flag == EmissionType::SYNCHROTRON) {
5263
const auto photon_species = params.get<spidx_t>(
5364
"radiation.emission.synchrotron.photon_species");
@@ -77,13 +88,14 @@ namespace ntt {
7788
Kokkos::parallel_for(
7889
"ParticlePusher",
7990
range,
80-
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy)>(
91+
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy), decltype(custom_prtl_bc)>(
8192
pusher_params,
8293
pusher_arrays,
8394
EB,
8495
metric,
8596
external_fields,
86-
emission_policy));
97+
emission_policy,
98+
custom_prtl_bc));
8799
const auto n_inj = emission_policy.numbers_injected();
88100
raise::ErrorIf(n_inj.size() != 1,
89101
"Synchrotron emission should only inject one species",
@@ -120,13 +132,14 @@ namespace ntt {
120132
Kokkos::parallel_for(
121133
"ParticlePusher",
122134
range,
123-
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy)>(
135+
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy), decltype(custom_prtl_bc)>(
124136
pusher_params,
125137
pusher_arrays,
126138
EB,
127139
metric,
128140
external_fields,
129-
emission_policy));
141+
emission_policy,
142+
custom_prtl_bc));
130143
const auto n_inj = emission_policy.numbers_injected();
131144
raise::ErrorIf(n_inj.size() != 1,
132145
"Compton emission should only inject one species",
@@ -148,13 +161,14 @@ namespace ntt {
148161
Kokkos::parallel_for(
149162
"ParticlePusher",
150163
range,
151-
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy)>(
164+
kernel::sr::Pusher_kernel<M, F, Atm, decltype(emission_policy), decltype(custom_prtl_bc)>(
152165
pusher_params,
153166
pusher_arrays,
154167
EB,
155168
metric,
156169
external_fields,
157-
emission_policy));
170+
emission_policy,
171+
custom_prtl_bc));
158172
const auto emitted_species = emission_policy.emitted_species_indices();
159173
const auto n_inj = emission_policy.numbers_injected();
160174
raise::ErrorIf(emitted_species.size() != n_inj.size(),

‎src/kernels/digital_filter.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -93,8 +93,8 @@ namespace kernel {
9393
? (boundaries[2].second == FldsBC::CONDUCTOR)
9494
: false }
9595
, i1_max { size_[0] + N_GHOSTS }
96-
, i2_max { (short)D > 1 ? (size_[1] + N_GHOSTS) : 0 }
97-
, i3_max { (short)D > 2 ? (size_[2] + N_GHOSTS) : 0 } {}
96+
, i2_max { (short)D > 1 ? (size_[(short)D > 1 ? 1 : 0] + N_GHOSTS) : 0 }
97+
, i3_max { (short)D > 2 ? (size_[(short)D > 2 ? 2 : 0] + N_GHOSTS) : 0 } {}
9898

9999
Inline void operator()(index_t i1) const {
100100
if constexpr ((D == Dim::_1D) && (C == Coord::Cart)) {

‎src/kernels/particle_pusher_sr.hpp‎

Lines changed: 36 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -101,11 +101,17 @@ namespace kernel::sr {
101101
array_t<short*> tag;
102102
};
103103

104+
template <SimEngine::type S, class M>
105+
struct NoCustomPrtlBC_t {
106+
template <class PusherKernel>
107+
Inline void operator()(index_t, int, bool, const PusherKernel&) const {}
108+
};
109+
104110
/**
105111
* @tparam M Metric
106112
* @tparam F Additional force
107113
*/
108-
template <class M, class F, bool Atm, class E>
114+
template <class M, class F, bool Atm, class E, class PBC = NoCustomPrtlBC_t<SimEngine::SRPIC, M>>
109115
requires metric::traits::HasD<M> && metric::traits::HasTransformXYZ<M> &&
110116
metric::traits::HasConvertXYZ<M> &&
111117
metric::traits::HasTransform_i<M> && metric::traits::HasConvert_i<M>
@@ -120,12 +126,12 @@ namespace kernel::sr {
120126
std::is_same<E, const NoEmissionPolicy_t<SimEngine::SRPIC, M>>::value,
121127
"Invalid emission policy E for Pusher_kernel");
122128

123-
private:
124129
const ParticlePusherFlags pusher_flags;
125130
const RadiativeDragFlags radiative_drag_flags;
126131

127132
const randacc_ndfield_t<D, 6> EB;
128133

134+
public:
129135
const spidx_t sp;
130136
array_t<int*> i1, i2, i3;
131137
array_t<int*> i1_prev, i2_prev, i3_prev;
@@ -136,14 +142,17 @@ namespace kernel::sr {
136142
array_t<real_t*> weight;
137143
array_t<short*> tag;
138144
const M metric;
145+
const int ni1, ni2, ni3;
146+
147+
private:
139148
const F external_fields;
140149

141150
const E emission_policy;
151+
const PBC custom_prtl_bc;
142152

143153
const simtime_t time;
144154
const real_t coeff, dt;
145155

146-
const int ni1, ni2, ni3;
147156
bool is_absorb_i1min { false }, is_absorb_i1max { false };
148157
bool is_absorb_i2min { false }, is_absorb_i2max { false };
149158
bool is_absorb_i3min { false }, is_absorb_i3max { false };
@@ -153,6 +162,9 @@ namespace kernel::sr {
153162
bool is_reflect_i1min { false }, is_reflect_i1max { false };
154163
bool is_reflect_i2min { false }, is_reflect_i2max { false };
155164
bool is_reflect_i3min { false }, is_reflect_i3max { false };
165+
bool is_custom_i1min { false }, is_custom_i1max { false };
166+
bool is_custom_i2min { false }, is_custom_i2max { false };
167+
bool is_custom_i3min { false }, is_custom_i3max { false };
156168
bool is_axis_i2min { false }, is_axis_i2max { false };
157169

158170
// gca parameters
@@ -169,7 +181,8 @@ namespace kernel::sr {
169181
const randacc_ndfield_t<D, 6>& EB,
170182
const M& metric,
171183
const F& external_fields,
172-
const E& emission_policy)
184+
const E& emission_policy,
185+
const PBC& custom_prtl_bc)
173186
: pusher_flags { pusher_params.pusher_flags }
174187
, radiative_drag_flags { pusher_params.radiative_drag_flags }
175188
, EB { EB }
@@ -195,6 +208,7 @@ namespace kernel::sr {
195208
, metric { metric }
196209
, external_fields { external_fields }
197210
, emission_policy { emission_policy }
211+
, custom_prtl_bc { custom_prtl_bc }
198212
, time { pusher_params.time }
199213
, coeff { HALF * (pusher_params.charge / pusher_params.mass) *
200214
pusher_params.omegaB0 * pusher_params.dt }
@@ -256,6 +270,8 @@ namespace kernel::sr {
256270
is_periodic_i1max = (pusher_params.boundaries[0].second == PrtlBC::PERIODIC);
257271
is_reflect_i1min = (pusher_params.boundaries[0].first == PrtlBC::REFLECT);
258272
is_reflect_i1max = (pusher_params.boundaries[0].second == PrtlBC::REFLECT);
273+
is_custom_i1min = (pusher_params.boundaries[0].first == PrtlBC::CUSTOM);
274+
is_custom_i1max = (pusher_params.boundaries[0].second == PrtlBC::CUSTOM);
259275
if constexpr ((D == Dim::_2D) || (D == Dim::_3D)) {
260276
raise::ErrorIf(pusher_params.boundaries.size() < 2,
261277
"pusher_params.boundaries defined incorrectly",
@@ -269,6 +285,8 @@ namespace kernel::sr {
269285
is_periodic_i2max = (pusher_params.boundaries[1].second == PrtlBC::PERIODIC);
270286
is_reflect_i2min = (pusher_params.boundaries[1].first == PrtlBC::REFLECT);
271287
is_reflect_i2max = (pusher_params.boundaries[1].second == PrtlBC::REFLECT);
288+
is_custom_i2min = (pusher_params.boundaries[1].first == PrtlBC::CUSTOM);
289+
is_custom_i2max = (pusher_params.boundaries[1].second == PrtlBC::CUSTOM);
272290
is_axis_i2min = (pusher_params.boundaries[1].first == PrtlBC::AXIS);
273291
is_axis_i2max = (pusher_params.boundaries[1].second == PrtlBC::AXIS);
274292
}
@@ -285,6 +303,8 @@ namespace kernel::sr {
285303
is_periodic_i3max = (pusher_params.boundaries[2].second == PrtlBC::PERIODIC);
286304
is_reflect_i3min = (pusher_params.boundaries[2].first == PrtlBC::REFLECT);
287305
is_reflect_i3max = (pusher_params.boundaries[2].second == PrtlBC::REFLECT);
306+
is_custom_i3min = (pusher_params.boundaries[2].first == PrtlBC::CUSTOM);
307+
is_custom_i3max = (pusher_params.boundaries[2].second == PrtlBC::CUSTOM);
288308
}
289309
}
290310

@@ -1315,6 +1335,8 @@ namespace kernel::sr {
13151335
i1(p) = 0;
13161336
dx1(p) = ONE - dx1(p);
13171337
invert_vel = true;
1338+
} else if (is_custom_i1min) {
1339+
custom_prtl_bc(p, 1, true, *this);
13181340
}
13191341
} else if (i1(p) >= ni1) {
13201342
if (is_periodic_i1max) {
@@ -1326,6 +1348,8 @@ namespace kernel::sr {
13261348
i1(p) = ni1 - 1;
13271349
dx1(p) = ONE - dx1(p);
13281350
invert_vel = true;
1351+
} else if (is_custom_i1max) {
1352+
custom_prtl_bc(p, 1, false, *this);
13291353
}
13301354
}
13311355
if (invert_vel) {
@@ -1357,6 +1381,8 @@ namespace kernel::sr {
13571381
i2(p) = 0;
13581382
dx2(p) = ONE - dx2(p);
13591383
invert_vel = true;
1384+
} else if (is_custom_i2min) {
1385+
custom_prtl_bc(p, 2, true, *this);
13601386
} else if (is_axis_i2min) {
13611387
i2(p) = 0;
13621388
dx2(p) = ONE - dx2(p);
@@ -1371,6 +1397,8 @@ namespace kernel::sr {
13711397
i2(p) = ni2 - 1;
13721398
dx2(p) = ONE - dx2(p);
13731399
invert_vel = true;
1400+
} else if (is_custom_i2max) {
1401+
custom_prtl_bc(p, 2, false, *this);
13741402
} else if (is_axis_i2max) {
13751403
i2(p) = ni2 - 1;
13761404
dx2(p) = ONE - dx2(p);
@@ -1405,6 +1433,8 @@ namespace kernel::sr {
14051433
i3(p) = 0;
14061434
dx3(p) = ONE - dx3(p);
14071435
invert_vel = true;
1436+
} else if (is_custom_i3min) {
1437+
custom_prtl_bc(p, 3, true, *this);
14081438
}
14091439
} else if (i3(p) >= ni3) {
14101440
if (is_periodic_i3max) {
@@ -1416,6 +1446,8 @@ namespace kernel::sr {
14161446
i3(p) = ni3 - 1;
14171447
dx3(p) = ONE - dx3(p);
14181448
invert_vel = true;
1449+
} else if (is_custom_i3max) {
1450+
custom_prtl_bc(p, 3, false, *this);
14191451
}
14201452
}
14211453
if (invert_vel) {

0 commit comments

Comments
 (0)