diff --git a/TODO.md b/TODO.md new file mode 100644 index 00000000..95278bf2 --- /dev/null +++ b/TODO.md @@ -0,0 +1,3 @@ +- [ ] optional weights in communications +- [ ] make weight buffers for computing moments optional +- [ ] make writing (outputting) of weights optional diff --git a/input.example.toml b/input.example.toml index 6444a16f..d0db7240 100644 --- a/input.example.toml +++ b/input.example.toml @@ -264,7 +264,7 @@ ppc0 = "" # Toggle for using particle weights # @type: bool - # @default: false + # @default: false [Cartesian]; true [otherwise] use_weights = "" # Timesteps between particle re-sorting (removing dead particles) # @type: uint diff --git a/src/framework/containers/particles.cpp b/src/framework/containers/particles.cpp index d2db9c49..39920942 100644 --- a/src/framework/containers/particles.cpp +++ b/src/framework/containers/particles.cpp @@ -22,10 +22,11 @@ namespace ntt { float ch, npart_t maxnpart, const PrtlPusher& pusher, + bool use_weights, bool use_gca, const Cooling& cooling, unsigned short npld) - : ParticleSpecies(index, label, m, ch, maxnpart, pusher, use_gca, cooling, npld) { + : ParticleSpecies(index, label, m, ch, maxnpart, pusher, use_weights, use_gca, cooling, npld) { if constexpr (D == Dim::_1D or D == Dim::_2D or D == Dim::_3D) { i1 = array_t { label + "_i1", maxnpart }; @@ -52,7 +53,9 @@ namespace ntt { ux2 = array_t { label + "_ux2", maxnpart }; ux3 = array_t { label + "_ux3", maxnpart }; - weight = array_t { label + "_w", maxnpart }; + if (this->use_weights()) { + weight = array_t { label + "_w", maxnpart }; + } tag = array_t { label + "_tag", maxnpart }; @@ -200,7 +203,10 @@ namespace ntt { RemoveDeadInArray(ux1, indices_alive); RemoveDeadInArray(ux2, indices_alive); RemoveDeadInArray(ux3, indices_alive); - RemoveDeadInArray(weight, indices_alive); + + if (this->use_weights()) { + RemoveDeadInArray(weight, indices_alive); + } if constexpr (D == Dim::_2D && C != Coord::Cart) { RemoveDeadInArray(phi, indices_alive); diff --git a/src/framework/containers/particles.h b/src/framework/containers/particles.h index 8ff74be3..4a5e7df6 100644 --- a/src/framework/containers/particles.h +++ b/src/framework/containers/particles.h @@ -76,6 +76,7 @@ namespace ntt { * @param ch The charge of the species * @param maxnpart The maximum number of allocated particles for the species * @param pusher The pusher assigned for the species + * @param use_weights Use weights for the species * @param use_gca Use hybrid GCA pusher for the species * @param cooling The cooling mechanism assigned for the species * @param npld The number of payloads for the species @@ -86,6 +87,7 @@ namespace ntt { float ch, npart_t maxnpart, const PrtlPusher& pusher, + bool use_weights, bool use_gca, const Cooling& cooling, unsigned short npld = 0); @@ -102,6 +104,7 @@ namespace ntt { spec.charge(), spec.maxnpart(), spec.pusher(), + spec.use_weights(), spec.use_gca(), spec.cooling(), spec.npld()) {} diff --git a/src/framework/containers/species.h b/src/framework/containers/species.h index ada0282e..5f8586ab 100644 --- a/src/framework/containers/species.h +++ b/src/framework/containers/species.h @@ -33,6 +33,8 @@ namespace ntt { // Pusher assigned for the species const PrtlPusher m_pusher; + // Use particle weights + const bool m_use_weights; // Use byrid gca pusher for the species const bool m_use_gca; @@ -50,6 +52,7 @@ namespace ntt { , m_charge { 0.0 } , m_maxnpart { 0 } , m_pusher { PrtlPusher::INVALID } + , m_use_weights { false } , m_use_gca { false } , m_cooling { Cooling::INVALID } , m_npld { 0 } {} @@ -70,6 +73,7 @@ namespace ntt { float ch, npart_t maxnpart, const PrtlPusher& pusher, + bool use_weights, bool use_gca, const Cooling& cooling, unsigned short npld = 0) @@ -79,6 +83,7 @@ namespace ntt { , m_charge { ch } , m_maxnpart { maxnpart } , m_pusher { pusher } + , m_use_weights { use_weights } , m_use_gca { use_gca } , m_cooling { cooling } , m_npld { npld } {} @@ -120,6 +125,11 @@ namespace ntt { return m_pusher; } + [[nodiscard]] + auto use_weights() const -> bool { + return m_use_weights; + } + [[nodiscard]] auto use_gca() const -> bool { return m_use_gca; diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp index 079ad615..53e5ac9c 100644 --- a/src/framework/parameters.cpp +++ b/src/framework/parameters.cpp @@ -173,8 +173,11 @@ namespace ntt { const auto ppc0 = toml::find(toml_data, "particles", "ppc0"); set("particles.ppc0", ppc0); raise::ErrorIf(ppc0 <= 0.0, "ppc0 must be positive", HERE); - set("particles.use_weights", - toml::find_or(toml_data, "particles", "use_weights", false)); + const auto use_weights = toml::find_or(toml_data, + "particles", + "use_weights", + coord_enum != Coord::Cart); + set("particles.use_weights", use_weights); /* [particles.species] -------------------------------------------------- */ std::vector species; @@ -241,6 +244,7 @@ namespace ntt { charge, maxnpart, pusher_enum, + use_weights, use_gca, cooling_enum, npayloads)); diff --git a/src/kernels/currents_deposit.hpp b/src/kernels/currents_deposit.hpp index bd84554f..3f470508 100644 --- a/src/kernels/currents_deposit.hpp +++ b/src/kernels/currents_deposit.hpp @@ -43,6 +43,8 @@ namespace kernel { const M metric; const real_t charge, inv_dt; + const bool use_weights; + public: /** * @brief explicit constructor. @@ -90,7 +92,8 @@ namespace kernel { , tag { tag } , metric { metric } , charge { charge } - , inv_dt { ONE / dt } {} + , inv_dt { ONE / dt } + , use_weights { weight.extent(0) > 0 } {} /** * @brief Iteration of the loop over particles. @@ -155,7 +158,7 @@ namespace kernel { vp[2] *= inv_energy; } - const real_t coeff { weight(p) * charge }; + const real_t coeff { (use_weights ? weight(p) : ONE) * charge }; const auto dxp_r_1 { static_cast(i1(p) == i1_prev(p)) * (dx1(p) + dx1_prev(p)) * static_cast(INV_2) }; diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index c321d18f..56c637d0 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -53,6 +53,8 @@ namespace kernel { const real_t inv_V0; random_number_pool_t random_pool; + const bool use_weights; + UniformInjector_kernel(spidx_t spidx1, spidx_t spidx2, Particles& species1, @@ -98,7 +100,13 @@ namespace kernel { , xi_max { xi_max } , energy_dist { energy_dist } , inv_V0 { inv_V0 } - , random_pool { random_pool } {} + , random_pool { random_pool } + , use_weights { (weights_1.extent(0) > 0) and (weights_2.extent(0) > 0) } { + raise::ErrorIf((weights_1.extent(0) > 0) != (weights_2.extent(0) > 0), + "Weights for species 1 and 2 must be either both defined " + "or both undefined", + HERE); + } Inline void operator()(index_t p) const { coord_t x_Cd { ZERO }; @@ -173,12 +181,16 @@ namespace kernel { 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; + if (use_weights) { + 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; + if (use_weights) { + 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 @@ -217,6 +229,8 @@ namespace kernel { const real_t inv_V0; random_number_pool_t random_pool; + const bool use_weights; + UniformInjector_kernel(spidx_t spidx1, spidx_t spidx2, Particles& species1, @@ -264,7 +278,14 @@ namespace kernel { , energy_dist_1 { energy_dist_1 } , energy_dist_2 { energy_dist_2 } , inv_V0 { inv_V0 } - , random_pool { random_pool } {} + , random_pool { random_pool } + , use_weights { (weights_1.extent(0) > 0) and (weights_2.extent(0) > 0) } { + raise::ErrorIf( + (weights_1.extent(0) > 0) != (weights_2.extent(0) > 0), + "Weights for species 1 and 2 must be either both defined " + "or both undefined", + HERE); + } // UniformInjector_kernel Inline void operator()(index_t p) const { coord_t x_Cd { ZERO }; @@ -341,12 +362,16 @@ namespace kernel { 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; + if (use_weights) { + 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; + if (use_weights) { + 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 @@ -384,13 +409,13 @@ namespace kernel { real_t x1_min, x1_max, x2_min, x2_max, x3_min, x3_max; ncells_t i1_offset, i2_offset, i3_offset; + const bool use_weights; + 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 } + const std::map>& data) + : i1s { species.i1 } , i2s { species.i2 } , i3s { species.i3 } , dx1s { species.dx1 } @@ -403,7 +428,8 @@ namespace kernel { , weights { species.weight } , tags { species.tag } , offset { species.npart() } - , global_metric { global_metric } { + , global_metric { global_metric } + , use_weights { weights.extent(0) > 0 } { const auto n_inject = data.at("x1").size(); x1_min = local_domain.mesh.extent(in::x1).first; @@ -477,9 +503,7 @@ namespace kernel { ux3s(index) = u_XYZ[2]; tags(index) = ParticleTag::alive; if (use_weights) { - weights(index) = weights(p); - } else { - weights(index) = ONE; + weights(index) = in_wei(p); } } } else if constexpr (D == Dim::_2D) { @@ -527,9 +551,7 @@ namespace kernel { } tags(index) = ParticleTag::alive; if (use_weights) { - weights(index) = weights(p); - } else { - weights(index) = ONE; + weights(index) = in_wei(p); } } } else { @@ -575,9 +597,7 @@ namespace kernel { ux3s(index) = u_Cd[2]; tags(index) = ParticleTag::alive; if (use_weights) { - weights(index) = weights(p); - } else { - weights(index) = ONE; + weights(index) = in_wei(p); } } } @@ -616,6 +636,8 @@ namespace kernel { const real_t inv_V0; random_number_pool_t random_pool; + const bool use_weights; + NonUniformInjector_kernel(real_t ppc0, spidx_t spidx1, spidx_t spidx2, @@ -661,7 +683,16 @@ namespace kernel { , energy_dist { energy_dist } , spatial_dist { spatial_dist } , inv_V0 { inv_V0 } - , random_pool { random_pool } {} + , random_pool { random_pool } + , use_weights { (weights_1.extent(0) > 0) and (weights_2.extent(0) > 0) } { + { + raise::ErrorIf( + (weights_1.extent(0) > 0) != (weights_2.extent(0) > 0), + "Weights for species 1 and 2 must be either both defined " + "or both undefined", + HERE); + } + } auto number_injected() const -> npart_t { auto idx_h = Kokkos::create_mirror_view(idx); @@ -704,12 +735,16 @@ namespace kernel { 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; + if (use_weights) { + 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; + if (use_weights) { + 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); @@ -774,12 +809,17 @@ namespace kernel { 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; + if (use_weights) { + 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; + if (use_weights) { + 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); @@ -847,14 +887,18 @@ namespace kernel { 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; + if (use_weights) { + 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; + if (use_weights) { + 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); diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index 9c66aed0..085f035b 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -56,7 +56,6 @@ namespace kernel { const array_t tag; const float mass; const float charge; - const bool use_weights; const M metric; const int ni2; const unsigned short window; @@ -65,6 +64,8 @@ namespace kernel { const real_t smooth; bool is_axis_i2min { false }, is_axis_i2max { false }; + const bool use_weights; + public: ParticleMoments_kernel(const std::vector& components, const scatter_ndfield_t& scatter_buff, @@ -83,7 +84,6 @@ namespace kernel { const array_t& tag, float mass, float charge, - bool use_weights, const M& metric, const boundaries_t& boundaries, ncells_t ni2, @@ -109,13 +109,13 @@ namespace kernel { , tag { tag } , mass { mass } , charge { charge } - , use_weights { use_weights } , metric { metric } , ni2 { static_cast(ni2) } , window { window } , contrib { get_contrib(mass, charge) } , smooth { inv_n0 / (real_t)(math::pow(TWO * (real_t)window + ONE, - static_cast(D))) } { + static_cast(D))) } + , use_weights { weight.extent(0) > 0 } { raise::ErrorIf(buff_idx >= N, "Invalid buffer index", HERE); raise::ErrorIf(window > N_GHOSTS, "Window size too large", HERE); raise::ErrorIf(((F == FldsID::Rho) || (F == FldsID::Charge)) && (mass == ZERO), diff --git a/src/kernels/prtls_to_phys.hpp b/src/kernels/prtls_to_phys.hpp index 4dd7d88b..4e694ddf 100644 --- a/src/kernels/prtls_to_phys.hpp +++ b/src/kernels/prtls_to_phys.hpp @@ -46,6 +46,8 @@ namespace kernel { const array_t weight; const M metric; + const bool use_weights; + public: PrtlToPhys_kernel(npart_t stride, array_t& buff_x1, @@ -86,7 +88,8 @@ namespace kernel { , ux3 { ux3 } , phi { phi } , weight { weight } - , metric { metric } { + , metric { metric } + , use_weights { weight.extent(0) > 0 } { if constexpr ((D == Dim::_1D) || (D == Dim::_2D) || (D == Dim::_3D)) { raise::ErrorIf(buff_x1.extent(0) == 0, "Invalid buffer size", HERE); } @@ -105,7 +108,9 @@ namespace kernel { Inline void operator()(index_t p) const { bufferX(p); bufferU(p); - buff_wei(p) = weight(p * stride); + if (use_weights) { + buff_wei(p) = weight(p * stride); + } } Inline void bufferX(index_t& p) const { diff --git a/src/kernels/reduced_stats.hpp b/src/kernels/reduced_stats.hpp index 68e6aa97..acef4d88 100644 --- a/src/kernels/reduced_stats.hpp +++ b/src/kernels/reduced_stats.hpp @@ -415,11 +415,12 @@ namespace kernel { const array_t tag; const float mass; const float charge; - const bool use_weights; const M metric; const real_t contrib; + const bool use_weights; + public: ReducedParticleMoments_kernel(const std::vector& components, const array_t& i1, @@ -436,8 +437,7 @@ namespace kernel { const array_t& tag, float mass, float charge, - bool use_weights, - const M& metric) + const M& metric) : c1 { (components.size() > 0) ? components[0] : static_cast(0) } , c2 { (components.size() == 2) ? components[1] @@ -456,9 +456,9 @@ namespace kernel { , tag { tag } , mass { mass } , charge { charge } - , use_weights { use_weights } , metric { metric } - , contrib { get_contrib

(mass, charge) } { + , contrib { get_contrib

(mass, charge) } + , use_weights { weight.extent(0) > 0 } { raise::ErrorIf(((P == StatsID::Rho) || (P == StatsID::Charge)) && (mass == ZERO), "Rho & Charge for massless particles not defined",