diff --git a/pgens/shock/pgen.hpp b/pgens/shock/pgen.hpp index d3082bad..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 { @@ -85,15 +82,16 @@ namespace user { 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 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, Metadomain& global_domain) : arch::ProblemGenerator { p } @@ -109,13 +107,14 @@ 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_frequency { p.template get("setup.injection_frequency", 100) } - , dt { p.template get("algorithms.timestep.dt") } {} + , injection_start { static_cast( + p.template get("setup.injection_start", 0.0)) } + , injection_frequency { p.template get("setup.injection_frequency", + 100) } {} inline PGen() {} - auto MatchFields(real_t time) const -> InitFields { + auto MatchFields(real_t) const -> InitFields { return init_flds; } @@ -139,57 +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); - } - } - - // define temperatures of species - const auto temperatures = std::make_pair(temperature, - temperature_ratio * temperature); - // define drift speed of species - const auto drifts = std::make_pair( - std::vector { -drift_ux, ZERO, ZERO }, - std::vector { -drift_ux, ZERO, ZERO }); - - // inject particles - arch::InjectUniformMaxwellians(params, - domain, - ONE, - temperatures, - { 1, 2 }, - drifts, - false, - box); - } + inline void InitPrtls(Domain&) {} void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { @@ -209,77 +158,80 @@ namespace user { */ // check if the injector should be active - if (step % injection_frequency != 0) { + 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 - 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 }); + if (step == 0) { + xmin = global_xmin; } - // define box to reset fields - boundaries_t purge_box; + /* + 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) { - 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::E | Comm::B); - - /* - 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]; + // 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 }); + + 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; 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; } @@ -288,42 +240,39 @@ 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) { + ux1(p) = TWO * inj_vel / + math::max(static_cast(1e-2), + math::sqrt(ONE - SQR(inj_vel))); + 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, + 1.0, + inj_box); } }; + } // namespace user #endif diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h index f92f6def..ea6f622c 100644 --- a/src/archetypes/particle_injector.h +++ b/src/archetypes/particle_injector.h @@ -360,67 +360,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 @@ -732,6 +732,119 @@ 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, + 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"); + + 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, + // density_tolerance }; + + const auto ppc = target_tot_ndens * + params.template get("particles.ppc0") * HALF; + + auto injector_kernel = + kernel::experimental::ConstNPPCInjector_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, + 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); + + 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 69642e2f..47f58c6d 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,56 +26,6 @@ namespace kernel { using namespace ntt; - template - Inline void InjectParticle(const npart_t& p, - const array_t& i1_arr, - const array_t& i2_arr, - const array_t& i3_arr, - const array_t& dx1_arr, - const array_t& dx2_arr, - const array_t& dx3_arr, - const array_t& ux1_arr, - const array_t& ux2_arr, - const array_t& ux3_arr, - const array_t& phi_arr, - const array_t& weight_arr, - const array_t& tag_arr, - const array_t& pld_i_arr, - const tuple_t& xi_Cd, - const tuple_t& dxi_Cd, - const vec_t& v_Cd, - const real_t& weight = ONE, - const real_t& phi = ZERO, - const npart_t& domain_idx = 0u, - const npart_t& species_cntr = 0u) { - if constexpr (D == Dim::_1D or D == Dim::_2D or D == Dim::_3D) { - i1_arr(p) = xi_Cd[0]; - dx1_arr(p) = dxi_Cd[0]; - } - if constexpr (D == Dim::_2D or D == Dim::_3D) { - i2_arr(p) = xi_Cd[1]; - dx2_arr(p) = dxi_Cd[1]; - } - if constexpr (D == Dim::_3D) { - i3_arr(p) = xi_Cd[2]; - dx3_arr(p) = dxi_Cd[2]; - } - if constexpr (D == Dim::_2D and C != Coord::Cart) { - phi_arr(p) = phi; - } - ux1_arr(p) = v_Cd[0]; - ux2_arr(p) = v_Cd[1]; - ux3_arr(p) = v_Cd[2]; - tag_arr(p) = ParticleTag::alive; - weight_arr(p) = weight; - if constexpr (T) { - pld_i_arr(p, pldi::spcCtr) = species_cntr; -#if defined(MPI_ENABLED) - pld_i_arr(p, pldi::domIdx) = domain_idx; -#endif - } - } - template struct UniformInjector_kernel { static_assert(ED1::is_energy_dist, "ED1 must be an energy distribution class"); @@ -298,240 +249,10 @@ namespace kernel { } }; // struct UniformInjector_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 - - 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"); @@ -556,7 +277,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; @@ -569,7 +291,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) @@ -603,7 +326,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 } {} @@ -620,11 +344,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); @@ -635,12 +360,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]; @@ -676,15 +401,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; @@ -697,7 +426,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) { @@ -706,7 +435,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) { @@ -727,7 +456,6 @@ namespace kernel { weights_2(index + offset2) = wei; } } - random_pool.free_state(rand_gen); } else { @@ -743,11 +471,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); @@ -770,7 +499,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) { @@ -779,7 +508,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) { @@ -809,6 +538,526 @@ 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 { + 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 + + 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