From 341f672056af4c917d70c40393320a6ba6dc4d80 Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 22 Sep 2025 15:24:09 -0400 Subject: [PATCH 1/8] more proper shock injector --- pgens/shock/pgen.hpp | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/pgens/shock/pgen.hpp b/pgens/shock/pgen.hpp index f93bbfd0..87fb7e2f 100644 --- a/pgens/shock/pgen.hpp +++ b/pgens/shock/pgen.hpp @@ -82,19 +82,24 @@ namespace user { using arch::ProblemGenerator::C; using arch::ProblemGenerator::params; + Metadomain& global_domain; + // domain properties - const real_t global_xmin, global_xmax; + const real_t global_xmin, global_xmax; // gas properties - const real_t drift_ux, temperature, temperature_ratio, filling_fraction; + const real_t drift_ux, temperature, temperature_ratio, filling_fraction; // injector properties - const real_t injector_velocity, injection_start, dt; - const int injection_frequency; + const real_t injector_velocity; + const simtime_t injection_start; + const real_t dt; + const int injection_frequency; // magnetic field properties - real_t Btheta, Bphi, Bmag; - InitFields init_flds; + real_t Btheta, Bphi, Bmag; + InitFields init_flds; - inline PGen(const SimulationParams& p, const Metadomain& global_domain) + inline PGen(const SimulationParams& p, 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 } , drift_ux { p.template get("setup.drift_ux") } @@ -106,7 +111,8 @@ namespace user { , init_flds { Bmag, Btheta, Bphi, drift_ux } , filling_fraction { p.template get("setup.filling_fraction", 1.0) } , injector_velocity { p.template get("setup.injector_velocity", 1.0) } - , injection_start { p.template get("setup.injection_start", 0.0) } + , injection_start { static_cast( + p.template get("setup.injection_start", 0.0)) } , injection_frequency { p.template get("setup.injection_frequency", 100) } , dt { p.template get("algorithms.timestep.dt") } {} @@ -206,7 +212,7 @@ namespace user { */ // check if the injector should be active - if (step % injection_frequency != 0) { + if (step % injection_frequency != 0 and time < injection_start) { return; } @@ -222,7 +228,7 @@ namespace user { } // compute the beginning of the injected region - auto xmin = xmax - injection_frequency * dt; + auto xmin = xmax - injector_velocity * injection_frequency * dt; if (xmin <= global_xmin) { xmin = global_xmin; } @@ -258,6 +264,8 @@ namespace user { init_flds, domain.mesh.metric }); + global_domain.CommunicateFields(domain, Comm::B | Comm::E); + /* tag particles inside the injection zone as dead */ From 8ba317a53e8606ef8fdf7bec9b6061247ddcf7f4 Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 23 Sep 2025 23:11:04 +0000 Subject: [PATCH 2/8] test of shock injector --- pgens/shock/pgen.hpp | 198 +++++++++++++------- src/archetypes/particle_injector.h | 218 ++++++++++++++++------ src/archetypes/spatial_dist.h | 75 +++++++- src/kernels/injectors.hpp | 287 +++++++++++++++++++++++++++++ 4 files changed, 643 insertions(+), 135 deletions(-) diff --git a/pgens/shock/pgen.hpp b/pgens/shock/pgen.hpp index 87fb7e2f..16f5197b 100644 --- a/pgens/shock/pgen.hpp +++ b/pgens/shock/pgen.hpp @@ -118,12 +118,12 @@ namespace user { inline PGen() {} - auto MatchFields(real_t time) const -> InitFields { + auto MatchFields(real_t) const -> InitFields { return init_flds; } - auto FixFieldsConst(const bc_in&, const em& comp) const - -> std::pair { + auto FixFieldsConst(const bc_in&, + const em& comp) const -> std::pair { if (comp == em::ex1) { return { init_flds.ex1({ ZERO }), true }; } else if (comp == em::ex2) { @@ -228,62 +228,122 @@ namespace user { } // compute the beginning of the injected region - auto xmin = xmax - injector_velocity * injection_frequency * dt; + auto xmin = xmax - 1.1 * injection_frequency * dt; if (xmin <= global_xmin) { xmin = global_xmin; } - // define indice range to reset fields - boundaries_t incl_ghosts; - for (auto d = 0; d < M::Dim; ++d) { - incl_ghosts.push_back({ false, false }); - } + // // define indice range to reset fields + // boundaries_t incl_ghosts; + // for (auto d = 0; d < M::Dim; ++d) { + // incl_ghosts.push_back({ false, false }); + // } + + // // define box to reset fields + // boundaries_t purge_box; + // // loop over all dimension + // for (auto d = 0u; d < M::Dim; ++d) { + // if (d == 0) { + // purge_box.push_back({ xmin, global_xmax }); + // } else { + // purge_box.push_back(Range::All); + // } + // } + // + // const auto extent = domain.mesh.ExtentToRange(purge_box, incl_ghosts); + // tuple_t x_min { 0 }, x_max { 0 }; + // for (auto d = 0; d < M::Dim; ++d) { + // x_min[d] = extent[d].first; + // x_max[d] = extent[d].second; + // } + // + // Kokkos::parallel_for("ResetFields", + // CreateRangePolicy(x_min, x_max), + // arch::SetEMFields_kernel { + // domain.fields.em, + // init_flds, + // domain.mesh.metric }); + // + // global_domain.CommunicateFields(domain, Comm::B | Comm::E); + + // /* + // tag particles inside the injection zone as dead + // */ + // const auto& mesh = domain.mesh; + // + // // loop over particle species + // for (auto s { 0u }; s < 2; ++s) { + // // get particle properties + // auto& species = domain.species[s]; + // auto i1 = species.i1; + // auto dx1 = species.dx1; + // auto tag = species.tag; + // + // Kokkos::parallel_for( + // "RemoveParticles", + // species.rangeActiveParticles(), + // Lambda(index_t p) { + // // check if the particle is already dead + // if (tag(p) == ParticleTag::dead) { + // return; + // } + // const auto x_Cd = static_cast(i1(p)) + + // static_cast(dx1(p)); + // const auto x_Ph = mesh.metric.template convert<1, Crd::Cd, Crd::XYZ>( + // x_Cd); + // + // if (x_Ph > xmin) { + // tag(p) = ParticleTag::dead; + // } + // }); + // } + + /* + Inject slab of fresh plasma + */ - // define box to reset fields - boundaries_t purge_box; + // define box to inject into + boundaries_t inj_box; // loop over all dimension for (auto d = 0u; d < M::Dim; ++d) { if (d == 0) { - purge_box.push_back({ xmin, global_xmax }); + inj_box.push_back({ xmin, xmax }); } else { - purge_box.push_back(Range::All); + inj_box.push_back(Range::All); } } - const auto extent = domain.mesh.ExtentToRange(purge_box, incl_ghosts); - tuple_t x_min { 0 }, x_max { 0 }; - for (auto d = 0; d < M::Dim; ++d) { - x_min[d] = extent[d].first; - x_max[d] = extent[d].second; - } - - Kokkos::parallel_for("ResetFields", - CreateRangePolicy(x_min, x_max), - arch::SetEMFields_kernel { - domain.fields.em, - init_flds, - domain.mesh.metric }); - - global_domain.CommunicateFields(domain, Comm::B | Comm::E); + // same maxwell distribution as above + const auto mass_1 = domain.species[0].mass(); + const auto mass_2 = domain.species[1].mass(); + const auto temperature_1 = temperature / mass_1; + const auto temperature_2 = temperature * temperature_ratio / mass_2; + + const auto maxwellian_1 = arch::experimental::Maxwellian( + domain.mesh.metric, + domain.random_pool, + temperature_1, + { -drift_ux, ZERO, ZERO }); + const auto maxwellian_2 = arch::experimental::Maxwellian( + domain.mesh.metric, + domain.random_pool, + temperature_2, + { -drift_ux, ZERO, ZERO }); - /* - tag particles inside the injection zone as dead - */ const auto& mesh = domain.mesh; - - // loop over particle species - for (auto s { 0u }; s < 2; ++s) { - // get particle properties - auto& species = domain.species[s]; + for (auto& s : { 1u, 2u }) { + auto& species = domain.species[s - 1]; auto i1 = species.i1; auto dx1 = species.dx1; + auto ux1 = species.ux1; + auto ux2 = species.ux2; + auto ux3 = species.ux3; auto tag = species.tag; Kokkos::parallel_for( - "RemoveParticles", + "ResetParticles", species.rangeActiveParticles(), Lambda(index_t p) { - // check if the particle is already dead if (tag(p) == ParticleTag::dead) { return; } @@ -292,42 +352,42 @@ namespace user { const auto x_Ph = mesh.metric.template convert<1, Crd::Cd, Crd::XYZ>( x_Cd); - if (x_Ph > xmin) { - tag(p) = ParticleTag::dead; + const coord_t x_dummy { ZERO }; + if (x_Ph > xmin and x_Ph <= xmax) { + vec_t v_T { ZERO }, v_Cd { ZERO }; + if (s == 1u) { + maxwellian_1(x_dummy, v_T, s); + } else { + maxwellian_2(x_dummy, v_T, s); + } + mesh.metric.template transform_xyz(x_dummy, v_T, v_Cd); + ux1(p) = v_Cd[0]; + ux2(p) = v_Cd[1]; + ux3(p) = v_Cd[2]; + } else if (x_Ph > xmax) { + if (injector_velocity < 0.9) { + ux1(p) = TWO * injector_velocity / + math::sqrt(ONE - SQR(injector_velocity)); + } else { + ux1(p) = TWO * injector_velocity / + math::sqrt(ONE - SQR(injector_velocity) + 0.001); + } + ux2(p) = ZERO; + ux3(p) = ZERO; } }); } - /* - Inject slab of fresh plasma - */ - - // define box to inject into - boundaries_t inj_box; - // loop over all dimension - for (auto d = 0u; d < M::Dim; ++d) { - if (d == 0) { - inj_box.push_back({ xmin, xmax }); - } else { - inj_box.push_back(Range::All); - } - } - - // same maxwell distribution as above - const auto temperatures = std::make_pair(temperature, - temperature_ratio * temperature); - const auto drifts = std::make_pair( - std::vector { -drift_ux, ZERO, ZERO }, - std::vector { -drift_ux, ZERO, ZERO }); - arch::InjectUniformMaxwellians(params, - domain, - ONE, - temperatures, - { 1, 2 }, - drifts, - false, - inj_box); + arch::InjectReplenishConst( + params, + domain, + ONE, + { 1, 2 }, + maxwellian_1, + maxwellian_2, + inj_box); } }; + } // namespace user #endif diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h index 6313031d..4ab3db6a 100644 --- a/src/archetypes/particle_injector.h +++ b/src/archetypes/particle_injector.h @@ -364,67 +364,67 @@ namespace arch { ~AtmosphereInjector() = default; }; - template - struct MovingInjector { - struct TargetDensityProfile { - const real_t nmax, xinj, xdrift; - - TargetDensityProfile(real_t xinj, real_t xdrift, real_t nmax) - : xinj { xinj } - , xdrift { xdrift } - , nmax { nmax } {} - - Inline auto operator()(const coord_t& x_Ph) const -> real_t { - if constexpr ((O == in::x1) or - (O == in::x2 and (M::Dim == Dim::_2D or M::Dim == Dim::_3D)) or - (O == in::x3 and M::Dim == Dim::_3D)) { - const auto xi = x_Ph[static_cast(O)]; - // + direction - if (xi < xdrift or xi >= xinj) { - return ZERO; - } else { - if constexpr (M::CoordType == Coord::Cart) { - return nmax; - } else { - raise::KernelError( - HERE, - "Moving injector in +x cannot be applied for non-cartesian"); - return ZERO; - } - } - } else { - raise::KernelError(HERE, "Wrong direction"); - return ZERO; - } - } - }; - - using energy_dist_t = Maxwellian; - using spatial_dist_t = Replenish; - static_assert(M::is_metric, "M must be a metric class"); - static constexpr bool is_nonuniform_injector { true }; - static constexpr Dimension D { M::Dim }; - static constexpr Coord C { M::CoordType }; - - const energy_dist_t energy_dist; - const TargetDensityProfile target_density; - const spatial_dist_t spatial_dist; - const std::pair species; - - MovingInjector(const M& metric, - const ndfield_t& density, - const energy_dist_t& energy_dist, - real_t xinj, - real_t xdrift, - real_t nmax, - const std::pair& species) - : energy_dist { energy_dist } - , target_density { xinj, xdrift, nmax } - , spatial_dist { metric, density, 0, target_density, nmax } - , species { species } {} - - ~MovingInjector() = default; - }; + // template + // struct MovingInjector { + // struct TargetDensityProfile { + // const real_t nmax, xinj, xdrift; + // + // TargetDensityProfile(real_t xinj, real_t xdrift, real_t nmax) + // : xinj { xinj } + // , xdrift { xdrift } + // , nmax { nmax } {} + // + // Inline auto operator()(const coord_t& x_Ph) const -> real_t { + // if constexpr ((O == in::x1) or + // (O == in::x2 and (M::Dim == Dim::_2D or M::Dim == + // Dim::_3D)) or (O == in::x3 and M::Dim == Dim::_3D)) { + // const auto xi = x_Ph[static_cast(O)]; + // // + direction + // if (xi < xdrift or xi >= xinj) { + // return ZERO; + // } else { + // if constexpr (M::CoordType == Coord::Cart) { + // return nmax; + // } else { + // raise::KernelError( + // HERE, + // "Moving injector in +x cannot be applied for non-cartesian"); + // return ZERO; + // } + // } + // } else { + // raise::KernelError(HERE, "Wrong direction"); + // return ZERO; + // } + // } + // }; + // + // using energy_dist_t = Maxwellian; + // using spatial_dist_t = Replenish; + // static_assert(M::is_metric, "M must be a metric class"); + // static constexpr bool is_nonuniform_injector { true }; + // static constexpr Dimension D { M::Dim }; + // static constexpr Coord C { M::CoordType }; + // + // const energy_dist_t energy_dist; + // const TargetDensityProfile target_density; + // const spatial_dist_t spatial_dist; + // const std::pair species; + // + // MovingInjector(const M& metric, + // const ndfield_t& density, + // const energy_dist_t& energy_dist, + // real_t xinj, + // real_t xdrift, + // real_t nmax, + // const std::pair& species) + // : energy_dist { energy_dist } + // , target_density { xinj, xdrift, nmax } + // , spatial_dist { metric, density, 0, target_density, nmax } + // , species { species } {} + // + // ~MovingInjector() = default; + // }; /** * @brief Injects uniform number density of particles everywhere in the domain @@ -738,6 +738,102 @@ namespace arch { } } + template + void InjectReplenishConst(const SimulationParams& params, + Domain& domain, + real_t target_tot_ndens, + const std::pair& species, + const ED1& energy_dist_1, + const ED2& energy_dist_2, + const boundaries_t& box = {}) { + static_assert(M::is_metric, "M must be a metric class"); + static_assert(ED1::is_energy_dist, "ED1 must be an energy distribution class"); + static_assert(ED2::is_energy_dist, "ED2 must be an energy distribution class"); + + range_t cell_range; + if (box.size() == 0) { + cell_range = domain.mesh.rangeActiveCells(); + } else { + raise::ErrorIf(box.size() != M::Dim, + "Box must have the same dimension as the mesh", + HERE); + boundaries_t incl_ghosts; + for (auto d = 0; d < M::Dim; ++d) { + incl_ghosts.push_back({ false, false }); + } + const auto extent = domain.mesh.ExtentToRange(box, incl_ghosts); + tuple_t x_min { 0 }, x_max { 0 }; + for (auto d = 0; d < M::Dim; ++d) { + x_min[d] = extent[d].first; + x_max[d] = extent[d].second; + } + cell_range = CreateRangePolicy(x_min, x_max); + } + + const idx_t buff_idx = 0u; + { + // compute total number density of the species + const auto use_weights = params.template get( + "particles.use_weights"); + const auto ni2 = domain.mesh.n_active(in::x2); + const auto inv_n0 = ONE / params.template get("scales.n0"); + + auto scatter_buff = Kokkos::Experimental::create_scatter_view( + domain.fields.buff); + Kokkos::deep_copy(domain.fields.buff, ZERO); + for (const auto sp : std::vector { species.first, species.second }) { + const auto& prtl_spec = domain.species[sp - 1]; + // clang-format off + Kokkos::parallel_for( + "ComputeMoments", + prtl_spec.rangeActiveParticles(), + kernel::ParticleMoments_kernel({}, scatter_buff, buff_idx, + prtl_spec.i1, prtl_spec.i2, prtl_spec.i3, + prtl_spec.dx1, prtl_spec.dx2, prtl_spec.dx3, + prtl_spec.ux1, prtl_spec.ux2, prtl_spec.ux3, + prtl_spec.phi, prtl_spec.weight, prtl_spec.tag, + prtl_spec.mass(), prtl_spec.charge(), + use_weights, + domain.mesh.metric, domain.mesh.flds_bc(), + ni2, inv_n0, 0u)); + // clang-format on + } + Kokkos::Experimental::contribute(domain.fields.buff, scatter_buff); + } + const auto spatial_dist = ReplenishConst { domain.mesh.metric, + domain.fields.buff, + buff_idx, + target_tot_ndens, + 1.0 }; + + const auto ppc = target_tot_ndens * + params.template get("particles.ppc0") * HALF; + + auto injector_kernel = + kernel::experimental::NonUniformInjector_kernel( + ppc, + species.first, + species.second, + domain.species[species.first - 1], + domain.species[species.second - 1], + domain.species[species.first - 1].npart(), + domain.species[species.second - 1].npart(), + domain.mesh.metric, + energy_dist_1, + energy_dist_2, + spatial_dist, + ONE / params.template get("scales.V0"), + domain.random_pool); + + Kokkos::parallel_for("InjectReplenishConst", cell_range, injector_kernel); + + const auto n_inj = injector_kernel.number_injected(); + domain.species[species.first - 1].set_npart( + domain.species[species.first - 1].npart() + n_inj); + domain.species[species.second - 1].set_npart( + domain.species[species.second - 1].npart() + n_inj); + } + } // namespace arch #endif // ARCHETYPES_PARTICLE_INJECTOR_H diff --git a/src/archetypes/spatial_dist.h b/src/archetypes/spatial_dist.h index 55c84ddf..676ab035 100644 --- a/src/archetypes/spatial_dist.h +++ b/src/archetypes/spatial_dist.h @@ -5,6 +5,7 @@ * - arch::SpatialDistribution<> * - arch::Uniform<> : arch::SpatialDistribution<> * - arch::Replenish<> : arch::SpatialDistribution<> + * - arch::ReplenishConst<> : arch::SpatialDistribution<> * @namespace * - arch:: * @note @@ -45,7 +46,7 @@ namespace arch { } }; - template + template struct Replenish : public SpatialDistribution { using SpatialDistribution::metric; const ndfield_t density; @@ -54,16 +55,28 @@ namespace arch { const T target_density; const real_t target_max_density; + const real_t fill_ratio; + Replenish(const M& metric, - const ndfield_t& density, + const ndfield_t& density, idx_t idx, const T& target_density, - real_t target_max_density) + real_t target_max_density, + real_t fill_ratio = 0.9) : SpatialDistribution { metric } , density { density } , idx { idx } , target_density { target_density } - , target_max_density { target_max_density } {} + , target_max_density { target_max_density } + , fill_ratio { fill_ratio } { + raise::ErrorIf(target_max_density <= ZERO, + "Target max density must be positive", + HERE); + raise::ErrorIf(fill_ratio <= ZERO or fill_ratio > ONE, + "Fill ratio must be in (0, 1]", + HERE); + raise::ErrorIf(idx >= N, "Index for density field is out of bounds", HERE); + } Inline auto operator()(const coord_t& x_Ph) const -> real_t { coord_t x_Cd { ZERO }; @@ -84,7 +97,7 @@ namespace arch { raise::KernelError(HERE, "Invalid dimension"); } const auto target = target_density(x_Ph); - if (0.9 * target > dens) { + if (fill_ratio * target > dens) { return (target - dens) / target_max_density; } else { return ZERO; @@ -92,6 +105,58 @@ namespace arch { } }; + template + struct ReplenishConst : public SpatialDistribution { + using SpatialDistribution::metric; + const ndfield_t density; + const idx_t idx; + + const real_t target_density; + const real_t fill_ratio; + + ReplenishConst(const M& metric, + const ndfield_t& density, + idx_t idx, + real_t target_density, + real_t fill_ratio = 0.9) + : SpatialDistribution { metric } + , density { density } + , idx { idx } + , target_density { target_density } + , fill_ratio { fill_ratio } { + raise::ErrorIf(target_density <= ZERO, "Target density must be positive", HERE); + raise::ErrorIf(fill_ratio <= ZERO or fill_ratio > ONE, + "Fill ratio must be in (0, 1]", + HERE); + raise::ErrorIf(idx >= N, "Index for density field is out of bounds", HERE); + } + + Inline auto operator()(const coord_t& x_Ph) const -> real_t { + coord_t x_Cd { ZERO }; + metric.template convert(x_Ph, x_Cd); + real_t dens { ZERO }; + if constexpr (M::Dim == Dim::_1D) { + dens = density(static_cast(x_Cd[0]) + N_GHOSTS, idx); + } else if constexpr (M::Dim == Dim::_2D) { + dens = density(static_cast(x_Cd[0]) + N_GHOSTS, + static_cast(x_Cd[1]) + N_GHOSTS, + idx); + } else if constexpr (M::Dim == Dim::_3D) { + dens = density(static_cast(x_Cd[0]) + N_GHOSTS, + static_cast(x_Cd[1]) + N_GHOSTS, + static_cast(x_Cd[2]) + N_GHOSTS, + idx); + } else { + raise::KernelError(HERE, "Invalid dimension"); + } + if (fill_ratio * target_density > dens) { + return (target_density - dens) / target_density; + } else { + return ZERO; + } + } + }; + } // namespace arch #endif // ARCHETYPES_SPATIAL_DIST_HPP diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index c321d18f..48fdf869 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -351,6 +351,293 @@ namespace kernel { } }; // struct UniformInjector_kernel + template + struct NonUniformInjector_kernel { + static_assert(ED1::is_energy_dist, + "ED1 must be an energy distribution class"); + static_assert(ED2::is_energy_dist, + "ED2 must be an energy distribution class"); + static_assert(SD::is_spatial_dist, "SD must be a spatial distribution class"); + static_assert(M::is_metric, "M must be a metric class"); + + const real_t ppc0; + const spidx_t spidx1, spidx2; + + array_t i1s_1, i2s_1, i3s_1; + array_t dx1s_1, dx2s_1, dx3s_1; + array_t ux1s_1, ux2s_1, ux3s_1; + array_t phis_1; + array_t weights_1; + array_t tags_1; + + array_t i1s_2, i2s_2, i3s_2; + array_t dx1s_2, dx2s_2, dx3s_2; + array_t ux1s_2, ux2s_2, ux3s_2; + array_t phis_2; + array_t weights_2; + array_t tags_2; + + array_t idx { "idx" }; + + npart_t offset1, offset2; + M metric; + const ED1 energy_dist_1; + const ED2 energy_dist_2; + const SD spatial_dist; + const real_t inv_V0; + random_number_pool_t random_pool; + + NonUniformInjector_kernel(real_t ppc0, + spidx_t spidx1, + spidx_t spidx2, + Particles& species1, + Particles& species2, + npart_t offset1, + npart_t offset2, + const M& metric, + const ED1& energy_dist_1, + const ED2& energy_dist_2, + const SD& spatial_dist, + real_t inv_V0, + random_number_pool_t& random_pool) + : ppc0 { ppc0 } + , spidx1 { spidx1 } + , spidx2 { spidx2 } + , i1s_1 { species1.i1 } + , i2s_1 { species1.i2 } + , i3s_1 { species1.i3 } + , dx1s_1 { species1.dx1 } + , dx2s_1 { species1.dx2 } + , dx3s_1 { species1.dx3 } + , ux1s_1 { species1.ux1 } + , ux2s_1 { species1.ux2 } + , ux3s_1 { species1.ux3 } + , phis_1 { species1.phi } + , weights_1 { species1.weight } + , tags_1 { species1.tag } + , i1s_2 { species2.i1 } + , i2s_2 { species2.i2 } + , i3s_2 { species2.i3 } + , dx1s_2 { species2.dx1 } + , dx2s_2 { species2.dx2 } + , dx3s_2 { species2.dx3 } + , ux1s_2 { species2.ux1 } + , ux2s_2 { species2.ux2 } + , ux3s_2 { species2.ux3 } + , phis_2 { species2.phi } + , weights_2 { species2.weight } + , tags_2 { species2.tag } + , offset1 { offset1 } + , offset2 { offset2 } + , metric { metric } + , energy_dist_1 { energy_dist_1 } + , energy_dist_2 { energy_dist_2 } + , spatial_dist { spatial_dist } + , inv_V0 { inv_V0 } + , random_pool { random_pool } {} + + auto number_injected() const -> npart_t { + auto idx_h = Kokkos::create_mirror_view(idx); + Kokkos::deep_copy(idx_h, idx); + return idx_h(); + } + + Inline void operator()(index_t i1) const { + if constexpr (M::Dim == Dim::_1D) { + const auto i1_ = COORD(i1); + coord_t x_Cd { i1_ + HALF }; + coord_t x_Ph { ZERO }; + metric.template convert(x_Cd, x_Ph); + const auto ppc = static_cast(ppc0 * spatial_dist(x_Ph)); + if (ppc == 0) { + return; + } + auto rand_gen = random_pool.get_state(); + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + const auto dx1 = Random(rand_gen); + + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; + + vec_t v_T { ZERO }, v_XYZ { ZERO }; + energy_dist_1(x_Ph, v_T, spidx1); + metric.template transform_xyz(x_Cd, v_T, v_XYZ); + ux1s_1(index + offset1) = v_XYZ[0]; + ux2s_1(index + offset1) = v_XYZ[1]; + ux3s_1(index + offset1) = v_XYZ[2]; + energy_dist_2(x_Ph, v_T, spidx2); + metric.template transform_xyz(x_Cd, v_T, v_XYZ); + ux1s_2(index + offset2) = v_XYZ[0]; + ux2s_2(index + offset2) = v_XYZ[1]; + ux3s_2(index + offset2) = v_XYZ[2]; + + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; + if (M::CoordType == Coord::Cart) { + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; + } else { + const auto wei = metric.sqrt_det_h({ i1_ + HALF }) * inv_V0; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; + } + } + random_pool.free_state(rand_gen); + } else { + raise::KernelError(HERE, "NonUniformInjector_kernel 1D called for 2D/3D"); + } + } + + Inline void operator()(index_t i1, index_t i2) const { + if constexpr (M::Dim == Dim::_2D) { + const auto i1_ = COORD(i1); + const auto i2_ = COORD(i2); + coord_t x_Cd { i1_ + HALF, i2_ + HALF }; + coord_t x_Ph { ZERO }; + coord_t x_Cd_ { ZERO }; + x_Cd_[0] = x_Cd[0]; + x_Cd_[1] = x_Cd[1]; + if constexpr (S == SimEngine::SRPIC and M::CoordType != Coord::Cart) { + x_Cd_[2] = ZERO; + } + metric.template convert(x_Cd, x_Ph); + const auto ppc = static_cast(ppc0 * spatial_dist(x_Ph)); + if (ppc == 0) { + return; + } + auto rand_gen = random_pool.get_state(); + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); + + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; + + i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; + dx2s_1(index + offset1) = dx2; + i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; + dx2s_2(index + offset2) = dx2; + + vec_t v_T { ZERO }, v_Cd { ZERO }; + energy_dist_1(x_Ph, v_T, spidx1); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd_, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd_, v_T, v_Cd); + } + ux1s_1(index + offset1) = v_Cd[0]; + ux2s_1(index + offset1) = v_Cd[1]; + ux3s_1(index + offset1) = v_Cd[2]; + energy_dist_2(x_Ph, v_T, spidx2); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd_, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd_, v_T, v_Cd); + } + ux1s_2(index + offset2) = v_Cd[0]; + ux2s_2(index + offset2) = v_Cd[1]; + ux3s_2(index + offset2) = v_Cd[2]; + + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; + if (M::CoordType == Coord::Cart) { + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; + } else { + const auto wei = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF }) * + inv_V0; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; + } + } + random_pool.free_state(rand_gen); + } + + else { + raise::KernelError(HERE, "NonUniformInjector_kernel 2D called for 1D/3D"); + } + } + + Inline void operator()(index_t i1, index_t i2, index_t i3) const { + if constexpr (M::Dim == Dim::_3D) { + const auto i1_ = COORD(i1); + const auto i2_ = COORD(i2); + const auto i3_ = COORD(i3); + coord_t x_Cd { i1_ + HALF, i2_ + HALF, i3_ + HALF }; + coord_t x_Ph { ZERO }; + metric.template convert(x_Cd, x_Ph); + const auto ppc = static_cast(ppc0 * spatial_dist(x_Ph)); + if (ppc == 0) { + return; + } + auto rand_gen = random_pool.get_state(); + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); + const auto dx3 = Random(rand_gen); + + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; + + i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; + dx2s_1(index + offset1) = dx2; + i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; + dx2s_2(index + offset2) = dx2; + + i3s_1(index + offset1) = static_cast(i3) - N_GHOSTS; + dx3s_1(index + offset1) = dx3; + i3s_2(index + offset2) = static_cast(i3) - N_GHOSTS; + dx3s_2(index + offset2) = dx3; + + vec_t v_T { ZERO }, v_Cd { ZERO }; + energy_dist_1(x_Ph, v_T, spidx1); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd, v_T, v_Cd); + } + ux1s_1(index + offset1) = v_Cd[0]; + ux2s_1(index + offset1) = v_Cd[1]; + ux3s_1(index + offset1) = v_Cd[2]; + energy_dist_2(x_Ph, v_T, spidx2); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd, v_T, v_Cd); + } + ux1s_2(index + offset2) = v_Cd[0]; + ux2s_2(index + offset2) = v_Cd[1]; + ux3s_2(index + offset2) = v_Cd[2]; + + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; + if (M::CoordType == Coord::Cart) { + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; + } else { + const auto wei = metric.sqrt_det_h( + { i1_ + HALF, i2_ + HALF, i3_ + HALF }) * + inv_V0; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; + } + } + random_pool.free_state(rand_gen); + } else { + raise::KernelError(HERE, "NonUniformInjector_kernel 3D called for 1D/2D"); + } + } + }; // struct NonUniformInjector_kernel + } // namespace experimental template From 2d2fcff005c34f8949b11ad029face2c489e92ef Mon Sep 17 00:00:00 2001 From: haykh Date: Mon, 13 Oct 2025 23:18:35 -0400 Subject: [PATCH 3/8] shock injector (WIP) --- CMakeLists.txt | 2 +- pgens/shock/pgen.hpp | 173 ++++++--------- src/archetypes/particle_injector.h | 35 ++- src/kernels/injectors.hpp | 340 +++++++++++++++++++++++++++-- 4 files changed, 416 insertions(+), 134 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 29566913..b0329435 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ set(PROJECT_NAME entity) project( ${PROJECT_NAME} - VERSION 1.2.3 + VERSION 1.2.5 LANGUAGES CXX C) add_compile_options("-D ENTITY_VERSION=\"${PROJECT_VERSION}\"") set(hash_cmd "git diff --quiet src/ && echo $(git rev-parse HEAD) ") diff --git a/pgens/shock/pgen.hpp b/pgens/shock/pgen.hpp index 16f5197b..bee1fabd 100644 --- a/pgens/shock/pgen.hpp +++ b/pgens/shock/pgen.hpp @@ -91,7 +91,6 @@ namespace user { // injector properties const real_t injector_velocity; const simtime_t injection_start; - const real_t dt; const int injection_frequency; // magnetic field properties real_t Btheta, Bphi, Bmag; @@ -113,8 +112,8 @@ namespace user { , injector_velocity { p.template get("setup.injector_velocity", 1.0) } , injection_start { static_cast( p.template get("setup.injection_start", 0.0)) } - , injection_frequency { p.template get("setup.injection_frequency", 100) } - , dt { p.template get("algorithms.timestep.dt") } {} + , injection_frequency { p.template get("setup.injection_frequency", + 100) } {} inline PGen() {} @@ -159,39 +158,58 @@ namespace user { */ // minimum and maximum position of particles - real_t xg_min = global_xmin; - real_t xg_max = global_xmin + filling_fraction * (global_xmax - global_xmin); - - // 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); - } - } + // real_t xg_min = global_xmin; + // real_t xg_max = global_xmin + filling_fraction * (global_xmax - global_xmin); + // + // // 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); + // } + // } // species #1 -> e^- // species #2 -> protons // energy distribution of the particles - const auto temperatures = std::make_pair(temperature, - temperature_ratio * temperature); - const auto drifts = std::make_pair( - std::vector { -drift_ux, ZERO, ZERO }, - std::vector { -drift_ux, ZERO, ZERO }); - arch::InjectUniformMaxwellians(params, - domain, - ONE, - temperatures, - { 1, 2 }, - drifts, - false, - box); + // const auto temperatures = std::make_pair(temperature, + // temperature_ratio * temperature); + // const auto drifts = std::make_pair( + // std::vector { -drift_ux, ZERO, ZERO }, + // std::vector { -drift_ux, ZERO, ZERO }); + // arch::InjectUniformMaxwellians(params, + // domain, + // ONE, + // temperatures, + // { 1, 2 }, + // drifts, + // false, + // box); + // const auto maxwellian_1 = arch::experimental::Maxwellian( + // domain.mesh.metric, + // domain.random_pool, + // temperatures.first, + // { -drift_ux, ZERO, ZERO }); + // const auto maxwellian_2 = arch::experimental::Maxwellian( + // domain.mesh.metric, + // domain.random_pool, + // temperatures.second, + // { -drift_ux, ZERO, ZERO }); + // arch::InjectReplenishConst( + // params, + // domain, + // ONE, + // { 1, 2 }, + // maxwellian_1, + // maxwellian_2, + // 0.95, + // box); } void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { @@ -212,91 +230,31 @@ namespace user { */ // check if the injector should be active - if (step % injection_frequency != 0 and time < injection_start) { + if ((step % injection_frequency) != 0 or + ((time < injection_start) and (step != 0))) { return; } + const auto dt = params.template get("algorithms.timestep.dt"); + // initial position of injector const auto x_init = global_xmin + filling_fraction * (global_xmax - global_xmin); // compute the position of the injector after the current timestep - auto xmax = x_init + injector_velocity * - (std::max(time - injection_start, ZERO) + dt); + auto xmax = x_init + injector_velocity * (time + dt - injection_start); if (xmax >= global_xmax) { xmax = global_xmax; } // compute the beginning of the injected region - auto xmin = xmax - 1.1 * injection_frequency * dt; + auto xmin = xmax - 2.2 * static_cast(injection_frequency) * dt; if (xmin <= global_xmin) { xmin = global_xmin; } - - // // define indice range to reset fields - // boundaries_t incl_ghosts; - // for (auto d = 0; d < M::Dim; ++d) { - // incl_ghosts.push_back({ false, false }); - // } - - // // define box to reset fields - // boundaries_t purge_box; - // // loop over all dimension - // for (auto d = 0u; d < M::Dim; ++d) { - // if (d == 0) { - // purge_box.push_back({ xmin, global_xmax }); - // } else { - // purge_box.push_back(Range::All); - // } - // } - // - // const auto extent = domain.mesh.ExtentToRange(purge_box, incl_ghosts); - // tuple_t x_min { 0 }, x_max { 0 }; - // for (auto d = 0; d < M::Dim; ++d) { - // x_min[d] = extent[d].first; - // x_max[d] = extent[d].second; - // } - // - // Kokkos::parallel_for("ResetFields", - // CreateRangePolicy(x_min, x_max), - // arch::SetEMFields_kernel { - // domain.fields.em, - // init_flds, - // domain.mesh.metric }); - // - // global_domain.CommunicateFields(domain, Comm::B | Comm::E); - - // /* - // tag particles inside the injection zone as dead - // */ - // const auto& mesh = domain.mesh; - // - // // loop over particle species - // for (auto s { 0u }; s < 2; ++s) { - // // get particle properties - // auto& species = domain.species[s]; - // auto i1 = species.i1; - // auto dx1 = species.dx1; - // auto tag = species.tag; - // - // Kokkos::parallel_for( - // "RemoveParticles", - // species.rangeActiveParticles(), - // Lambda(index_t p) { - // // check if the particle is already dead - // if (tag(p) == ParticleTag::dead) { - // return; - // } - // const auto x_Cd = static_cast(i1(p)) + - // static_cast(dx1(p)); - // const auto x_Ph = mesh.metric.template convert<1, Crd::Cd, Crd::XYZ>( - // x_Cd); - // - // if (x_Ph > xmin) { - // tag(p) = ParticleTag::dead; - // } - // }); - // } + if (step == 0) { + xmin = global_xmin; + } /* Inject slab of fresh plasma @@ -330,7 +288,9 @@ namespace user { temperature_2, { -drift_ux, ZERO, ZERO }); - const auto& mesh = domain.mesh; + const auto& mesh = domain.mesh; + const auto inj_vel = injector_velocity; + for (auto& s : { 1u, 2u }) { auto& species = domain.species[s - 1]; auto i1 = species.i1; @@ -365,13 +325,9 @@ namespace user { ux2(p) = v_Cd[1]; ux3(p) = v_Cd[2]; } else if (x_Ph > xmax) { - if (injector_velocity < 0.9) { - ux1(p) = TWO * injector_velocity / - math::sqrt(ONE - SQR(injector_velocity)); - } else { - ux1(p) = TWO * injector_velocity / - math::sqrt(ONE - SQR(injector_velocity) + 0.001); - } + ux1(p) = TWO * inj_vel / + math::max(static_cast(1e-2), + math::sqrt(ONE - SQR(inj_vel))); ux2(p) = ZERO; ux3(p) = ZERO; } @@ -385,6 +341,7 @@ namespace user { { 1, 2 }, maxwellian_1, maxwellian_2, + 1.0, inj_box); } }; diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h index 4ab3db6a..514e992e 100644 --- a/src/archetypes/particle_injector.h +++ b/src/archetypes/particle_injector.h @@ -745,7 +745,8 @@ namespace arch { const std::pair& species, const ED1& energy_dist_1, const ED2& energy_dist_2, - const boundaries_t& box = {}) { + real_t density_tolerance = 0.9, + const boundaries_t& box = {}) { static_assert(M::is_metric, "M must be a metric class"); static_assert(ED1::is_energy_dist, "ED1 must be an energy distribution class"); static_assert(ED2::is_energy_dist, "ED2 must be an energy distribution class"); @@ -787,7 +788,7 @@ namespace arch { Kokkos::parallel_for( "ComputeMoments", prtl_spec.rangeActiveParticles(), - kernel::ParticleMoments_kernel({}, scatter_buff, buff_idx, + kernel::ParticleMoments_kernel({}, scatter_buff, buff_idx, prtl_spec.i1, prtl_spec.i2, prtl_spec.i3, prtl_spec.dx1, prtl_spec.dx2, prtl_spec.dx3, prtl_spec.ux1, prtl_spec.ux2, prtl_spec.ux3, @@ -800,17 +801,17 @@ namespace arch { } Kokkos::Experimental::contribute(domain.fields.buff, scatter_buff); } - const auto spatial_dist = ReplenishConst { domain.mesh.metric, - domain.fields.buff, - buff_idx, - target_tot_ndens, - 1.0 }; + // const auto spatial_dist = ReplenishConst { domain.mesh.metric, + // domain.fields.buff, + // buff_idx, + // target_tot_ndens, + // density_tolerance }; const auto ppc = target_tot_ndens * params.template get("particles.ppc0") * HALF; auto injector_kernel = - kernel::experimental::NonUniformInjector_kernel( + kernel::experimental::ConstNPPCInjector_kernel( ppc, species.first, species.second, @@ -821,9 +822,25 @@ namespace arch { domain.mesh.metric, energy_dist_1, energy_dist_2, - spatial_dist, + domain.fields.buff, + buff_idx, ONE / params.template get("scales.V0"), domain.random_pool); + // auto injector_kernel = + // kernel::experimental::NonUniformInjector_kernel( + // ppc, + // species.first, + // species.second, + // domain.species[species.first - 1], + // domain.species[species.second - 1], + // domain.species[species.first - 1].npart(), + // domain.species[species.second - 1].npart(), + // domain.mesh.metric, + // energy_dist_1, + // energy_dist_2, + // spatial_dist, + // ONE / params.template get("scales.V0"), + // domain.random_pool); Kokkos::parallel_for("InjectReplenishConst", cell_range, injector_kernel); diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index 48fdf869..14d94a72 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -448,11 +448,12 @@ namespace kernel { coord_t x_Cd { i1_ + HALF }; coord_t x_Ph { ZERO }; metric.template convert(x_Cd, x_Ph); - const auto ppc = static_cast(ppc0 * spatial_dist(x_Ph)); - if (ppc == 0) { - return; + const auto ppc_real = static_cast(ppc0) * spatial_dist(x_Ph); + auto ppc = static_cast(ppc_real); + auto rand_gen = random_pool.get_state(); + if (Random(rand_gen) < (ppc_real - math::floor(ppc_real))) { + ++ppc; } - auto rand_gen = random_pool.get_state(); for (auto p { 0u }; p < ppc; ++p) { const auto index = Kokkos::atomic_fetch_add(&idx(), 1); const auto dx1 = Random(rand_gen); @@ -504,15 +505,19 @@ namespace kernel { x_Cd_[2] = ZERO; } metric.template convert(x_Cd, x_Ph); - const auto ppc = static_cast(ppc0 * spatial_dist(x_Ph)); - if (ppc == 0) { - return; + const auto ppc_real = static_cast(ppc0) * spatial_dist(x_Ph); + auto ppc = static_cast(ppc_real); + auto rand_gen = random_pool.get_state(); + if (Random(rand_gen) < (ppc_real - math::floor(ppc_real))) { + ++ppc; } - auto rand_gen = random_pool.get_state(); + random_pool.free_state(rand_gen); for (auto p { 0u }; p < ppc; ++p) { - const auto index = Kokkos::atomic_fetch_add(&idx(), 1); - const auto dx1 = Random(rand_gen); - const auto dx2 = Random(rand_gen); + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + auto rand_gen = random_pool.get_state(); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); + random_pool.free_state(rand_gen); i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; dx1s_1(index + offset1) = dx1; @@ -556,7 +561,6 @@ namespace kernel { weights_2(index + offset2) = wei; } } - random_pool.free_state(rand_gen); } else { @@ -572,11 +576,12 @@ namespace kernel { coord_t x_Cd { i1_ + HALF, i2_ + HALF, i3_ + HALF }; coord_t x_Ph { ZERO }; metric.template convert(x_Cd, x_Ph); - const auto ppc = static_cast(ppc0 * spatial_dist(x_Ph)); - if (ppc == 0) { - return; + const auto ppc_real = static_cast(ppc0) * spatial_dist(x_Ph); + auto ppc = static_cast(ppc_real); + auto rand_gen = random_pool.get_state(); + if (Random(rand_gen) < (ppc_real - math::floor(ppc_real))) { + ++ppc; } - auto rand_gen = random_pool.get_state(); for (auto p { 0u }; p < ppc; ++p) { const auto index = Kokkos::atomic_fetch_add(&idx(), 1); const auto dx1 = Random(rand_gen); @@ -638,6 +643,309 @@ namespace kernel { } }; // struct NonUniformInjector_kernel + template + struct ConstNPPCInjector_kernel { + static_assert(ED1::is_energy_dist, + "ED1 must be an energy distribution class"); + static_assert(ED2::is_energy_dist, + "ED2 must be an energy distribution class"); + static_assert(M::is_metric, "M must be a metric class"); + + const real_t ppc0; + const spidx_t spidx1, spidx2; + + array_t i1s_1, i2s_1, i3s_1; + array_t dx1s_1, dx2s_1, dx3s_1; + array_t ux1s_1, ux2s_1, ux3s_1; + array_t phis_1; + array_t weights_1; + array_t tags_1; + + array_t i1s_2, i2s_2, i3s_2; + array_t dx1s_2, dx2s_2, dx3s_2; + array_t ux1s_2, ux2s_2, ux3s_2; + array_t phis_2; + array_t weights_2; + array_t tags_2; + + array_t idx { "idx" }; + + npart_t offset1, offset2; + M metric; + const ED1 energy_dist_1; + const ED2 energy_dist_2; + const ndfield_t nppc_array; + const idx_t nppc_idx; + const real_t inv_V0; + random_number_pool_t random_pool; + + ConstNPPCInjector_kernel(real_t ppc0, + spidx_t spidx1, + spidx_t spidx2, + Particles& species1, + Particles& species2, + npart_t offset1, + npart_t offset2, + const M& metric, + const ED1& energy_dist_1, + const ED2& energy_dist_2, + const ndfield_t& nppc_array, + idx_t nppc_idx, + real_t inv_V0, + random_number_pool_t& random_pool) + : ppc0 { ppc0 } + , spidx1 { spidx1 } + , spidx2 { spidx2 } + , i1s_1 { species1.i1 } + , i2s_1 { species1.i2 } + , i3s_1 { species1.i3 } + , dx1s_1 { species1.dx1 } + , dx2s_1 { species1.dx2 } + , dx3s_1 { species1.dx3 } + , ux1s_1 { species1.ux1 } + , ux2s_1 { species1.ux2 } + , ux3s_1 { species1.ux3 } + , phis_1 { species1.phi } + , weights_1 { species1.weight } + , tags_1 { species1.tag } + , i1s_2 { species2.i1 } + , i2s_2 { species2.i2 } + , i3s_2 { species2.i3 } + , dx1s_2 { species2.dx1 } + , dx2s_2 { species2.dx2 } + , dx3s_2 { species2.dx3 } + , ux1s_2 { species2.ux1 } + , ux2s_2 { species2.ux2 } + , ux3s_2 { species2.ux3 } + , phis_2 { species2.phi } + , weights_2 { species2.weight } + , tags_2 { species2.tag } + , offset1 { offset1 } + , offset2 { offset2 } + , metric { metric } + , energy_dist_1 { energy_dist_1 } + , energy_dist_2 { energy_dist_2 } // , spatial_dist { spatial_dist } + , nppc_array { nppc_array } + , nppc_idx { nppc_idx } + , inv_V0 { inv_V0 } + , random_pool { random_pool } {} + + auto number_injected() const -> npart_t { + auto idx_h = Kokkos::create_mirror_view(idx); + Kokkos::deep_copy(idx_h, idx); + return idx_h(); + } + + Inline auto get_inj_ppc(real_t& ppc_real) const -> npart_t { + // const auto delta_ppc = ; + // if (delta_ppc > 100.0) { + // Kokkos::printf( + // "Warning: ConstNPPCInjector_kernel::get_inj_ppc(): " + // "ppc_real too small compared to ppc0: ppc0=%f, ppc_real=%f\n", + // ppc0, + // ppc_real); + // Kokkos::abort("Error: ConstNPPCInjector_kernel::get_inj_ppc(): " + // "ppc_real too small compared to ppc0"); + // } + return static_cast(math::max(ppc0 - ppc_real, ZERO)); + } + + Inline void operator()(index_t i1) const { + if constexpr (M::Dim == Dim::_1D) { + // coord_t x_Cd { i1_ + HALF }; + // coord_t x_Ph { ZERO }; + // metric.template convert(x_Cd, x_Ph); + const auto ppc = get_inj_ppc(nppc_array(i1, nppc_idx)); + const auto i1_ = COORD(i1); + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + auto rand_gen = random_pool.get_state(); + const auto dx1 = Random(rand_gen); + random_pool.free_state(rand_gen); + + coord_t x_Cd { i1_ + dx1 }; + coord_t x_Ph { ZERO }; + metric.template convert(x_Cd, x_Ph); + + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; + + vec_t v_T { ZERO }, v_XYZ { ZERO }; + energy_dist_1(x_Ph, v_T, spidx1); + metric.template transform_xyz(x_Cd, v_T, v_XYZ); + ux1s_1(index + offset1) = v_XYZ[0]; + ux2s_1(index + offset1) = v_XYZ[1]; + ux3s_1(index + offset1) = v_XYZ[2]; + energy_dist_2(x_Ph, v_T, spidx2); + metric.template transform_xyz(x_Cd, v_T, v_XYZ); + ux1s_2(index + offset2) = v_XYZ[0]; + ux2s_2(index + offset2) = v_XYZ[1]; + ux3s_2(index + offset2) = v_XYZ[2]; + + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; + if (M::CoordType == Coord::Cart) { + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; + } else { + const auto wei = metric.sqrt_det_h({ i1_ + HALF }) * inv_V0; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; + } + } + } else { + raise::KernelError(HERE, "NonUniformInjector_kernel 1D called for 2D/3D"); + } + } + + Inline void operator()(index_t i1, index_t i2) const { + if constexpr (M::Dim == Dim::_2D) { + const auto ppc = get_inj_ppc(nppc_array(i1, i2, nppc_idx)); + const auto i1_ = COORD(i1); + const auto i2_ = COORD(i2); + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + auto rand_gen = random_pool.get_state(); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); + random_pool.free_state(rand_gen); + + coord_t x_Cd { i1_ + dx1, i2_ + dx2 }; + coord_t x_Ph { ZERO }; + coord_t x_Cd_ { ZERO }; + x_Cd_[0] = x_Cd[0]; + x_Cd_[1] = x_Cd[1]; + if constexpr (S == SimEngine::SRPIC and M::CoordType != Coord::Cart) { + x_Cd_[2] = ZERO; + } + metric.template convert(x_Cd, x_Ph); + + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; + + i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; + dx2s_1(index + offset1) = dx2; + i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; + dx2s_2(index + offset2) = dx2; + + vec_t v_T { ZERO }, v_Cd { ZERO }; + energy_dist_1(x_Ph, v_T, spidx1); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd_, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd_, v_T, v_Cd); + } + ux1s_1(index + offset1) = v_Cd[0]; + ux2s_1(index + offset1) = v_Cd[1]; + ux3s_1(index + offset1) = v_Cd[2]; + energy_dist_2(x_Ph, v_T, spidx2); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd_, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd_, v_T, v_Cd); + } + ux1s_2(index + offset2) = v_Cd[0]; + ux2s_2(index + offset2) = v_Cd[1]; + ux3s_2(index + offset2) = v_Cd[2]; + + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; + if (M::CoordType == Coord::Cart) { + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; + } else { + const auto wei = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF }) * + inv_V0; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; + } + } + } + + else { + raise::KernelError(HERE, "NonUniformInjector_kernel 2D called for 1D/3D"); + } + } + + Inline void operator()(index_t i1, index_t i2, index_t i3) const { + // if constexpr (M::Dim == Dim::_3D) { + // const auto i1_ = COORD(i1); + // const auto i2_ = COORD(i2); + // const auto i3_ = COORD(i3); + // coord_t x_Cd { i1_ + HALF, i2_ + HALF, i3_ + HALF }; + // coord_t x_Ph { ZERO }; + // metric.template convert(x_Cd, x_Ph); + // const auto ppc_real = static_cast(ppc0) * spatial_dist(x_Ph); + // auto ppc = static_cast(ppc_real); + // auto rand_gen = random_pool.get_state(); + // if (Random(rand_gen) < (ppc_real - math::floor(ppc_real))) { + // ++ppc; + // } + // for (auto p { 0u }; p < ppc; ++p) { + // const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + // const auto dx1 = Random(rand_gen); + // const auto dx2 = Random(rand_gen); + // const auto dx3 = Random(rand_gen); + // + // i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + // dx1s_1(index + offset1) = dx1; + // i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + // dx1s_2(index + offset2) = dx1; + // + // i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; + // dx2s_1(index + offset1) = dx2; + // i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; + // dx2s_2(index + offset2) = dx2; + // + // i3s_1(index + offset1) = static_cast(i3) - N_GHOSTS; + // dx3s_1(index + offset1) = dx3; + // i3s_2(index + offset2) = static_cast(i3) - N_GHOSTS; + // dx3s_2(index + offset2) = dx3; + // + // vec_t v_T { ZERO }, v_Cd { ZERO }; + // energy_dist_1(x_Ph, v_T, spidx1); + // if constexpr (S == SimEngine::SRPIC) { + // metric.template transform_xyz(x_Cd, v_T, v_Cd); + // } else if constexpr (S == SimEngine::GRPIC) { + // metric.template transform(x_Cd, v_T, v_Cd); + // } + // ux1s_1(index + offset1) = v_Cd[0]; + // ux2s_1(index + offset1) = v_Cd[1]; + // ux3s_1(index + offset1) = v_Cd[2]; + // energy_dist_2(x_Ph, v_T, spidx2); + // if constexpr (S == SimEngine::SRPIC) { + // metric.template transform_xyz(x_Cd, v_T, v_Cd); + // } else if constexpr (S == SimEngine::GRPIC) { + // metric.template transform(x_Cd, v_T, v_Cd); + // } + // ux1s_2(index + offset2) = v_Cd[0]; + // ux2s_2(index + offset2) = v_Cd[1]; + // ux3s_2(index + offset2) = v_Cd[2]; + // + // tags_1(index + offset1) = ParticleTag::alive; + // tags_2(index + offset2) = ParticleTag::alive; + // if (M::CoordType == Coord::Cart) { + // weights_1(index + offset1) = ONE; + // weights_2(index + offset2) = ONE; + // } else { + // const auto wei = metric.sqrt_det_h( + // { i1_ + HALF, i2_ + HALF, i3_ + HALF }) * + // inv_V0; + // weights_1(index + offset1) = wei; + // weights_2(index + offset2) = wei; + // } + // } + // random_pool.free_state(rand_gen); + // } else { + // raise::KernelError(HERE, "NonUniformInjector_kernel 3D called for 1D/2D"); + // } + } + }; // struct NonUniformInjector_kernel + } // namespace experimental template From 3a13be199baa129b8d227cc0ea52d1798ea65fdb Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 14 Oct 2025 10:36:46 -0500 Subject: [PATCH 4/8] bugfix --- src/kernels/injectors.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index 14d94a72..ee103c27 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -747,7 +747,7 @@ namespace kernel { // Kokkos::abort("Error: ConstNPPCInjector_kernel::get_inj_ppc(): " // "ppc_real too small compared to ppc0"); // } - return static_cast(math::max(ppc0 - ppc_real, ZERO)); + return static_cast(math::max(ppc0 - ppc_real, static_cast(0.0))); } Inline void operator()(index_t i1) const { From e5b663103163a14ff317b742edb0f3e0de9ea584 Mon Sep 17 00:00:00 2001 From: haykh Date: Sat, 8 Nov 2025 10:44:52 -0800 Subject: [PATCH 5/8] finalized shock pgen --- pgens/shock/pgen.hpp | 74 +------------------------------------------- 1 file changed, 1 insertion(+), 73 deletions(-) diff --git a/pgens/shock/pgen.hpp b/pgens/shock/pgen.hpp index bee1fabd..8bafc273 100644 --- a/pgens/shock/pgen.hpp +++ b/pgens/shock/pgen.hpp @@ -9,13 +9,10 @@ #include "utils/numeric.h" #include "archetypes/energy_dist.h" -#include "archetypes/field_setter.h" #include "archetypes/particle_injector.h" #include "archetypes/problem_generator.h" -#include "archetypes/utils.h" #include "framework/domain/metadomain.h" -#include #include namespace user { @@ -141,76 +138,7 @@ namespace user { } } - inline void InitPrtls(Domain& domain) { - - /* - * Plasma setup as partially filled box - * - * Plasma setup: - * - * global_xmin global_xmax - * | | - * V V - * |:::::::::::|..........................| - * ^ - * | - * filling_fraction - */ - - // minimum and maximum position of particles - // real_t xg_min = global_xmin; - // real_t xg_max = global_xmin + filling_fraction * (global_xmax - global_xmin); - // - // // 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); - // } - // } - - // species #1 -> e^- - // species #2 -> protons - - // energy distribution of the particles - // const auto temperatures = std::make_pair(temperature, - // temperature_ratio * temperature); - // const auto drifts = std::make_pair( - // std::vector { -drift_ux, ZERO, ZERO }, - // std::vector { -drift_ux, ZERO, ZERO }); - // arch::InjectUniformMaxwellians(params, - // domain, - // ONE, - // temperatures, - // { 1, 2 }, - // drifts, - // false, - // box); - // const auto maxwellian_1 = arch::experimental::Maxwellian( - // domain.mesh.metric, - // domain.random_pool, - // temperatures.first, - // { -drift_ux, ZERO, ZERO }); - // const auto maxwellian_2 = arch::experimental::Maxwellian( - // domain.mesh.metric, - // domain.random_pool, - // temperatures.second, - // { -drift_ux, ZERO, ZERO }); - // arch::InjectReplenishConst( - // params, - // domain, - // ONE, - // { 1, 2 }, - // maxwellian_1, - // maxwellian_2, - // 0.95, - // box); - } + inline void InitPrtls(Domain&) {} void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { From 9f3bb058bdda17c386bc45f34488e92d5873a9f3 Mon Sep 17 00:00:00 2001 From: haykh Date: Sat, 8 Nov 2025 10:54:51 -0800 Subject: [PATCH 6/8] nonuniform inj kernel 3d --- src/kernels/injectors.hpp | 159 ++++++++++++++++++-------------------- 1 file changed, 74 insertions(+), 85 deletions(-) diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index ee103c27..1f860d04 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -737,24 +737,12 @@ namespace kernel { } Inline auto get_inj_ppc(real_t& ppc_real) const -> npart_t { - // const auto delta_ppc = ; - // if (delta_ppc > 100.0) { - // Kokkos::printf( - // "Warning: ConstNPPCInjector_kernel::get_inj_ppc(): " - // "ppc_real too small compared to ppc0: ppc0=%f, ppc_real=%f\n", - // ppc0, - // ppc_real); - // Kokkos::abort("Error: ConstNPPCInjector_kernel::get_inj_ppc(): " - // "ppc_real too small compared to ppc0"); - // } - return static_cast(math::max(ppc0 - ppc_real, static_cast(0.0))); + return static_cast( + math::max(ppc0 - ppc_real, static_cast(0.0))); } Inline void operator()(index_t i1) const { if constexpr (M::Dim == Dim::_1D) { - // coord_t x_Cd { i1_ + HALF }; - // coord_t x_Ph { ZERO }; - // metric.template convert(x_Cd, x_Ph); const auto ppc = get_inj_ppc(nppc_array(i1, nppc_idx)); const auto i1_ = COORD(i1); for (auto p { 0u }; p < ppc; ++p) { @@ -872,77 +860,78 @@ namespace kernel { } Inline void operator()(index_t i1, index_t i2, index_t i3) const { - // if constexpr (M::Dim == Dim::_3D) { - // const auto i1_ = COORD(i1); - // const auto i2_ = COORD(i2); - // const auto i3_ = COORD(i3); - // coord_t x_Cd { i1_ + HALF, i2_ + HALF, i3_ + HALF }; - // coord_t x_Ph { ZERO }; - // metric.template convert(x_Cd, x_Ph); - // const auto ppc_real = static_cast(ppc0) * spatial_dist(x_Ph); - // auto ppc = static_cast(ppc_real); - // auto rand_gen = random_pool.get_state(); - // if (Random(rand_gen) < (ppc_real - math::floor(ppc_real))) { - // ++ppc; - // } - // for (auto p { 0u }; p < ppc; ++p) { - // const auto index = Kokkos::atomic_fetch_add(&idx(), 1); - // const auto dx1 = Random(rand_gen); - // const auto dx2 = Random(rand_gen); - // const auto dx3 = Random(rand_gen); - // - // i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; - // dx1s_1(index + offset1) = dx1; - // i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; - // dx1s_2(index + offset2) = dx1; - // - // i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; - // dx2s_1(index + offset1) = dx2; - // i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; - // dx2s_2(index + offset2) = dx2; - // - // i3s_1(index + offset1) = static_cast(i3) - N_GHOSTS; - // dx3s_1(index + offset1) = dx3; - // i3s_2(index + offset2) = static_cast(i3) - N_GHOSTS; - // dx3s_2(index + offset2) = dx3; - // - // vec_t v_T { ZERO }, v_Cd { ZERO }; - // energy_dist_1(x_Ph, v_T, spidx1); - // if constexpr (S == SimEngine::SRPIC) { - // metric.template transform_xyz(x_Cd, v_T, v_Cd); - // } else if constexpr (S == SimEngine::GRPIC) { - // metric.template transform(x_Cd, v_T, v_Cd); - // } - // ux1s_1(index + offset1) = v_Cd[0]; - // ux2s_1(index + offset1) = v_Cd[1]; - // ux3s_1(index + offset1) = v_Cd[2]; - // energy_dist_2(x_Ph, v_T, spidx2); - // if constexpr (S == SimEngine::SRPIC) { - // metric.template transform_xyz(x_Cd, v_T, v_Cd); - // } else if constexpr (S == SimEngine::GRPIC) { - // metric.template transform(x_Cd, v_T, v_Cd); - // } - // ux1s_2(index + offset2) = v_Cd[0]; - // ux2s_2(index + offset2) = v_Cd[1]; - // ux3s_2(index + offset2) = v_Cd[2]; - // - // tags_1(index + offset1) = ParticleTag::alive; - // tags_2(index + offset2) = ParticleTag::alive; - // if (M::CoordType == Coord::Cart) { - // weights_1(index + offset1) = ONE; - // weights_2(index + offset2) = ONE; - // } else { - // const auto wei = metric.sqrt_det_h( - // { i1_ + HALF, i2_ + HALF, i3_ + HALF }) * - // inv_V0; - // weights_1(index + offset1) = wei; - // weights_2(index + offset2) = wei; - // } - // } - // random_pool.free_state(rand_gen); - // } else { - // raise::KernelError(HERE, "NonUniformInjector_kernel 3D called for 1D/2D"); - // } + if constexpr (M::Dim == Dim::_3D) { + const auto ppc = get_inj_ppc(nppc_array(i1, i2, i3, nppc_idx)); + const auto i1_ = COORD(i1); + const auto i2_ = COORD(i2); + const auto i3_ = COORD(i3); + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + auto rand_gen = random_pool.get_state(); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); + const auto dx3 = Random(rand_gen); + random_pool.free_state(rand_gen); + + coord_t x_Cd { i1_ + dx1, i2_ + dx2, i3_ + dx3 }; + coord_t x_Ph { ZERO }; + coord_t x_Cd_ { ZERO }; + x_Cd_[0] = x_Cd[0]; + x_Cd_[1] = x_Cd[1]; + x_Cd_[2] = x_Cd[2]; + metric.template convert(x_Cd, x_Ph); + + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; + + i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; + dx2s_1(index + offset1) = dx2; + i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; + dx2s_2(index + offset2) = dx2; + + i3s_1(index + offset1) = static_cast(i3) - N_GHOSTS; + dx3s_1(index + offset1) = dx3; + i3s_2(index + offset2) = static_cast(i3) - N_GHOSTS; + dx3s_2(index + offset2) = dx3; + + vec_t v_T { ZERO }, v_Cd { ZERO }; + energy_dist_1(x_Ph, v_T, spidx1); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd, v_T, v_Cd); + } + ux1s_1(index + offset1) = v_Cd[0]; + ux2s_1(index + offset1) = v_Cd[1]; + ux3s_1(index + offset1) = v_Cd[2]; + energy_dist_2(x_Ph, v_T, spidx2); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd, v_T, v_Cd); + } + ux1s_2(index + offset2) = v_Cd[0]; + ux2s_2(index + offset2) = v_Cd[1]; + ux3s_2(index + offset2) = v_Cd[2]; + + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; + if (M::CoordType == Coord::Cart) { + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; + } else { + const auto wei = metric.sqrt_det_h( + { i1_ + HALF, i2_ + HALF, i3_ + HALF }) * + inv_V0; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; + } + } + } else { + raise::KernelError(HERE, "NonUniformInjector_kernel 3D called for 1D/2D"); + } } }; // struct NonUniformInjector_kernel From 56c8249f564a6541c9bd94e931e22ee7eb5ab420 Mon Sep 17 00:00:00 2001 From: haykh Date: Sat, 8 Nov 2025 10:55:04 -0800 Subject: [PATCH 7/8] version bump --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b0329435..4cff5b41 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ set(PROJECT_NAME entity) project( ${PROJECT_NAME} - VERSION 1.2.5 + VERSION 1.3.0 LANGUAGES CXX C) add_compile_options("-D ENTITY_VERSION=\"${PROJECT_VERSION}\"") set(hash_cmd "git diff --quiet src/ && echo $(git rev-parse HEAD) ") From 5054ec813ff054ee74feba1f301d8c503f000d0f Mon Sep 17 00:00:00 2001 From: haykh Date: Sat, 8 Nov 2025 10:59:00 -0800 Subject: [PATCH 8/8] experiment rm --- src/kernels/injectors.hpp | 1613 +++++++++++++------------------------ 1 file changed, 581 insertions(+), 1032 deletions(-) diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index 1f860d04..0d30839c 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -5,6 +5,7 @@ * - kernel::UniformInjector_kernel<> * - kernel::GlobalInjector_kernel<> * - kernel::NonUniformInjector_kernel<> + * - kernel::ConstNPPCInjector_kernel<> * @namespaces: * - kernel:: */ @@ -25,9 +26,10 @@ namespace kernel { using namespace ntt; - template + template struct UniformInjector_kernel { - static_assert(ED::is_energy_dist, "ED must be an energy distribution class"); + static_assert(ED1::is_energy_dist, "ED1 must be an energy distribution class"); + static_assert(ED2::is_energy_dist, "ED2 must be an energy distribution class"); static_assert(M::is_metric, "M must be a metric class"); const spidx_t spidx1, spidx2; @@ -49,7 +51,8 @@ namespace kernel { npart_t offset1, offset2; const M metric; const array_t xi_min, xi_max; - const ED energy_dist; + const ED1 energy_dist_1; + const ED2 energy_dist_2; const real_t inv_V0; random_number_pool_t random_pool; @@ -62,7 +65,8 @@ namespace kernel { const M& metric, const array_t& xi_min, const array_t& xi_max, - const ED& energy_dist, + const ED1& energy_dist_1, + const ED2& energy_dist_2, real_t inv_V0, random_number_pool_t& random_pool) : spidx1 { spidx1 } @@ -96,7 +100,8 @@ namespace kernel { , metric { metric } , xi_min { xi_min } , xi_max { xi_max } - , energy_dist { energy_dist } + , energy_dist_1 { energy_dist_1 } + , energy_dist_2 { energy_dist_2 } , inv_V0 { inv_V0 } , random_pool { random_pool } {} @@ -118,23 +123,23 @@ namespace kernel { coord_t x_Ph { ZERO }; metric.template convert(x_Cd, x_Ph); if constexpr (M::CoordType == Coord::Cart) { - energy_dist(x_Ph, v1, spidx1); - energy_dist(x_Ph, v2, spidx2); + energy_dist_1(x_Ph, v1, spidx1); + energy_dist_2(x_Ph, v2, spidx2); } else if constexpr (S == SimEngine::SRPIC) { coord_t x_Cd_ { ZERO }; x_Cd_[0] = x_Cd[0]; x_Cd_[1] = x_Cd[1]; x_Cd_[2] = ZERO; // phi = 0 vec_t v_Ph { ZERO }; - energy_dist(x_Ph, v_Ph, spidx1); + energy_dist_1(x_Ph, v_Ph, spidx1); metric.template transform_xyz(x_Cd_, v_Ph, v1); - energy_dist(x_Ph, v_Ph, spidx2); + energy_dist_2(x_Ph, v_Ph, spidx2); metric.template transform_xyz(x_Cd_, v_Ph, v2); } else if constexpr (S == SimEngine::GRPIC) { vec_t v_Ph { ZERO }; - energy_dist(x_Ph, v_Ph, spidx1); + energy_dist_1(x_Ph, v_Ph, spidx1); metric.template transform(x_Cd, v_Ph, v1); - energy_dist(x_Ph, v_Ph, spidx2); + energy_dist_2(x_Ph, v_Ph, spidx2); metric.template transform(x_Cd, v_Ph, v2); } else { raise::KernelError(HERE, "Unknown simulation engine"); @@ -183,994 +188,10 @@ namespace kernel { } }; // struct UniformInjector_kernel - namespace experimental { - - template - struct UniformInjector_kernel { - static_assert(ED1::is_energy_dist, - "ED1 must be an energy distribution class"); - static_assert(ED2::is_energy_dist, - "ED2 must be an energy distribution class"); - static_assert(M::is_metric, "M must be a metric class"); - - const spidx_t spidx1, spidx2; - - array_t i1s_1, i2s_1, i3s_1; - array_t dx1s_1, dx2s_1, dx3s_1; - array_t ux1s_1, ux2s_1, ux3s_1; - array_t phis_1; - array_t weights_1; - array_t tags_1; - - array_t i1s_2, i2s_2, i3s_2; - array_t dx1s_2, dx2s_2, dx3s_2; - array_t ux1s_2, ux2s_2, ux3s_2; - array_t phis_2; - array_t weights_2; - array_t tags_2; - - npart_t offset1, offset2; - const M metric; - const array_t xi_min, xi_max; - const ED1 energy_dist_1; - const ED2 energy_dist_2; - const real_t inv_V0; - random_number_pool_t random_pool; - - UniformInjector_kernel(spidx_t spidx1, - spidx_t spidx2, - Particles& species1, - Particles& species2, - npart_t offset1, - npart_t offset2, - const M& metric, - const array_t& xi_min, - const array_t& xi_max, - const ED1& energy_dist_1, - const ED2& energy_dist_2, - real_t inv_V0, - random_number_pool_t& random_pool) - : spidx1 { spidx1 } - , spidx2 { spidx2 } - , i1s_1 { species1.i1 } - , i2s_1 { species1.i2 } - , i3s_1 { species1.i3 } - , dx1s_1 { species1.dx1 } - , dx2s_1 { species1.dx2 } - , dx3s_1 { species1.dx3 } - , ux1s_1 { species1.ux1 } - , ux2s_1 { species1.ux2 } - , ux3s_1 { species1.ux3 } - , phis_1 { species1.phi } - , weights_1 { species1.weight } - , tags_1 { species1.tag } - , i1s_2 { species2.i1 } - , i2s_2 { species2.i2 } - , i3s_2 { species2.i3 } - , dx1s_2 { species2.dx1 } - , dx2s_2 { species2.dx2 } - , dx3s_2 { species2.dx3 } - , ux1s_2 { species2.ux1 } - , ux2s_2 { species2.ux2 } - , ux3s_2 { species2.ux3 } - , phis_2 { species2.phi } - , weights_2 { species2.weight } - , tags_2 { species2.tag } - , offset1 { offset1 } - , offset2 { offset2 } - , metric { metric } - , xi_min { xi_min } - , xi_max { xi_max } - , energy_dist_1 { energy_dist_1 } - , energy_dist_2 { energy_dist_2 } - , inv_V0 { inv_V0 } - , random_pool { random_pool } {} - - Inline void operator()(index_t p) const { - coord_t x_Cd { ZERO }; - vec_t v1 { ZERO }, v2 { ZERO }; - { // generate a random coordinate - auto rand_gen = random_pool.get_state(); - x_Cd[0] = xi_min(0) + Random(rand_gen) * (xi_max(0) - xi_min(0)); - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - x_Cd[1] = xi_min(1) + - Random(rand_gen) * (xi_max(1) - xi_min(1)); - } - if constexpr (M::Dim == Dim::_3D) { - x_Cd[2] = xi_min(2) + - Random(rand_gen) * (xi_max(2) - xi_min(2)); - } - random_pool.free_state(rand_gen); - } - { // generate the velocity - coord_t x_Ph { ZERO }; - metric.template convert(x_Cd, x_Ph); - if constexpr (M::CoordType == Coord::Cart) { - energy_dist_1(x_Ph, v1, spidx1); - energy_dist_2(x_Ph, v2, spidx2); - } else if constexpr (S == SimEngine::SRPIC) { - coord_t x_Cd_ { ZERO }; - x_Cd_[0] = x_Cd[0]; - x_Cd_[1] = x_Cd[1]; - x_Cd_[2] = ZERO; // phi = 0 - vec_t v_Ph { ZERO }; - energy_dist_1(x_Ph, v_Ph, spidx1); - metric.template transform_xyz(x_Cd_, v_Ph, v1); - energy_dist_2(x_Ph, v_Ph, spidx2); - metric.template transform_xyz(x_Cd_, v_Ph, v2); - } else if constexpr (S == SimEngine::GRPIC) { - vec_t v_Ph { ZERO }; - energy_dist_1(x_Ph, v_Ph, spidx1); - metric.template transform(x_Cd, v_Ph, v1); - energy_dist_2(x_Ph, v_Ph, spidx2); - metric.template transform(x_Cd, v_Ph, v2); - } else { - raise::KernelError(HERE, "Unknown simulation engine"); - } - } - // inject - i1s_1(p + offset1) = static_cast(x_Cd[0]); - dx1s_1(p + offset1) = static_cast( - x_Cd[0] - static_cast(i1s_1(p + offset1))); - i1s_2(p + offset2) = i1s_1(p + offset1); - dx1s_2(p + offset2) = dx1s_1(p + offset1); - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - i2s_1(p + offset1) = static_cast(x_Cd[1]); - dx2s_1(p + offset1) = static_cast( - x_Cd[1] - static_cast(i2s_1(p + offset1))); - i2s_2(p + offset2) = i2s_1(p + offset1); - dx2s_2(p + offset2) = dx2s_1(p + offset1); - if constexpr (S == SimEngine::SRPIC && M::CoordType != Coord::Cart) { - phis_1(p + offset1) = ZERO; - phis_2(p + offset2) = ZERO; - } - } - if constexpr (M::Dim == Dim::_3D) { - i3s_1(p + offset1) = static_cast(x_Cd[2]); - dx3s_1(p + offset1) = static_cast( - x_Cd[2] - static_cast(i3s_1(p + offset1))); - i3s_2(p + offset2) = i3s_1(p + offset1); - dx3s_2(p + offset2) = dx3s_1(p + offset1); - } - ux1s_1(p + offset1) = v1[0]; - ux2s_1(p + offset1) = v1[1]; - ux3s_1(p + offset1) = v1[2]; - ux1s_2(p + offset2) = v2[0]; - ux2s_2(p + offset2) = v2[1]; - ux3s_2(p + offset2) = v2[2]; - tags_1(p + offset1) = ParticleTag::alive; - tags_2(p + offset2) = ParticleTag::alive; - if constexpr (M::CoordType == Coord::Cart) { - weights_1(p + offset1) = ONE; - weights_2(p + offset2) = ONE; - } else { - const auto sqrt_det_h = metric.sqrt_det_h(x_Cd); - weights_1(p + offset1) = sqrt_det_h * inv_V0; - weights_2(p + offset2) = sqrt_det_h * inv_V0; - } - } - }; // struct UniformInjector_kernel - - template - struct NonUniformInjector_kernel { - static_assert(ED1::is_energy_dist, - "ED1 must be an energy distribution class"); - static_assert(ED2::is_energy_dist, - "ED2 must be an energy distribution class"); - static_assert(SD::is_spatial_dist, "SD must be a spatial distribution class"); - static_assert(M::is_metric, "M must be a metric class"); - - const real_t ppc0; - const spidx_t spidx1, spidx2; - - array_t i1s_1, i2s_1, i3s_1; - array_t dx1s_1, dx2s_1, dx3s_1; - array_t ux1s_1, ux2s_1, ux3s_1; - array_t phis_1; - array_t weights_1; - array_t tags_1; - - array_t i1s_2, i2s_2, i3s_2; - array_t dx1s_2, dx2s_2, dx3s_2; - array_t ux1s_2, ux2s_2, ux3s_2; - array_t phis_2; - array_t weights_2; - array_t tags_2; - - array_t idx { "idx" }; - - npart_t offset1, offset2; - M metric; - const ED1 energy_dist_1; - const ED2 energy_dist_2; - const SD spatial_dist; - const real_t inv_V0; - random_number_pool_t random_pool; - - NonUniformInjector_kernel(real_t ppc0, - spidx_t spidx1, - spidx_t spidx2, - Particles& species1, - Particles& species2, - npart_t offset1, - npart_t offset2, - const M& metric, - const ED1& energy_dist_1, - const ED2& energy_dist_2, - const SD& spatial_dist, - real_t inv_V0, - random_number_pool_t& random_pool) - : ppc0 { ppc0 } - , spidx1 { spidx1 } - , spidx2 { spidx2 } - , i1s_1 { species1.i1 } - , i2s_1 { species1.i2 } - , i3s_1 { species1.i3 } - , dx1s_1 { species1.dx1 } - , dx2s_1 { species1.dx2 } - , dx3s_1 { species1.dx3 } - , ux1s_1 { species1.ux1 } - , ux2s_1 { species1.ux2 } - , ux3s_1 { species1.ux3 } - , phis_1 { species1.phi } - , weights_1 { species1.weight } - , tags_1 { species1.tag } - , i1s_2 { species2.i1 } - , i2s_2 { species2.i2 } - , i3s_2 { species2.i3 } - , dx1s_2 { species2.dx1 } - , dx2s_2 { species2.dx2 } - , dx3s_2 { species2.dx3 } - , ux1s_2 { species2.ux1 } - , ux2s_2 { species2.ux2 } - , ux3s_2 { species2.ux3 } - , phis_2 { species2.phi } - , weights_2 { species2.weight } - , tags_2 { species2.tag } - , offset1 { offset1 } - , offset2 { offset2 } - , metric { metric } - , energy_dist_1 { energy_dist_1 } - , energy_dist_2 { energy_dist_2 } - , spatial_dist { spatial_dist } - , inv_V0 { inv_V0 } - , random_pool { random_pool } {} - - auto number_injected() const -> npart_t { - auto idx_h = Kokkos::create_mirror_view(idx); - Kokkos::deep_copy(idx_h, idx); - return idx_h(); - } - - Inline void operator()(index_t i1) const { - if constexpr (M::Dim == Dim::_1D) { - const auto i1_ = COORD(i1); - coord_t x_Cd { i1_ + HALF }; - coord_t x_Ph { ZERO }; - metric.template convert(x_Cd, x_Ph); - const auto ppc_real = static_cast(ppc0) * spatial_dist(x_Ph); - auto ppc = static_cast(ppc_real); - auto rand_gen = random_pool.get_state(); - if (Random(rand_gen) < (ppc_real - math::floor(ppc_real))) { - ++ppc; - } - for (auto p { 0u }; p < ppc; ++p) { - const auto index = Kokkos::atomic_fetch_add(&idx(), 1); - const auto dx1 = Random(rand_gen); - - i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; - dx1s_1(index + offset1) = dx1; - i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; - dx1s_2(index + offset2) = dx1; - - vec_t v_T { ZERO }, v_XYZ { ZERO }; - energy_dist_1(x_Ph, v_T, spidx1); - metric.template transform_xyz(x_Cd, v_T, v_XYZ); - ux1s_1(index + offset1) = v_XYZ[0]; - ux2s_1(index + offset1) = v_XYZ[1]; - ux3s_1(index + offset1) = v_XYZ[2]; - energy_dist_2(x_Ph, v_T, spidx2); - metric.template transform_xyz(x_Cd, v_T, v_XYZ); - ux1s_2(index + offset2) = v_XYZ[0]; - ux2s_2(index + offset2) = v_XYZ[1]; - ux3s_2(index + offset2) = v_XYZ[2]; - - tags_1(index + offset1) = ParticleTag::alive; - tags_2(index + offset2) = ParticleTag::alive; - if (M::CoordType == Coord::Cart) { - weights_1(index + offset1) = ONE; - weights_2(index + offset2) = ONE; - } else { - const auto wei = metric.sqrt_det_h({ i1_ + HALF }) * inv_V0; - weights_1(index + offset1) = wei; - weights_2(index + offset2) = wei; - } - } - random_pool.free_state(rand_gen); - } else { - raise::KernelError(HERE, "NonUniformInjector_kernel 1D called for 2D/3D"); - } - } - - Inline void operator()(index_t i1, index_t i2) const { - if constexpr (M::Dim == Dim::_2D) { - const auto i1_ = COORD(i1); - const auto i2_ = COORD(i2); - coord_t x_Cd { i1_ + HALF, i2_ + HALF }; - coord_t x_Ph { ZERO }; - coord_t x_Cd_ { ZERO }; - x_Cd_[0] = x_Cd[0]; - x_Cd_[1] = x_Cd[1]; - if constexpr (S == SimEngine::SRPIC and M::CoordType != Coord::Cart) { - x_Cd_[2] = ZERO; - } - metric.template convert(x_Cd, x_Ph); - const auto ppc_real = static_cast(ppc0) * spatial_dist(x_Ph); - auto ppc = static_cast(ppc_real); - auto rand_gen = random_pool.get_state(); - if (Random(rand_gen) < (ppc_real - math::floor(ppc_real))) { - ++ppc; - } - random_pool.free_state(rand_gen); - for (auto p { 0u }; p < ppc; ++p) { - const auto index = Kokkos::atomic_fetch_add(&idx(), 1); - auto rand_gen = random_pool.get_state(); - const auto dx1 = Random(rand_gen); - const auto dx2 = Random(rand_gen); - random_pool.free_state(rand_gen); - - i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; - dx1s_1(index + offset1) = dx1; - i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; - dx1s_2(index + offset2) = dx1; - - i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; - dx2s_1(index + offset1) = dx2; - i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; - dx2s_2(index + offset2) = dx2; - - vec_t v_T { ZERO }, v_Cd { ZERO }; - energy_dist_1(x_Ph, v_T, spidx1); - if constexpr (S == SimEngine::SRPIC) { - metric.template transform_xyz(x_Cd_, v_T, v_Cd); - } else if constexpr (S == SimEngine::GRPIC) { - metric.template transform(x_Cd_, v_T, v_Cd); - } - ux1s_1(index + offset1) = v_Cd[0]; - ux2s_1(index + offset1) = v_Cd[1]; - ux3s_1(index + offset1) = v_Cd[2]; - energy_dist_2(x_Ph, v_T, spidx2); - if constexpr (S == SimEngine::SRPIC) { - metric.template transform_xyz(x_Cd_, v_T, v_Cd); - } else if constexpr (S == SimEngine::GRPIC) { - metric.template transform(x_Cd_, v_T, v_Cd); - } - ux1s_2(index + offset2) = v_Cd[0]; - ux2s_2(index + offset2) = v_Cd[1]; - ux3s_2(index + offset2) = v_Cd[2]; - - tags_1(index + offset1) = ParticleTag::alive; - tags_2(index + offset2) = ParticleTag::alive; - if (M::CoordType == Coord::Cart) { - weights_1(index + offset1) = ONE; - weights_2(index + offset2) = ONE; - } else { - const auto wei = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF }) * - inv_V0; - weights_1(index + offset1) = wei; - weights_2(index + offset2) = wei; - } - } - } - - else { - raise::KernelError(HERE, "NonUniformInjector_kernel 2D called for 1D/3D"); - } - } - - Inline void operator()(index_t i1, index_t i2, index_t i3) const { - if constexpr (M::Dim == Dim::_3D) { - const auto i1_ = COORD(i1); - const auto i2_ = COORD(i2); - const auto i3_ = COORD(i3); - coord_t x_Cd { i1_ + HALF, i2_ + HALF, i3_ + HALF }; - coord_t x_Ph { ZERO }; - metric.template convert(x_Cd, x_Ph); - const auto ppc_real = static_cast(ppc0) * spatial_dist(x_Ph); - auto ppc = static_cast(ppc_real); - auto rand_gen = random_pool.get_state(); - if (Random(rand_gen) < (ppc_real - math::floor(ppc_real))) { - ++ppc; - } - for (auto p { 0u }; p < ppc; ++p) { - const auto index = Kokkos::atomic_fetch_add(&idx(), 1); - const auto dx1 = Random(rand_gen); - const auto dx2 = Random(rand_gen); - const auto dx3 = Random(rand_gen); - - i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; - dx1s_1(index + offset1) = dx1; - i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; - dx1s_2(index + offset2) = dx1; - - i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; - dx2s_1(index + offset1) = dx2; - i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; - dx2s_2(index + offset2) = dx2; - - i3s_1(index + offset1) = static_cast(i3) - N_GHOSTS; - dx3s_1(index + offset1) = dx3; - i3s_2(index + offset2) = static_cast(i3) - N_GHOSTS; - dx3s_2(index + offset2) = dx3; - - vec_t v_T { ZERO }, v_Cd { ZERO }; - energy_dist_1(x_Ph, v_T, spidx1); - if constexpr (S == SimEngine::SRPIC) { - metric.template transform_xyz(x_Cd, v_T, v_Cd); - } else if constexpr (S == SimEngine::GRPIC) { - metric.template transform(x_Cd, v_T, v_Cd); - } - ux1s_1(index + offset1) = v_Cd[0]; - ux2s_1(index + offset1) = v_Cd[1]; - ux3s_1(index + offset1) = v_Cd[2]; - energy_dist_2(x_Ph, v_T, spidx2); - if constexpr (S == SimEngine::SRPIC) { - metric.template transform_xyz(x_Cd, v_T, v_Cd); - } else if constexpr (S == SimEngine::GRPIC) { - metric.template transform(x_Cd, v_T, v_Cd); - } - ux1s_2(index + offset2) = v_Cd[0]; - ux2s_2(index + offset2) = v_Cd[1]; - ux3s_2(index + offset2) = v_Cd[2]; - - tags_1(index + offset1) = ParticleTag::alive; - tags_2(index + offset2) = ParticleTag::alive; - if (M::CoordType == Coord::Cart) { - weights_1(index + offset1) = ONE; - weights_2(index + offset2) = ONE; - } else { - const auto wei = metric.sqrt_det_h( - { i1_ + HALF, i2_ + HALF, i3_ + HALF }) * - inv_V0; - weights_1(index + offset1) = wei; - weights_2(index + offset2) = wei; - } - } - random_pool.free_state(rand_gen); - } else { - raise::KernelError(HERE, "NonUniformInjector_kernel 3D called for 1D/2D"); - } - } - }; // struct NonUniformInjector_kernel - - template - struct ConstNPPCInjector_kernel { - static_assert(ED1::is_energy_dist, - "ED1 must be an energy distribution class"); - static_assert(ED2::is_energy_dist, - "ED2 must be an energy distribution class"); - static_assert(M::is_metric, "M must be a metric class"); - - const real_t ppc0; - const spidx_t spidx1, spidx2; - - array_t i1s_1, i2s_1, i3s_1; - array_t dx1s_1, dx2s_1, dx3s_1; - array_t ux1s_1, ux2s_1, ux3s_1; - array_t phis_1; - array_t weights_1; - array_t tags_1; - - array_t i1s_2, i2s_2, i3s_2; - array_t dx1s_2, dx2s_2, dx3s_2; - array_t ux1s_2, ux2s_2, ux3s_2; - array_t phis_2; - array_t weights_2; - array_t tags_2; - - array_t idx { "idx" }; - - npart_t offset1, offset2; - M metric; - const ED1 energy_dist_1; - const ED2 energy_dist_2; - const ndfield_t nppc_array; - const idx_t nppc_idx; - const real_t inv_V0; - random_number_pool_t random_pool; - - ConstNPPCInjector_kernel(real_t ppc0, - spidx_t spidx1, - spidx_t spidx2, - Particles& species1, - Particles& species2, - npart_t offset1, - npart_t offset2, - const M& metric, - const ED1& energy_dist_1, - const ED2& energy_dist_2, - const ndfield_t& nppc_array, - idx_t nppc_idx, - real_t inv_V0, - random_number_pool_t& random_pool) - : ppc0 { ppc0 } - , spidx1 { spidx1 } - , spidx2 { spidx2 } - , i1s_1 { species1.i1 } - , i2s_1 { species1.i2 } - , i3s_1 { species1.i3 } - , dx1s_1 { species1.dx1 } - , dx2s_1 { species1.dx2 } - , dx3s_1 { species1.dx3 } - , ux1s_1 { species1.ux1 } - , ux2s_1 { species1.ux2 } - , ux3s_1 { species1.ux3 } - , phis_1 { species1.phi } - , weights_1 { species1.weight } - , tags_1 { species1.tag } - , i1s_2 { species2.i1 } - , i2s_2 { species2.i2 } - , i3s_2 { species2.i3 } - , dx1s_2 { species2.dx1 } - , dx2s_2 { species2.dx2 } - , dx3s_2 { species2.dx3 } - , ux1s_2 { species2.ux1 } - , ux2s_2 { species2.ux2 } - , ux3s_2 { species2.ux3 } - , phis_2 { species2.phi } - , weights_2 { species2.weight } - , tags_2 { species2.tag } - , offset1 { offset1 } - , offset2 { offset2 } - , metric { metric } - , energy_dist_1 { energy_dist_1 } - , energy_dist_2 { energy_dist_2 } // , spatial_dist { spatial_dist } - , nppc_array { nppc_array } - , nppc_idx { nppc_idx } - , inv_V0 { inv_V0 } - , random_pool { random_pool } {} - - auto number_injected() const -> npart_t { - auto idx_h = Kokkos::create_mirror_view(idx); - Kokkos::deep_copy(idx_h, idx); - return idx_h(); - } - - Inline auto get_inj_ppc(real_t& ppc_real) const -> npart_t { - return static_cast( - math::max(ppc0 - ppc_real, static_cast(0.0))); - } - - Inline void operator()(index_t i1) const { - if constexpr (M::Dim == Dim::_1D) { - const auto ppc = get_inj_ppc(nppc_array(i1, nppc_idx)); - const auto i1_ = COORD(i1); - for (auto p { 0u }; p < ppc; ++p) { - const auto index = Kokkos::atomic_fetch_add(&idx(), 1); - auto rand_gen = random_pool.get_state(); - const auto dx1 = Random(rand_gen); - random_pool.free_state(rand_gen); - - coord_t x_Cd { i1_ + dx1 }; - coord_t x_Ph { ZERO }; - metric.template convert(x_Cd, x_Ph); - - i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; - dx1s_1(index + offset1) = dx1; - i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; - dx1s_2(index + offset2) = dx1; - - vec_t v_T { ZERO }, v_XYZ { ZERO }; - energy_dist_1(x_Ph, v_T, spidx1); - metric.template transform_xyz(x_Cd, v_T, v_XYZ); - ux1s_1(index + offset1) = v_XYZ[0]; - ux2s_1(index + offset1) = v_XYZ[1]; - ux3s_1(index + offset1) = v_XYZ[2]; - energy_dist_2(x_Ph, v_T, spidx2); - metric.template transform_xyz(x_Cd, v_T, v_XYZ); - ux1s_2(index + offset2) = v_XYZ[0]; - ux2s_2(index + offset2) = v_XYZ[1]; - ux3s_2(index + offset2) = v_XYZ[2]; - - tags_1(index + offset1) = ParticleTag::alive; - tags_2(index + offset2) = ParticleTag::alive; - if (M::CoordType == Coord::Cart) { - weights_1(index + offset1) = ONE; - weights_2(index + offset2) = ONE; - } else { - const auto wei = metric.sqrt_det_h({ i1_ + HALF }) * inv_V0; - weights_1(index + offset1) = wei; - weights_2(index + offset2) = wei; - } - } - } else { - raise::KernelError(HERE, "NonUniformInjector_kernel 1D called for 2D/3D"); - } - } - - Inline void operator()(index_t i1, index_t i2) const { - if constexpr (M::Dim == Dim::_2D) { - const auto ppc = get_inj_ppc(nppc_array(i1, i2, nppc_idx)); - const auto i1_ = COORD(i1); - const auto i2_ = COORD(i2); - for (auto p { 0u }; p < ppc; ++p) { - const auto index = Kokkos::atomic_fetch_add(&idx(), 1); - auto rand_gen = random_pool.get_state(); - const auto dx1 = Random(rand_gen); - const auto dx2 = Random(rand_gen); - random_pool.free_state(rand_gen); - - coord_t x_Cd { i1_ + dx1, i2_ + dx2 }; - coord_t x_Ph { ZERO }; - coord_t x_Cd_ { ZERO }; - x_Cd_[0] = x_Cd[0]; - x_Cd_[1] = x_Cd[1]; - if constexpr (S == SimEngine::SRPIC and M::CoordType != Coord::Cart) { - x_Cd_[2] = ZERO; - } - metric.template convert(x_Cd, x_Ph); - - i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; - dx1s_1(index + offset1) = dx1; - i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; - dx1s_2(index + offset2) = dx1; - - i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; - dx2s_1(index + offset1) = dx2; - i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; - dx2s_2(index + offset2) = dx2; - - vec_t v_T { ZERO }, v_Cd { ZERO }; - energy_dist_1(x_Ph, v_T, spidx1); - if constexpr (S == SimEngine::SRPIC) { - metric.template transform_xyz(x_Cd_, v_T, v_Cd); - } else if constexpr (S == SimEngine::GRPIC) { - metric.template transform(x_Cd_, v_T, v_Cd); - } - ux1s_1(index + offset1) = v_Cd[0]; - ux2s_1(index + offset1) = v_Cd[1]; - ux3s_1(index + offset1) = v_Cd[2]; - energy_dist_2(x_Ph, v_T, spidx2); - if constexpr (S == SimEngine::SRPIC) { - metric.template transform_xyz(x_Cd_, v_T, v_Cd); - } else if constexpr (S == SimEngine::GRPIC) { - metric.template transform(x_Cd_, v_T, v_Cd); - } - ux1s_2(index + offset2) = v_Cd[0]; - ux2s_2(index + offset2) = v_Cd[1]; - ux3s_2(index + offset2) = v_Cd[2]; - - tags_1(index + offset1) = ParticleTag::alive; - tags_2(index + offset2) = ParticleTag::alive; - if (M::CoordType == Coord::Cart) { - weights_1(index + offset1) = ONE; - weights_2(index + offset2) = ONE; - } else { - const auto wei = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF }) * - inv_V0; - weights_1(index + offset1) = wei; - weights_2(index + offset2) = wei; - } - } - } - - else { - raise::KernelError(HERE, "NonUniformInjector_kernel 2D called for 1D/3D"); - } - } - - Inline void operator()(index_t i1, index_t i2, index_t i3) const { - if constexpr (M::Dim == Dim::_3D) { - const auto ppc = get_inj_ppc(nppc_array(i1, i2, i3, nppc_idx)); - const auto i1_ = COORD(i1); - const auto i2_ = COORD(i2); - const auto i3_ = COORD(i3); - for (auto p { 0u }; p < ppc; ++p) { - const auto index = Kokkos::atomic_fetch_add(&idx(), 1); - auto rand_gen = random_pool.get_state(); - const auto dx1 = Random(rand_gen); - const auto dx2 = Random(rand_gen); - const auto dx3 = Random(rand_gen); - random_pool.free_state(rand_gen); - - coord_t x_Cd { i1_ + dx1, i2_ + dx2, i3_ + dx3 }; - coord_t x_Ph { ZERO }; - coord_t x_Cd_ { ZERO }; - x_Cd_[0] = x_Cd[0]; - x_Cd_[1] = x_Cd[1]; - x_Cd_[2] = x_Cd[2]; - metric.template convert(x_Cd, x_Ph); - - i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; - dx1s_1(index + offset1) = dx1; - i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; - dx1s_2(index + offset2) = dx1; - - i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; - dx2s_1(index + offset1) = dx2; - i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; - dx2s_2(index + offset2) = dx2; - - i3s_1(index + offset1) = static_cast(i3) - N_GHOSTS; - dx3s_1(index + offset1) = dx3; - i3s_2(index + offset2) = static_cast(i3) - N_GHOSTS; - dx3s_2(index + offset2) = dx3; - - vec_t v_T { ZERO }, v_Cd { ZERO }; - energy_dist_1(x_Ph, v_T, spidx1); - if constexpr (S == SimEngine::SRPIC) { - metric.template transform_xyz(x_Cd, v_T, v_Cd); - } else if constexpr (S == SimEngine::GRPIC) { - metric.template transform(x_Cd, v_T, v_Cd); - } - ux1s_1(index + offset1) = v_Cd[0]; - ux2s_1(index + offset1) = v_Cd[1]; - ux3s_1(index + offset1) = v_Cd[2]; - energy_dist_2(x_Ph, v_T, spidx2); - if constexpr (S == SimEngine::SRPIC) { - metric.template transform_xyz(x_Cd, v_T, v_Cd); - } else if constexpr (S == SimEngine::GRPIC) { - metric.template transform(x_Cd, v_T, v_Cd); - } - ux1s_2(index + offset2) = v_Cd[0]; - ux2s_2(index + offset2) = v_Cd[1]; - ux3s_2(index + offset2) = v_Cd[2]; - - tags_1(index + offset1) = ParticleTag::alive; - tags_2(index + offset2) = ParticleTag::alive; - if (M::CoordType == Coord::Cart) { - weights_1(index + offset1) = ONE; - weights_2(index + offset2) = ONE; - } else { - const auto wei = metric.sqrt_det_h( - { i1_ + HALF, i2_ + HALF, i3_ + HALF }) * - inv_V0; - weights_1(index + offset1) = wei; - weights_2(index + offset2) = wei; - } - } - } else { - raise::KernelError(HERE, "NonUniformInjector_kernel 3D called for 1D/2D"); - } - } - }; // struct NonUniformInjector_kernel - - } // namespace experimental - - template - struct GlobalInjector_kernel { - static_assert(M::is_metric, "M must be a metric class"); - static constexpr auto D = M::Dim; - - const bool use_weights; - - array_t in_ux1; - array_t in_ux2; - array_t in_ux3; - array_t in_x1; - array_t in_x2; - array_t in_x3; - array_t in_phi; - array_t in_wei; - - array_t idx { "idx" }; - array_t i1s, i2s, i3s; - array_t dx1s, dx2s, dx3s; - array_t ux1s, ux2s, ux3s; - array_t phis; - array_t weights; - array_t tags; - - const npart_t offset; - - M global_metric; - - real_t x1_min, x1_max, x2_min, x2_max, x3_min, x3_max; - ncells_t i1_offset, i2_offset, i3_offset; - - GlobalInjector_kernel(Particles& species, - const M& global_metric, - const Domain& local_domain, - const std::map>& data, - bool use_weights) - : use_weights { use_weights } - , i1s { species.i1 } - , i2s { species.i2 } - , i3s { species.i3 } - , dx1s { species.dx1 } - , dx2s { species.dx2 } - , dx3s { species.dx3 } - , ux1s { species.ux1 } - , ux2s { species.ux2 } - , ux3s { species.ux3 } - , phis { species.phi } - , weights { species.weight } - , tags { species.tag } - , offset { species.npart() } - , global_metric { global_metric } { - const auto n_inject = data.at("x1").size(); - - x1_min = local_domain.mesh.extent(in::x1).first; - x1_max = local_domain.mesh.extent(in::x1).second; - i1_offset = local_domain.offset_ncells()[0]; - - copy_from_vector("x1", in_x1, data, n_inject); - copy_from_vector("ux1", in_ux1, data, n_inject); - copy_from_vector("ux2", in_ux2, data, n_inject); - copy_from_vector("ux3", in_ux3, data, n_inject); - if (use_weights) { - copy_from_vector("weights", in_wei, data, n_inject); - } - if constexpr (D == Dim::_2D or D == Dim::_3D) { - x2_min = local_domain.mesh.extent(in::x2).first; - x2_max = local_domain.mesh.extent(in::x2).second; - i2_offset = local_domain.offset_ncells()[1]; - copy_from_vector("x2", in_x2, data, n_inject); - } - if constexpr (D == Dim::_2D and M::CoordType != Coord::Cart) { - copy_from_vector("phi", in_phi, data, n_inject); - } - if constexpr (D == Dim::_3D) { - x3_min = local_domain.mesh.extent(in::x3).first; - x3_max = local_domain.mesh.extent(in::x3).second; - i3_offset = local_domain.offset_ncells()[2]; - copy_from_vector("x3", in_x3, data, n_inject); - } - } - - void copy_from_vector(const std::string& name, - array_t& arr, - const std::map>& data, - npart_t n_inject) { - raise::ErrorIf(data.find(name) == data.end(), name + " not found in data", HERE); - raise::ErrorIf(data.at(name).size() != n_inject, "Inconsistent data size", HERE); - arr = array_t { name, n_inject }; - auto arr_h = Kokkos::create_mirror_view(arr); - for (auto i = 0u; i < data.at(name).size(); ++i) { - arr_h(i) = data.at(name)[i]; - } - Kokkos::deep_copy(arr, arr_h); - } - - auto number_injected() const -> npart_t { - auto idx_h = Kokkos::create_mirror_view(idx); - Kokkos::deep_copy(idx_h, idx); - return idx_h(); - } - - Inline void operator()(index_t p) const { - if constexpr (D == Dim::_1D) { - if (in_x1(p) >= x1_min and in_x1(p) < x1_max) { - coord_t x_Cd { ZERO }; - vec_t u_XYZ { ZERO }; - const vec_t u_Ph { in_ux1(p), in_ux2(p), in_ux3(p) }; - - auto index { offset + Kokkos::atomic_fetch_add(&idx(), 1) }; - global_metric.template convert({ in_x1(p) }, x_Cd); - global_metric.template transform_xyz(x_Cd, u_Ph, u_XYZ); - - const auto i1 = static_cast( - static_cast(x_Cd[0]) - i1_offset); - const auto dx1 = static_cast( - x_Cd[0] - static_cast(i1 + i1_offset)); - - i1s(index) = i1; - dx1s(index) = dx1; - ux1s(index) = u_XYZ[0]; - ux2s(index) = u_XYZ[1]; - ux3s(index) = u_XYZ[2]; - tags(index) = ParticleTag::alive; - if (use_weights) { - weights(index) = weights(p); - } else { - weights(index) = ONE; - } - } - } else if constexpr (D == Dim::_2D) { - if ((in_x1(p) >= x1_min and in_x1(p) < x1_max) and - (in_x2(p) >= x2_min and in_x2(p) < x2_max)) { - coord_t x_Cd { ZERO }; - vec_t u_Cd { ZERO }; - vec_t u_Ph { in_ux1(p), in_ux2(p), in_ux3(p) }; - coord_t x_Cd_ { ZERO }; - - auto index { offset + - Kokkos::atomic_fetch_add(&idx(), static_cast(1)) }; - global_metric.template convert({ in_x1(p), in_x2(p) }, - x_Cd); - x_Cd_[0] = x_Cd[0]; - x_Cd_[1] = x_Cd[1]; - if constexpr (S == SimEngine::SRPIC and M::CoordType != Coord::Cart) { - x_Cd_[2] = in_phi(p); - } - if constexpr (S == SimEngine::SRPIC) { - global_metric.template transform_xyz(x_Cd_, u_Ph, u_Cd); - } else if constexpr (S == SimEngine::GRPIC) { - global_metric.template transform(x_Cd, u_Ph, u_Cd); - } else { - raise::KernelError(HERE, "Unknown simulation engine"); - } - const auto i1 = static_cast( - static_cast(x_Cd[0]) - i1_offset); - const auto dx1 = static_cast( - x_Cd[0] - static_cast(i1 + i1_offset)); - const auto i2 = static_cast( - static_cast(x_Cd[1]) - i2_offset); - const auto dx2 = static_cast( - x_Cd[1] - static_cast(i2 + i2_offset)); - - i1s(index) = i1; - dx1s(index) = dx1; - i2s(index) = i2; - dx2s(index) = dx2; - ux1s(index) = u_Cd[0]; - ux2s(index) = u_Cd[1]; - ux3s(index) = u_Cd[2]; - if (M::CoordType != Coord::Cart) { - phis(index) = in_phi(p); - } - tags(index) = ParticleTag::alive; - if (use_weights) { - weights(index) = weights(p); - } else { - weights(index) = ONE; - } - } - } else { - if ((in_x1(p) >= x1_min and in_x1(p) < x1_max) and - (in_x2(p) >= x2_min and in_x2(p) < x2_max) and - (in_x3(p) >= x3_min and in_x3(p) < x3_max)) { - coord_t x_Cd { ZERO }; - vec_t u_Cd { ZERO }; - vec_t u_Ph { in_ux1(p), in_ux2(p), in_ux3(p) }; - - auto index { offset + Kokkos::atomic_fetch_add(&idx(), 1) }; - global_metric.template convert( - { in_x1(p), in_x2(p), in_x3(p) }, - x_Cd); - if constexpr (S == SimEngine::SRPIC) { - global_metric.template transform_xyz(x_Cd, u_Ph, u_Cd); - } else if constexpr (S == SimEngine::GRPIC) { - global_metric.template transform(x_Cd, u_Ph, u_Cd); - } else { - raise::KernelError(HERE, "Unknown simulation engine"); - } - const auto i1 = static_cast( - static_cast(x_Cd[0]) - i1_offset); - const auto dx1 = static_cast( - x_Cd[0] - static_cast(i1 + i1_offset)); - const auto i2 = static_cast( - static_cast(x_Cd[1]) - i2_offset); - const auto dx2 = static_cast( - x_Cd[1] - static_cast(i2 + i2_offset)); - const auto i3 = static_cast( - static_cast(x_Cd[2]) - i3_offset); - const auto dx3 = static_cast( - x_Cd[2] - static_cast(i3 + i3_offset)); - - i1s(index) = i1; - dx1s(index) = dx1; - i2s(index) = i2; - dx2s(index) = dx2; - i3s(index) = i3; - dx3s(index) = dx3; - ux1s(index) = u_Cd[0]; - ux2s(index) = u_Cd[1]; - ux3s(index) = u_Cd[2]; - tags(index) = ParticleTag::alive; - if (use_weights) { - weights(index) = weights(p); - } else { - weights(index) = ONE; - } - } - } - } - }; // struct GlobalInjector_kernel - - template + template struct NonUniformInjector_kernel { - static_assert(ED::is_energy_dist, "ED must be an energy distribution class"); + static_assert(ED1::is_energy_dist, "ED1 must be an energy distribution class"); + static_assert(ED2::is_energy_dist, "ED2 must be an energy distribution class"); static_assert(SD::is_spatial_dist, "SD must be a spatial distribution class"); static_assert(M::is_metric, "M must be a metric class"); @@ -1195,7 +216,8 @@ namespace kernel { npart_t offset1, offset2; M metric; - const ED energy_dist; + const ED1 energy_dist_1; + const ED2 energy_dist_2; const SD spatial_dist; const real_t inv_V0; random_number_pool_t random_pool; @@ -1208,7 +230,8 @@ namespace kernel { npart_t offset1, npart_t offset2, const M& metric, - const ED& energy_dist, + const ED1& energy_dist_1, + const ED2& energy_dist_2, const SD& spatial_dist, real_t inv_V0, random_number_pool_t& random_pool) @@ -1242,7 +265,8 @@ namespace kernel { , offset1 { offset1 } , offset2 { offset2 } , metric { metric } - , energy_dist { energy_dist } + , energy_dist_1 { energy_dist_1 } + , energy_dist_2 { energy_dist_2 } , spatial_dist { spatial_dist } , inv_V0 { inv_V0 } , random_pool { random_pool } {} @@ -1259,11 +283,12 @@ namespace kernel { coord_t x_Cd { i1_ + HALF }; coord_t x_Ph { ZERO }; metric.template convert(x_Cd, x_Ph); - const auto ppc = static_cast(ppc0 * spatial_dist(x_Ph)); - if (ppc == 0) { - return; + const auto ppc_real = static_cast(ppc0) * spatial_dist(x_Ph); + auto ppc = static_cast(ppc_real); + auto rand_gen = random_pool.get_state(); + if (Random(rand_gen) < (ppc_real - math::floor(ppc_real))) { + ++ppc; } - auto rand_gen = random_pool.get_state(); for (auto p { 0u }; p < ppc; ++p) { const auto index = Kokkos::atomic_fetch_add(&idx(), 1); const auto dx1 = Random(rand_gen); @@ -1274,12 +299,12 @@ namespace kernel { dx1s_2(index + offset2) = dx1; vec_t v_T { ZERO }, v_XYZ { ZERO }; - energy_dist(x_Ph, v_T, spidx1); + energy_dist_1(x_Ph, v_T, spidx1); metric.template transform_xyz(x_Cd, v_T, v_XYZ); ux1s_1(index + offset1) = v_XYZ[0]; ux2s_1(index + offset1) = v_XYZ[1]; ux3s_1(index + offset1) = v_XYZ[2]; - energy_dist(x_Ph, v_T, spidx2); + energy_dist_2(x_Ph, v_T, spidx2); metric.template transform_xyz(x_Cd, v_T, v_XYZ); ux1s_2(index + offset2) = v_XYZ[0]; ux2s_2(index + offset2) = v_XYZ[1]; @@ -1315,15 +340,307 @@ namespace kernel { x_Cd_[2] = ZERO; } metric.template convert(x_Cd, x_Ph); - const auto ppc = static_cast(ppc0 * spatial_dist(x_Ph)); - if (ppc == 0) { - return; + const auto ppc_real = static_cast(ppc0) * spatial_dist(x_Ph); + auto ppc = static_cast(ppc_real); + auto rand_gen = random_pool.get_state(); + if (Random(rand_gen) < (ppc_real - math::floor(ppc_real))) { + ++ppc; + } + random_pool.free_state(rand_gen); + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + auto rand_gen = random_pool.get_state(); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); + random_pool.free_state(rand_gen); + + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; + + i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; + dx2s_1(index + offset1) = dx2; + i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; + dx2s_2(index + offset2) = dx2; + + vec_t v_T { ZERO }, v_Cd { ZERO }; + energy_dist_1(x_Ph, v_T, spidx1); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd_, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd_, v_T, v_Cd); + } + ux1s_1(index + offset1) = v_Cd[0]; + ux2s_1(index + offset1) = v_Cd[1]; + ux3s_1(index + offset1) = v_Cd[2]; + energy_dist_2(x_Ph, v_T, spidx2); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd_, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd_, v_T, v_Cd); + } + ux1s_2(index + offset2) = v_Cd[0]; + ux2s_2(index + offset2) = v_Cd[1]; + ux3s_2(index + offset2) = v_Cd[2]; + + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; + if (M::CoordType == Coord::Cart) { + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; + } else { + const auto wei = metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF }) * inv_V0; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; + } + } + } + + else { + raise::KernelError(HERE, "NonUniformInjector_kernel 2D called for 1D/3D"); + } + } + + Inline void operator()(index_t i1, index_t i2, index_t i3) const { + if constexpr (M::Dim == Dim::_3D) { + const auto i1_ = COORD(i1); + const auto i2_ = COORD(i2); + const auto i3_ = COORD(i3); + coord_t x_Cd { i1_ + HALF, i2_ + HALF, i3_ + HALF }; + coord_t x_Ph { ZERO }; + metric.template convert(x_Cd, x_Ph); + const auto ppc_real = static_cast(ppc0) * spatial_dist(x_Ph); + auto ppc = static_cast(ppc_real); + auto rand_gen = random_pool.get_state(); + if (Random(rand_gen) < (ppc_real - math::floor(ppc_real))) { + ++ppc; } - auto rand_gen = random_pool.get_state(); for (auto p { 0u }; p < ppc; ++p) { const auto index = Kokkos::atomic_fetch_add(&idx(), 1); const auto dx1 = Random(rand_gen); const auto dx2 = Random(rand_gen); + const auto dx3 = Random(rand_gen); + + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; + + i2s_1(index + offset1) = static_cast(i2) - N_GHOSTS; + dx2s_1(index + offset1) = dx2; + i2s_2(index + offset2) = static_cast(i2) - N_GHOSTS; + dx2s_2(index + offset2) = dx2; + + i3s_1(index + offset1) = static_cast(i3) - N_GHOSTS; + dx3s_1(index + offset1) = dx3; + i3s_2(index + offset2) = static_cast(i3) - N_GHOSTS; + dx3s_2(index + offset2) = dx3; + + vec_t v_T { ZERO }, v_Cd { ZERO }; + energy_dist_1(x_Ph, v_T, spidx1); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd, v_T, v_Cd); + } + ux1s_1(index + offset1) = v_Cd[0]; + ux2s_1(index + offset1) = v_Cd[1]; + ux3s_1(index + offset1) = v_Cd[2]; + energy_dist_2(x_Ph, v_T, spidx2); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd, v_T, v_Cd); + } + ux1s_2(index + offset2) = v_Cd[0]; + ux2s_2(index + offset2) = v_Cd[1]; + ux3s_2(index + offset2) = v_Cd[2]; + + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; + if (M::CoordType == Coord::Cart) { + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; + } else { + const auto wei = metric.sqrt_det_h( + { i1_ + HALF, i2_ + HALF, i3_ + HALF }) * + inv_V0; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; + } + } + random_pool.free_state(rand_gen); + } else { + raise::KernelError(HERE, "NonUniformInjector_kernel 3D called for 1D/2D"); + } + } + }; // struct NonUniformInjector_kernel + + template + struct ConstNPPCInjector_kernel { + static_assert(ED1::is_energy_dist, "ED1 must be an energy distribution class"); + static_assert(ED2::is_energy_dist, "ED2 must be an energy distribution class"); + static_assert(M::is_metric, "M must be a metric class"); + + const real_t ppc0; + const spidx_t spidx1, spidx2; + + array_t i1s_1, i2s_1, i3s_1; + array_t dx1s_1, dx2s_1, dx3s_1; + array_t ux1s_1, ux2s_1, ux3s_1; + array_t phis_1; + array_t weights_1; + array_t tags_1; + + array_t i1s_2, i2s_2, i3s_2; + array_t dx1s_2, dx2s_2, dx3s_2; + array_t ux1s_2, ux2s_2, ux3s_2; + array_t phis_2; + array_t weights_2; + array_t tags_2; + + array_t idx { "idx" }; + + npart_t offset1, offset2; + M metric; + const ED1 energy_dist_1; + const ED2 energy_dist_2; + const ndfield_t nppc_array; + const idx_t nppc_idx; + const real_t inv_V0; + random_number_pool_t random_pool; + + ConstNPPCInjector_kernel(real_t ppc0, + spidx_t spidx1, + spidx_t spidx2, + Particles& species1, + Particles& species2, + npart_t offset1, + npart_t offset2, + const M& metric, + const ED1& energy_dist_1, + const ED2& energy_dist_2, + const ndfield_t& nppc_array, + idx_t nppc_idx, + real_t inv_V0, + random_number_pool_t& random_pool) + : ppc0 { ppc0 } + , spidx1 { spidx1 } + , spidx2 { spidx2 } + , i1s_1 { species1.i1 } + , i2s_1 { species1.i2 } + , i3s_1 { species1.i3 } + , dx1s_1 { species1.dx1 } + , dx2s_1 { species1.dx2 } + , dx3s_1 { species1.dx3 } + , ux1s_1 { species1.ux1 } + , ux2s_1 { species1.ux2 } + , ux3s_1 { species1.ux3 } + , phis_1 { species1.phi } + , weights_1 { species1.weight } + , tags_1 { species1.tag } + , i1s_2 { species2.i1 } + , i2s_2 { species2.i2 } + , i3s_2 { species2.i3 } + , dx1s_2 { species2.dx1 } + , dx2s_2 { species2.dx2 } + , dx3s_2 { species2.dx3 } + , ux1s_2 { species2.ux1 } + , ux2s_2 { species2.ux2 } + , ux3s_2 { species2.ux3 } + , phis_2 { species2.phi } + , weights_2 { species2.weight } + , tags_2 { species2.tag } + , offset1 { offset1 } + , offset2 { offset2 } + , metric { metric } + , energy_dist_1 { energy_dist_1 } + , energy_dist_2 { energy_dist_2 } // , spatial_dist { spatial_dist } + , nppc_array { nppc_array } + , nppc_idx { nppc_idx } + , inv_V0 { inv_V0 } + , random_pool { random_pool } {} + + auto number_injected() const -> npart_t { + auto idx_h = Kokkos::create_mirror_view(idx); + Kokkos::deep_copy(idx_h, idx); + return idx_h(); + } + + Inline auto get_inj_ppc(real_t& ppc_real) const -> npart_t { + return static_cast( + math::max(ppc0 - ppc_real, static_cast(0.0))); + } + + Inline void operator()(index_t i1) const { + if constexpr (M::Dim == Dim::_1D) { + const auto ppc = get_inj_ppc(nppc_array(i1, nppc_idx)); + const auto i1_ = COORD(i1); + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + auto rand_gen = random_pool.get_state(); + const auto dx1 = Random(rand_gen); + random_pool.free_state(rand_gen); + + coord_t x_Cd { i1_ + dx1 }; + coord_t x_Ph { ZERO }; + metric.template convert(x_Cd, x_Ph); + + i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; + dx1s_1(index + offset1) = dx1; + i1s_2(index + offset2) = static_cast(i1) - N_GHOSTS; + dx1s_2(index + offset2) = dx1; + + vec_t v_T { ZERO }, v_XYZ { ZERO }; + energy_dist_1(x_Ph, v_T, spidx1); + metric.template transform_xyz(x_Cd, v_T, v_XYZ); + ux1s_1(index + offset1) = v_XYZ[0]; + ux2s_1(index + offset1) = v_XYZ[1]; + ux3s_1(index + offset1) = v_XYZ[2]; + energy_dist_2(x_Ph, v_T, spidx2); + metric.template transform_xyz(x_Cd, v_T, v_XYZ); + ux1s_2(index + offset2) = v_XYZ[0]; + ux2s_2(index + offset2) = v_XYZ[1]; + ux3s_2(index + offset2) = v_XYZ[2]; + + tags_1(index + offset1) = ParticleTag::alive; + tags_2(index + offset2) = ParticleTag::alive; + if (M::CoordType == Coord::Cart) { + weights_1(index + offset1) = ONE; + weights_2(index + offset2) = ONE; + } else { + const auto wei = metric.sqrt_det_h({ i1_ + HALF }) * inv_V0; + weights_1(index + offset1) = wei; + weights_2(index + offset2) = wei; + } + } + } else { + raise::KernelError(HERE, "NonUniformInjector_kernel 1D called for 2D/3D"); + } + } + + Inline void operator()(index_t i1, index_t i2) const { + if constexpr (M::Dim == Dim::_2D) { + const auto ppc = get_inj_ppc(nppc_array(i1, i2, nppc_idx)); + const auto i1_ = COORD(i1); + const auto i2_ = COORD(i2); + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + auto rand_gen = random_pool.get_state(); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); + random_pool.free_state(rand_gen); + + coord_t x_Cd { i1_ + dx1, i2_ + dx2 }; + coord_t x_Ph { ZERO }; + coord_t x_Cd_ { ZERO }; + x_Cd_[0] = x_Cd[0]; + x_Cd_[1] = x_Cd[1]; + if constexpr (S == SimEngine::SRPIC and M::CoordType != Coord::Cart) { + x_Cd_[2] = ZERO; + } + metric.template convert(x_Cd, x_Ph); i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; dx1s_1(index + offset1) = dx1; @@ -1336,7 +653,7 @@ namespace kernel { dx2s_2(index + offset2) = dx2; vec_t v_T { ZERO }, v_Cd { ZERO }; - energy_dist(x_Ph, v_T, spidx1); + energy_dist_1(x_Ph, v_T, spidx1); if constexpr (S == SimEngine::SRPIC) { metric.template transform_xyz(x_Cd_, v_T, v_Cd); } else if constexpr (S == SimEngine::GRPIC) { @@ -1345,7 +662,7 @@ namespace kernel { ux1s_1(index + offset1) = v_Cd[0]; ux2s_1(index + offset1) = v_Cd[1]; ux3s_1(index + offset1) = v_Cd[2]; - energy_dist(x_Ph, v_T, spidx2); + energy_dist_2(x_Ph, v_T, spidx2); if constexpr (S == SimEngine::SRPIC) { metric.template transform_xyz(x_Cd_, v_T, v_Cd); } else if constexpr (S == SimEngine::GRPIC) { @@ -1366,7 +683,6 @@ namespace kernel { weights_2(index + offset2) = wei; } } - random_pool.free_state(rand_gen); } else { @@ -1376,22 +692,25 @@ namespace kernel { Inline void operator()(index_t i1, index_t i2, index_t i3) const { if constexpr (M::Dim == Dim::_3D) { - const auto i1_ = COORD(i1); - const auto i2_ = COORD(i2); - const auto i3_ = COORD(i3); - coord_t x_Cd { i1_ + HALF, i2_ + HALF, i3_ + HALF }; - coord_t x_Ph { ZERO }; - metric.template convert(x_Cd, x_Ph); - const auto ppc = static_cast(ppc0 * spatial_dist(x_Ph)); - if (ppc == 0) { - return; - } - auto rand_gen = random_pool.get_state(); + const auto ppc = get_inj_ppc(nppc_array(i1, i2, i3, nppc_idx)); + const auto i1_ = COORD(i1); + const auto i2_ = COORD(i2); + const auto i3_ = COORD(i3); for (auto p { 0u }; p < ppc; ++p) { - const auto index = Kokkos::atomic_fetch_add(&idx(), 1); - const auto dx1 = Random(rand_gen); - const auto dx2 = Random(rand_gen); - const auto dx3 = Random(rand_gen); + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + auto rand_gen = random_pool.get_state(); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); + const auto dx3 = Random(rand_gen); + random_pool.free_state(rand_gen); + + coord_t x_Cd { i1_ + dx1, i2_ + dx2, i3_ + dx3 }; + coord_t x_Ph { ZERO }; + coord_t x_Cd_ { ZERO }; + x_Cd_[0] = x_Cd[0]; + x_Cd_[1] = x_Cd[1]; + x_Cd_[2] = x_Cd[2]; + metric.template convert(x_Cd, x_Ph); i1s_1(index + offset1) = static_cast(i1) - N_GHOSTS; dx1s_1(index + offset1) = dx1; @@ -1409,7 +728,7 @@ namespace kernel { dx3s_2(index + offset2) = dx3; vec_t v_T { ZERO }, v_Cd { ZERO }; - energy_dist(x_Ph, v_T, spidx1); + energy_dist_1(x_Ph, v_T, spidx1); if constexpr (S == SimEngine::SRPIC) { metric.template transform_xyz(x_Cd, v_T, v_Cd); } else if constexpr (S == SimEngine::GRPIC) { @@ -1418,7 +737,7 @@ namespace kernel { ux1s_1(index + offset1) = v_Cd[0]; ux2s_1(index + offset1) = v_Cd[1]; ux3s_1(index + offset1) = v_Cd[2]; - energy_dist(x_Ph, v_T, spidx2); + energy_dist_2(x_Ph, v_T, spidx2); if constexpr (S == SimEngine::SRPIC) { metric.template transform_xyz(x_Cd, v_T, v_Cd); } else if constexpr (S == SimEngine::GRPIC) { @@ -1441,13 +760,243 @@ namespace kernel { weights_2(index + offset2) = wei; } } - random_pool.free_state(rand_gen); } else { raise::KernelError(HERE, "NonUniformInjector_kernel 3D called for 1D/2D"); } } }; // struct NonUniformInjector_kernel + template + struct GlobalInjector_kernel { + static_assert(M::is_metric, "M must be a metric class"); + static constexpr auto D = M::Dim; + + const bool use_weights; + + array_t in_ux1; + array_t in_ux2; + array_t in_ux3; + array_t in_x1; + array_t in_x2; + array_t in_x3; + array_t in_phi; + array_t in_wei; + + array_t idx { "idx" }; + array_t i1s, i2s, i3s; + array_t dx1s, dx2s, dx3s; + array_t ux1s, ux2s, ux3s; + array_t phis; + array_t weights; + array_t tags; + + const npart_t offset; + + M global_metric; + + real_t x1_min, x1_max, x2_min, x2_max, x3_min, x3_max; + ncells_t i1_offset, i2_offset, i3_offset; + + GlobalInjector_kernel(Particles& species, + const M& global_metric, + const Domain& local_domain, + const std::map>& data, + bool use_weights) + : use_weights { use_weights } + , i1s { species.i1 } + , i2s { species.i2 } + , i3s { species.i3 } + , dx1s { species.dx1 } + , dx2s { species.dx2 } + , dx3s { species.dx3 } + , ux1s { species.ux1 } + , ux2s { species.ux2 } + , ux3s { species.ux3 } + , phis { species.phi } + , weights { species.weight } + , tags { species.tag } + , offset { species.npart() } + , global_metric { global_metric } { + const auto n_inject = data.at("x1").size(); + + x1_min = local_domain.mesh.extent(in::x1).first; + x1_max = local_domain.mesh.extent(in::x1).second; + i1_offset = local_domain.offset_ncells()[0]; + + copy_from_vector("x1", in_x1, data, n_inject); + copy_from_vector("ux1", in_ux1, data, n_inject); + copy_from_vector("ux2", in_ux2, data, n_inject); + copy_from_vector("ux3", in_ux3, data, n_inject); + if (use_weights) { + copy_from_vector("weights", in_wei, data, n_inject); + } + if constexpr (D == Dim::_2D or D == Dim::_3D) { + x2_min = local_domain.mesh.extent(in::x2).first; + x2_max = local_domain.mesh.extent(in::x2).second; + i2_offset = local_domain.offset_ncells()[1]; + copy_from_vector("x2", in_x2, data, n_inject); + } + if constexpr (D == Dim::_2D and M::CoordType != Coord::Cart) { + copy_from_vector("phi", in_phi, data, n_inject); + } + if constexpr (D == Dim::_3D) { + x3_min = local_domain.mesh.extent(in::x3).first; + x3_max = local_domain.mesh.extent(in::x3).second; + i3_offset = local_domain.offset_ncells()[2]; + copy_from_vector("x3", in_x3, data, n_inject); + } + } + + void copy_from_vector(const std::string& name, + array_t& arr, + const std::map>& data, + npart_t n_inject) { + raise::ErrorIf(data.find(name) == data.end(), name + " not found in data", HERE); + raise::ErrorIf(data.at(name).size() != n_inject, "Inconsistent data size", HERE); + arr = array_t { name, n_inject }; + auto arr_h = Kokkos::create_mirror_view(arr); + for (auto i = 0u; i < data.at(name).size(); ++i) { + arr_h(i) = data.at(name)[i]; + } + Kokkos::deep_copy(arr, arr_h); + } + + auto number_injected() const -> npart_t { + auto idx_h = Kokkos::create_mirror_view(idx); + Kokkos::deep_copy(idx_h, idx); + return idx_h(); + } + + Inline void operator()(index_t p) const { + if constexpr (D == Dim::_1D) { + if (in_x1(p) >= x1_min and in_x1(p) < x1_max) { + coord_t x_Cd { ZERO }; + vec_t u_XYZ { ZERO }; + const vec_t u_Ph { in_ux1(p), in_ux2(p), in_ux3(p) }; + + auto index { offset + Kokkos::atomic_fetch_add(&idx(), 1) }; + global_metric.template convert({ in_x1(p) }, x_Cd); + global_metric.template transform_xyz(x_Cd, u_Ph, u_XYZ); + + const auto i1 = static_cast( + static_cast(x_Cd[0]) - i1_offset); + const auto dx1 = static_cast( + x_Cd[0] - static_cast(i1 + i1_offset)); + + i1s(index) = i1; + dx1s(index) = dx1; + ux1s(index) = u_XYZ[0]; + ux2s(index) = u_XYZ[1]; + ux3s(index) = u_XYZ[2]; + tags(index) = ParticleTag::alive; + if (use_weights) { + weights(index) = weights(p); + } else { + weights(index) = ONE; + } + } + } else if constexpr (D == Dim::_2D) { + if ((in_x1(p) >= x1_min and in_x1(p) < x1_max) and + (in_x2(p) >= x2_min and in_x2(p) < x2_max)) { + coord_t x_Cd { ZERO }; + vec_t u_Cd { ZERO }; + vec_t u_Ph { in_ux1(p), in_ux2(p), in_ux3(p) }; + coord_t x_Cd_ { ZERO }; + + auto index { offset + + Kokkos::atomic_fetch_add(&idx(), static_cast(1)) }; + global_metric.template convert({ in_x1(p), in_x2(p) }, + x_Cd); + x_Cd_[0] = x_Cd[0]; + x_Cd_[1] = x_Cd[1]; + if constexpr (S == SimEngine::SRPIC and M::CoordType != Coord::Cart) { + x_Cd_[2] = in_phi(p); + } + if constexpr (S == SimEngine::SRPIC) { + global_metric.template transform_xyz(x_Cd_, u_Ph, u_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + global_metric.template transform(x_Cd, u_Ph, u_Cd); + } else { + raise::KernelError(HERE, "Unknown simulation engine"); + } + const auto i1 = static_cast( + static_cast(x_Cd[0]) - i1_offset); + const auto dx1 = static_cast( + x_Cd[0] - static_cast(i1 + i1_offset)); + const auto i2 = static_cast( + static_cast(x_Cd[1]) - i2_offset); + const auto dx2 = static_cast( + x_Cd[1] - static_cast(i2 + i2_offset)); + + i1s(index) = i1; + dx1s(index) = dx1; + i2s(index) = i2; + dx2s(index) = dx2; + ux1s(index) = u_Cd[0]; + ux2s(index) = u_Cd[1]; + ux3s(index) = u_Cd[2]; + if (M::CoordType != Coord::Cart) { + phis(index) = in_phi(p); + } + tags(index) = ParticleTag::alive; + if (use_weights) { + weights(index) = weights(p); + } else { + weights(index) = ONE; + } + } + } else { + if ((in_x1(p) >= x1_min and in_x1(p) < x1_max) and + (in_x2(p) >= x2_min and in_x2(p) < x2_max) and + (in_x3(p) >= x3_min and in_x3(p) < x3_max)) { + coord_t x_Cd { ZERO }; + vec_t u_Cd { ZERO }; + vec_t u_Ph { in_ux1(p), in_ux2(p), in_ux3(p) }; + + auto index { offset + Kokkos::atomic_fetch_add(&idx(), 1) }; + global_metric.template convert( + { in_x1(p), in_x2(p), in_x3(p) }, + x_Cd); + if constexpr (S == SimEngine::SRPIC) { + global_metric.template transform_xyz(x_Cd, u_Ph, u_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + global_metric.template transform(x_Cd, u_Ph, u_Cd); + } else { + raise::KernelError(HERE, "Unknown simulation engine"); + } + const auto i1 = static_cast( + static_cast(x_Cd[0]) - i1_offset); + const auto dx1 = static_cast( + x_Cd[0] - static_cast(i1 + i1_offset)); + const auto i2 = static_cast( + static_cast(x_Cd[1]) - i2_offset); + const auto dx2 = static_cast( + x_Cd[1] - static_cast(i2 + i2_offset)); + const auto i3 = static_cast( + static_cast(x_Cd[2]) - i3_offset); + const auto dx3 = static_cast( + x_Cd[2] - static_cast(i3 + i3_offset)); + + i1s(index) = i1; + dx1s(index) = dx1; + i2s(index) = i2; + dx2s(index) = dx2; + i3s(index) = i3; + dx3s(index) = dx3; + ux1s(index) = u_Cd[0]; + ux2s(index) = u_Cd[1]; + ux3s(index) = u_Cd[2]; + tags(index) = ParticleTag::alive; + if (use_weights) { + weights(index) = weights(p); + } else { + weights(index) = ONE; + } + } + } + } + }; // struct GlobalInjector_kernel + } // namespace kernel #endif // KERNELS_INJECTORS_HPP