Skip to content

Commit e9a124a

Browse files
StaticObserverStaticObserver
committed
injection with weights
1 parent cd1470a commit e9a124a

3 files changed

Lines changed: 420 additions & 8 deletions

File tree

‎pgens/BZjet/pgen.hpp‎

Lines changed: 26 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -189,6 +189,7 @@ namespace user {
189189
const real_t sigma_thr,
190190
const real_t inj_coeff,
191191
const real_t db_thr,
192+
const bool is_weight,
192193
const SimulationParams& params,
193194
Domain<S, M>* domain_ptr)
194195
: arch::SpatialDistribution<S, M> { domain_ptr->mesh.metric }
@@ -200,7 +201,9 @@ namespace user {
200201
, d0 { params.template get<real_t>("scales.skindepth0") }
201202
, rho0 { params.template get<real_t>("scales.larmor0") }
202203
, inv_n0 { ONE / params.template get<real_t>("scales.n0") }
203-
, inj_coeff { inj_coeff } {
204+
, inv_ppc0 { ONE / params.template get<real_t>("particles.ppc0") }
205+
, inj_coeff { inj_coeff }
206+
, is_weight { is_weight } {
204207
std::copy(xi_min.begin(), xi_min.end(), x_min);
205208
std::copy(xi_max.begin(), xi_max.end(), x_max);
206209

@@ -283,7 +286,11 @@ namespace user {
283286
DOT(B_cntrv[0], B_cntrv[1], B_cntrv[2], B_cov[0], B_cov[1], B_cov[2]);
284287
const auto db = DOT(D_cntrv[0], D_cntrv[1], D_cntrv[2], B_cov[0], B_cov[1], B_cov[2]);
285288
const real_t inj_n = inj_coeff * db * SIGN(db) / math::sqrt(bsqr) * SQR(d0) / rho0;
286-
return fill ? inj_n : ZERO;
289+
if (is_weight) {
290+
return fill ? inj_n : ZERO;
291+
} else {
292+
return fill ? TWO * inv_ppc0 * 1.01 : ZERO;
293+
}
287294
}
288295

289296
private:
@@ -295,10 +302,12 @@ namespace user {
295302
const real_t inv_n0;
296303
const real_t d0;
297304
const real_t rho0;
305+
const real_t inv_ppc0;
298306
Domain<S, M>* domain_ptr;
299307
ndfield_t<M::Dim, 3> density;
300308
ndfield_t<M::Dim, 6> EM;
301309
const M metric;
310+
const bool is_weight;
302311
};
303312

304313

@@ -341,24 +350,33 @@ namespace user {
341350
const auto energy_dist = arch::Maxwellian<S, M>(local_domain.mesh.metric,
342351
local_domain.random_pool,
343352
temperature);
344-
const auto spatial_dist = PointDistribution<S, M>(xi_min,
353+
const auto spatial_dist1 = PointDistribution<S, M>(xi_min,
354+
xi_max,
355+
sigma_max / sigma0,
356+
inj_coeff,
357+
db_thr,
358+
false,
359+
params,
360+
&local_domain);
361+
const auto spatial_dist2 = PointDistribution<S, M>(xi_min,
345362
xi_max,
346363
sigma_max / sigma0,
347364
inj_coeff,
348365
db_thr,
366+
true,
349367
params,
350368
&local_domain);
351369

352370
const auto injector =
353-
arch::NonUniformInjector<S, M, arch::Maxwellian, PointDistribution>(
371+
arch::experimental::Injector_with_weights<S, M, arch::Maxwellian, PointDistribution, PointDistribution>(
354372
energy_dist,
355-
spatial_dist,
373+
spatial_dist1,
374+
spatial_dist2,
356375
{ 1, 2 });
357-
arch::InjectNonUniform<S, M, decltype(injector)>(params,
376+
arch::experimental::InjectWithWeights<S, M, decltype(injector)>(params,
358377
local_domain,
359378
injector,
360-
1.0,
361-
true);
379+
1.0);
362380
}
363381

364382
void CustomFieldOutput(const std::string& name,

‎src/archetypes/particle_injector.h‎

Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -284,6 +284,8 @@ namespace arch {
284284
~NonUniformInjector() = default;
285285
};
286286

287+
288+
287289
template <SimEngine::type S, class M, bool P, in O>
288290
struct AtmosphereInjector {
289291
struct TargetDensityProfile {
@@ -734,6 +736,112 @@ namespace arch {
734736
}
735737
}
736738

739+
namespace experimental {
740+
template <SimEngine::type S,
741+
class M,
742+
template <SimEngine::type, class> class ED,
743+
template <SimEngine::type, class> class SD1,
744+
template <SimEngine::type, class> class SD2>
745+
struct Injector_with_weights {
746+
using energy_dist_t = ED<S, M>;
747+
using spatial_dist_t1 = SD1<S, M>;
748+
using spatial_dist_t2 = SD2<S, M>;
749+
static_assert(M::is_metric, "M must be a metric class");
750+
static_assert(energy_dist_t::is_energy_dist,
751+
"E must be an energy distribution class");
752+
static_assert(spatial_dist_t1::is_spatial_dist,
753+
"SD must be a spatial distribution class");
754+
static_assert(spatial_dist_t2::is_spatial_dist,
755+
"SD must be a spatial distribution class");
756+
static constexpr bool is_nonuniform_injector { true };
757+
static constexpr Dimension D { M::Dim };
758+
static constexpr Coord C { M::CoordType };
759+
760+
const energy_dist_t energy_dist;
761+
const spatial_dist_t1 spatial_dist1;
762+
const spatial_dist_t2 spatial_dist2;
763+
const std::pair<spidx_t, spidx_t> species;
764+
765+
Injector_with_weights(const energy_dist_t& energy_dist,
766+
const spatial_dist_t1& spatial_dist1,
767+
const spatial_dist_t2& spatial_dist2,
768+
const std::pair<spidx_t, spidx_t>& species)
769+
: energy_dist { energy_dist }
770+
, spatial_dist1 { spatial_dist1 }
771+
, spatial_dist2 { spatial_dist2 }
772+
, species { species } {}
773+
774+
~Injector_with_weights() = default;
775+
}; // struct Injector_with_weights
776+
777+
template <SimEngine::type S, class M, class I>
778+
inline void InjectWithWeights(const SimulationParams& params,
779+
Domain<S, M>& domain,
780+
const I& injector,
781+
real_t number_density,
782+
const boundaries_t<real_t>& box = {}) {
783+
static_assert(M::is_metric, "M must be a metric class");
784+
static_assert(I::is_nonuniform_injector,
785+
"I must be a nonuniform injector class");
786+
raise::ErrorIf(
787+
not params.template get<bool>("particles.use_weights"),
788+
"Weights must be enabled for the injector",
789+
HERE);
790+
if (domain.species[injector.species.first - 1].charge() +
791+
domain.species[injector.species.second - 1].charge() !=
792+
0.0f) {
793+
raise::Warning("Total charge of the injected species is non-zero", HERE);
794+
}
795+
{
796+
range_t<M::Dim> cell_range;
797+
if (box.size() == 0) {
798+
cell_range = domain.mesh.rangeActiveCells();
799+
} else {
800+
raise::ErrorIf(box.size() != M::Dim,
801+
"Box must have the same dimension as the mesh",
802+
HERE);
803+
boundaries_t<bool> incl_ghosts;
804+
for (auto d = 0; d < M::Dim; ++d) {
805+
incl_ghosts.push_back({ false, false });
806+
}
807+
const auto extent = domain.mesh.ExtentToRange(box, incl_ghosts);
808+
tuple_t<ncells_t, M::Dim> x_min { 0 }, x_max { 0 };
809+
for (auto d = 0; d < M::Dim; ++d) {
810+
x_min[d] = extent[d].first;
811+
x_max[d] = extent[d].second;
812+
}
813+
cell_range = CreateRangePolicy<M::Dim>(x_min, x_max);
814+
}
815+
const auto ppc = number_density *
816+
params.template get<real_t>("particles.ppc0") * HALF;
817+
auto injector_kernel =
818+
kernel::experimental::Injector_with_weights_kernel<S, M, typename I::energy_dist_t, typename I::spatial_dist_t1, typename I::spatial_dist_t2>(
819+
ppc,
820+
injector.species.first,
821+
injector.species.second,
822+
domain.species[injector.species.first - 1],
823+
domain.species[injector.species.second - 1],
824+
domain.species[injector.species.first - 1].npart(),
825+
domain.species[injector.species.second - 1].npart(),
826+
domain.mesh.metric,
827+
injector.energy_dist,
828+
injector.spatial_dist1,
829+
injector.spatial_dist2,
830+
ONE / params.template get<real_t>("scales.V0"),
831+
domain.random_pool);
832+
Kokkos::parallel_for("InjectWithWeights",
833+
cell_range,
834+
injector_kernel);
835+
const auto n_inj = injector_kernel.number_injected();
836+
domain.species[injector.species.first - 1].set_npart(
837+
domain.species[injector.species.first - 1].npart() + n_inj);
838+
domain.species[injector.species.second - 1].set_npart(
839+
domain.species[injector.species.second - 1].npart() + n_inj);
840+
}
841+
} // function InjectWithWeights
842+
843+
} // namespace experimental
844+
737845
} // namespace arch
738846

739847
#endif // ARCHETYPES_PARTICLE_INJECTOR_H

0 commit comments

Comments
 (0)