From acf6cac4ef0e60e0e56045d160cd55fa465642ab Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Thu, 6 Aug 2026 08:56:34 +0200 Subject: [PATCH 1/6] =?UTF-8?q?*refactor=20`ParticleAttrib::scatter()`=20s?= =?UTF-8?q?o=20the=20=E2=80=9Chashed=20scatter=E2=80=9D=20vs=20=E2=80=9Cpl?= =?UTF-8?q?ain=20scatter=E2=80=9D=20decision=20is=20made=20on=20the=20host?= =?UTF-8?q?=20before=20launching=20the=20Kokkos=20kernel.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit **Purpose** The old code did this inside every particle iteration: ```cpp size_t mapped_idx = useHashView ? hash_array(idx) : idx; ``` `useHashView` is a runtime boolean captured into the kernel. That means every particle pays for a branch, and the kernel body always mentions `hash_array(idx)` even when no hash view is used. The new code splits this into two template instantiations: ```cpp scatterImpl(...) scatterImpl(...) ``` Inside the kernel it becomes: ```cpp size_t mapped_idx = idx; if constexpr (UseHashView) { mapped_idx = hash_array(idx); } ``` So for the non-hashed path, the hash-array access is compiled out entirely. **Secondary Purpose** The change also moves the hash extent check out of the kernel implementation wrapper and normalizes the type: ```cpp const auto hashExtent = static_cast(hash_array.extent(0)); ``` This avoids the Kokkos 5.2 OpenMP/GCC signed/unsigned warning from comparing `iteration_policy.end()` with `hash_array.extent(0)` directly. **Behavioral Impact** The public scatter API does not change. Expected behavior stays the same: - no hash array: `mapped_idx = idx` - hash array present: `mapped_idx = hash_array(idx)` - hash array too small for the requested iteration policy: abort with the existing diagnostic This is mainly a compile-time dispatch/performance/cleanup change: - removes one per-particle runtime branch - avoids compiling hash access into the non-hash kernel - fixes the signed/unsigned comparison warning - keeps the existing bounds check and scatter semantics ** ToDo - Test on GPUs, on CPU (MAC) no improvement - The current cast is probably fine for realistic particle counts, but a more defensive version could use a named policy index type and possibly check representability before casting. For IPPL’s current usage, this is not a practicaxl concern. --- src/Particle/ParticleAttrib.h | 4 ++++ src/Particle/ParticleAttrib.hpp | 34 ++++++++++++++++++++++++--------- 2 files changed, 29 insertions(+), 9 deletions(-) diff --git a/src/Particle/ParticleAttrib.h b/src/Particle/ParticleAttrib.h index 80b9dad20a..d2725bea2c 100644 --- a/src/Particle/ParticleAttrib.h +++ b/src/Particle/ParticleAttrib.h @@ -272,6 +272,10 @@ namespace ippl { void internalCopy(const hash_type& indices) override; private: + template + void scatterImpl(Field& f, const ParticleAttrib, Properties...>& pp, + policy_type iteration_policy, hash_type hash_array) const; + view_type dview_m{"ParticleAttrib::dview", 0}; view_type buf_m{"ParticleAttrib::buf", 0}; }; diff --git a/src/Particle/ParticleAttrib.hpp b/src/Particle/ParticleAttrib.hpp index 15cb155d9c..944ae9950f 100644 --- a/src/Particle/ParticleAttrib.hpp +++ b/src/Particle/ParticleAttrib.hpp @@ -131,6 +131,27 @@ namespace ippl { template template void ParticleAttrib::scatter( + Field& f, const ParticleAttrib, Properties...>& pp, + policy_type iteration_policy, hash_type hash_array) const { + const auto hashExtent = static_cast(hash_array.extent(0)); + const bool useHashView = hashExtent > 0; + if (useHashView && (iteration_policy.end() > hashExtent)) { + Inform m("scatter"); + m << "Hash array was passed to scatter, but size does not match iteration policy." + << endl; + ippl::Comm->abort(); + } + + if (useHashView) { + scatterImpl(f, pp, iteration_policy, hash_array); + } else { + scatterImpl(f, pp, iteration_policy, hash_array); + } + } + + template + template + void ParticleAttrib::scatterImpl( Field& f, const ParticleAttrib, Properties...>& pp, policy_type iteration_policy, hash_type hash_array) const { constexpr unsigned Dim = Field::dim; @@ -155,20 +176,15 @@ namespace ippl { const NDIndex& lDom = layout.getLocalNDIndex(); const int nghost = f.getNghost(); - // using policy_type = Kokkos::RangePolicy; - const bool useHashView = hash_array.extent(0) > 0; - if (useHashView && (iteration_policy.end() > hash_array.extent(0))) { - Inform m("scatter"); - m << "Hash array was passed to scatter, but size does not match iteration policy." - << endl; - ippl::Comm->abort(); - } auto dview = dview_m; auto ppview = pp.getView(); Kokkos::parallel_for( "ParticleAttrib::scatter", iteration_policy, KOKKOS_LAMBDA(const size_t idx) { // map index to possible hash_map - size_t mapped_idx = useHashView ? hash_array(idx) : idx; + size_t mapped_idx = idx; + if constexpr (UseHashView) { + mapped_idx = hash_array(idx); + } vector_type l = (ppview(mapped_idx) - origin) * invdx + 0.5; Vector index = l; From 2a4669cee4f62fc503fce59561cf60fdb08619c2 Mon Sep 17 00:00:00 2001 From: Andreas Adelmann Date: Thu, 6 Aug 2026 20:50:29 +0200 Subject: [PATCH 2/6] Make hashed ParticleAttrib scatter CUDA-safe Port the CUDA-safe hashed scatter implementation from fixissue/#415 onto the hashed-scatter branch. Move the kernel-bearing scatter implementation out of the private ParticleAttrib member scope and into a namespace-scope detail helper. CUDA extended host-device lambdas cannot be enclosed by private or protected class member functions, so keeping the KOKKOS_LAMBDA inside the former private scatterImpl caused GH200/NVCC builds to fail. Keep ParticleAttrib::scatter as the public host-side dispatcher. It validates the optional hash view and selects the hashed or non-hashed implementation via: particleAttribScatterImpl particleAttribScatterImpl Route mapped-index selection through a small KOKKOS_INLINE_FUNCTION helper so the hash lookup is compiled only for the hashed instantiation and is not first captured inside an if constexpr context. Also allow CIC scatter to accept a value type that differs from the field value type. This enables mixed-value scatter such as: ParticleAttrib -> Field by casting the scattered value once to the field view value type before accumulating into the grid. Add regression coverage for: - plain ParticleAttrib -> Field scatter - hashed ParticleAttrib -> Field scatter - dimensions 1, 2, and 3 through GatherScatterTest Tested locally with: cmake --build build-fixissue-415-merge --target GatherScatterTest -j 8 ctest --test-dir build-fixissue-415-merge --output-on-failure -R '^GatherScatterTest$' mpirun -np 2 build-fixissue-415-merge/unit_tests/Particle/GatherScatterTest --gtest_filter='*MixedValueType*' mpirun -np 4 build-fixissue-415-merge/unit_tests/Particle/GatherScatterTest --gtest_filter='*MixedValueType*' All tests passed. --- src/Interpolation/CIC.h | 5 +- src/Interpolation/CIC.hpp | 11 +- src/Particle/ParticleAttrib.h | 16 +- src/Particle/ParticleAttrib.hpp | 133 +++++------ test/particle/TestScatter.cpp | 42 +++- unit_tests/Particle/GatherScatterTest.cpp | 257 ++++++++++++++++------ 6 files changed, 314 insertions(+), 150 deletions(-) diff --git a/src/Interpolation/CIC.h b/src/Interpolation/CIC.h index 5adac510d3..80e7b3a05f 100644 --- a/src/Interpolation/CIC.h +++ b/src/Interpolation/CIC.h @@ -41,6 +41,7 @@ namespace ippl { * @tparam View the field view type * @tparam T the field data type * @tparam IndexType the index type for accessing the field (default size_t) + * @tparam Val the type of the value to interpolate to the grid * @param view the field view on which to scatter * @param wlo lower weights for interpolation * @param whi upper weights for interpolation @@ -75,11 +76,11 @@ namespace ippl { * @param val the value to interpolate */ template + typename IndexType = size_t, typename Val = T> KOKKOS_INLINE_FUNCTION constexpr void scatterToField( const std::index_sequence&, const View& view, const Vector& wlo, const Vector& whi, - const Vector& args, T val = 1); + const Vector& args, Val val = T(1)); /*! * Gathers from a field at a single point diff --git a/src/Interpolation/CIC.hpp b/src/Interpolation/CIC.hpp index 0f8d1732c6..dd6e75e0c3 100644 --- a/src/Interpolation/CIC.hpp +++ b/src/Interpolation/CIC.hpp @@ -1,5 +1,6 @@ #include +#include namespace ippl { namespace detail { @@ -33,14 +34,18 @@ namespace ippl { val * (interpolationWeight(wlo, whi) * ...)); } - template + template KOKKOS_INLINE_FUNCTION constexpr void scatterToField( const std::index_sequence&, const View& view, const Vector& wlo, const Vector& whi, - const Vector& args, T val) { + const Vector& args, Val val) { // The number of indices is equal to the view rank + using out_type = + std::remove_cv_t>; + const out_type field_value = static_cast(val); (scatterToPoint(std::make_index_sequence{}, view, wlo, whi, - args, val), + args, field_value), ...); } diff --git a/src/Particle/ParticleAttrib.h b/src/Particle/ParticleAttrib.h index d2725bea2c..56e1b85353 100644 --- a/src/Particle/ParticleAttrib.h +++ b/src/Particle/ParticleAttrib.h @@ -137,14 +137,14 @@ namespace ippl { * particle range [0, size()). */ view_type getView() { - return Kokkos::subview(dview_m, - Kokkos::make_pair(size_type(0), - static_cast(*(this->localNum_mp)))); + return Kokkos::subview( + dview_m, + Kokkos::make_pair(size_type(0), static_cast(*(this->localNum_mp)))); } const view_type getView() const { - return Kokkos::subview(dview_m, - Kokkos::make_pair(size_type(0), - static_cast(*(this->localNum_mp)))); + return Kokkos::subview( + dview_m, + Kokkos::make_pair(size_type(0), static_cast(*(this->localNum_mp)))); } host_mirror_type getHostMirror() const { return Kokkos::create_mirror(getView()); } @@ -272,10 +272,6 @@ namespace ippl { void internalCopy(const hash_type& indices) override; private: - template - void scatterImpl(Field& f, const ParticleAttrib, Properties...>& pp, - policy_type iteration_policy, hash_type hash_array) const; - view_type dview_m{"ParticleAttrib::dview", 0}; view_type buf_m{"ParticleAttrib::buf", 0}; }; diff --git a/src/Particle/ParticleAttrib.hpp b/src/Particle/ParticleAttrib.hpp index 944ae9950f..ca477f21d0 100644 --- a/src/Particle/ParticleAttrib.hpp +++ b/src/Particle/ParticleAttrib.hpp @@ -16,7 +16,6 @@ #include "Ippl.h" #include - #include #include "Communicate/DataTypes.h" @@ -31,6 +30,72 @@ #include "Particle/SortBuffer.h" namespace ippl { + namespace detail { + template + KOKKOS_INLINE_FUNCTION size_t scatterMappedIndex(const size_t idx, + const HashView& hashView) { + if constexpr (UseHashView) { + return static_cast(hashView(idx)); + } else { + (void)hashView; + return idx; + } + } + + template + void particleAttribScatterImpl(Field& f, const ValuesView& dview, const PositionAttrib& pp, + policy_type iteration_policy, HashView hash_array) { + constexpr unsigned Dim = Field::dim; + using PositionType = typename Field::Mesh_t::value_type; + + static IpplTimings::TimerRef scatterTimer = IpplTimings::getTimer("scatter"); + IpplTimings::startTimer(scatterTimer); + using view_type = typename Field::view_type; + view_type view = f.getView(); + + using mesh_type = typename Field::Mesh_t; + const mesh_type& mesh = f.get_mesh(); + + using vector_type = typename mesh_type::vector_type; + using value_type = typename ValuesView::non_const_value_type; + + const vector_type& dx = mesh.getMeshSpacing(); + const vector_type& origin = mesh.getOrigin(); + const vector_type invdx = 1.0 / dx; + + const FieldLayout& layout = f.getLayout(); + const NDIndex& lDom = layout.getLocalNDIndex(); + const int nghost = f.getNghost(); + + auto ppview = pp.getView(); + auto hashView = hash_array; + Kokkos::parallel_for( + "ParticleAttrib::scatter", iteration_policy, KOKKOS_LAMBDA(const size_t idx) { + // map index to possible hash_map + const size_t mapped_idx = + detail::scatterMappedIndex(idx, hashView); + + vector_type l = (ppview(mapped_idx) - origin) * invdx + 0.5; + Vector index = l; + Vector whi = l - index; + Vector wlo = 1.0 - whi; + + Vector args = index - lDom.first() + nghost; + + const value_type& val = dview(mapped_idx); + detail::scatterToField(std::make_index_sequence<1 << Field::dim>{}, view, wlo, + whi, args, val); + }); + IpplTimings::stopTimer(scatterTimer); + + static IpplTimings::TimerRef accumulateHaloTimer = + IpplTimings::getTimer("accumulateHalo"); + IpplTimings::startTimer(accumulateHaloTimer); + f.accumulateHalo(); + IpplTimings::stopTimer(accumulateHaloTimer); + } + } // namespace detail template void ParticleAttrib::create(size_type n, bool non_destructive) { @@ -143,68 +208,12 @@ namespace ippl { } if (useHashView) { - scatterImpl(f, pp, iteration_policy, hash_array); + detail::particleAttribScatterImpl(f, dview_m, pp, iteration_policy, hash_array); } else { - scatterImpl(f, pp, iteration_policy, hash_array); + detail::particleAttribScatterImpl(f, dview_m, pp, iteration_policy, hash_array); } } - template - template - void ParticleAttrib::scatterImpl( - Field& f, const ParticleAttrib, Properties...>& pp, - policy_type iteration_policy, hash_type hash_array) const { - constexpr unsigned Dim = Field::dim; - using PositionType = typename Field::Mesh_t::value_type; - - static IpplTimings::TimerRef scatterTimer = IpplTimings::getTimer("scatter"); - IpplTimings::startTimer(scatterTimer); - using view_type = typename Field::view_type; - view_type view = f.getView(); - - using mesh_type = typename Field::Mesh_t; - const mesh_type& mesh = f.get_mesh(); - - using vector_type = typename mesh_type::vector_type; - using value_type = typename ParticleAttrib::value_type; - - const vector_type& dx = mesh.getMeshSpacing(); - const vector_type& origin = mesh.getOrigin(); - const vector_type invdx = 1.0 / dx; - - const FieldLayout& layout = f.getLayout(); - const NDIndex& lDom = layout.getLocalNDIndex(); - const int nghost = f.getNghost(); - - auto dview = dview_m; - auto ppview = pp.getView(); - Kokkos::parallel_for( - "ParticleAttrib::scatter", iteration_policy, KOKKOS_LAMBDA(const size_t idx) { - // map index to possible hash_map - size_t mapped_idx = idx; - if constexpr (UseHashView) { - mapped_idx = hash_array(idx); - } - - vector_type l = (ppview(mapped_idx) - origin) * invdx + 0.5; - Vector index = l; - Vector whi = l - index; - Vector wlo = 1.0 - whi; - - Vector args = index - lDom.first() + nghost; - - const value_type& val = dview(mapped_idx); - detail::scatterToField(std::make_index_sequence<1 << Field::dim>{}, view, wlo, whi, - args, val); - }); - IpplTimings::stopTimer(scatterTimer); - - static IpplTimings::TimerRef accumulateHaloTimer = IpplTimings::getTimer("accumulateHalo"); - IpplTimings::startTimer(accumulateHaloTimer); - f.accumulateHalo(); - IpplTimings::stopTimer(accumulateHaloTimer); - } - template template void ParticleAttrib::gather( @@ -426,8 +435,7 @@ namespace ippl { Field& f, Field& Sk, const ParticleAttrib, Properties...>& pp, FFT>* nufft, ParticleAttrib& q) { - static IpplTimings::TimerRef gatherPIFNUFFTTimer = - IpplTimings::getTimer("GatherPIFNUFFT"); + static IpplTimings::TimerRef gatherPIFNUFFTTimer = IpplTimings::getTimer("GatherPIFNUFFT"); IpplTimings::startTimer(gatherPIFNUFFTTimer); typename Field::uniform_type tempField; @@ -468,8 +476,7 @@ namespace ippl { "Gather NUFFT", mdrange_type( {nghost, nghost, nghost}, - {fview.extent(0) - nghost, fview.extent(1) - nghost, - fview.extent(2) - nghost}), + {fview.extent(0) - nghost, fview.extent(1) - nghost, fview.extent(2) - nghost}), KOKKOS_LAMBDA(const int i, const int j, const int k) { Vector iVec = {i, j, k}; for (unsigned d = 0; d < Dim; ++d) { diff --git a/test/particle/TestScatter.cpp b/test/particle/TestScatter.cpp index bfa11e9daf..8ca2403a13 100644 --- a/test/particle/TestScatter.cpp +++ b/test/particle/TestScatter.cpp @@ -6,13 +6,17 @@ template struct Bunch : public ippl::ParticleBase { Bunch(PLayout& playout) : ippl::ParticleBase(playout) { - this->addAttribute(Q); + this->addAttribute(Q1); + this->addAttribute(Q2); } ~Bunch() {} - typedef ippl::ParticleAttrib charge_container_type; - charge_container_type Q; + typedef ippl::ParticleAttrib charge_container_typeF; + charge_container_typeF Q1; + + typedef ippl::ParticleAttrib charge_container_typeD; + charge_container_typeD Q2; }; int main(int argc, char* argv[]) { @@ -85,26 +89,46 @@ int main(int argc, char* argv[]) { std::cout << "Sum coord: " << global_sum_coord << std::endl; } - bunch.Q = 1.0; + bunch.Q1 = 1.0; + + bunch.update(); + + field = 0.0; + + scatter(bunch.Q1, field, bunch.R); + + // Check charge conservation + try { + double Total_charge_field = field.sum(); + + std::cout << "Float:: Total charge in the field:" << Total_charge_field << std::endl; + std::cout << "Float:: Total charge of the particles:" << bunch.Q1.sum() << std::endl; + std::cout << "Float:: Error:" << std::fabs(bunch.Q1.sum() - Total_charge_field) + << std::endl; + } catch (const std::exception& e) { + std::cout << e.what() << std::endl; + } + + bunch.Q2 = 1.0; bunch.update(); field = 0.0; - scatter(bunch.Q, field, bunch.R); + scatter(bunch.Q2, field, bunch.R); // Check charge conservation try { double Total_charge_field = field.sum(); - std::cout << "Total charge in the field:" << Total_charge_field << std::endl; - std::cout << "Total charge of the particles:" << bunch.Q.sum() << std::endl; - std::cout << "Error:" << std::fabs(bunch.Q.sum() - Total_charge_field) << std::endl; + std::cout << "Double:: Total charge in the field:" << Total_charge_field << std::endl; + std::cout << "Double:: Total charge of the particles:" << bunch.Q2.sum() << std::endl; + std::cout << "Double:: Error:" << std::fabs(bunch.Q2.sum() - Total_charge_field) + << std::endl; } catch (const std::exception& e) { std::cout << e.what() << std::endl; } } ippl::finalize(); - return 0; } diff --git a/unit_tests/Particle/GatherScatterTest.cpp b/unit_tests/Particle/GatherScatterTest.cpp index 5a05f16189..387ce03ec0 100644 --- a/unit_tests/Particle/GatherScatterTest.cpp +++ b/unit_tests/Particle/GatherScatterTest.cpp @@ -9,41 +9,60 @@ // #include "Ippl.h" -#include "TestUtils.h" -#include "gtest/gtest.h" -#include -#include -#include #include +#include #include +#include +#include +#include +#include + +#include "TestUtils.h" +#include "gtest/gtest.h" // A helper needed to reduce over a hash_type. // This is needed, since Kokkos kernels apparently // cannot be called inside a TYPED_TEST on device. +template struct ComputeTotalChargeLambda { - Kokkos::View viewQ; - Kokkos::View hash; + ViewType viewQ; + HashType hash; - ComputeTotalChargeLambda(Kokkos::View viewQ_, Kokkos::View hash_) - : viewQ(viewQ_), hash(hash_) {} + ComputeTotalChargeLambda(ViewType viewQ_, HashType hash_) + : viewQ(viewQ_) + , hash(hash_) {} KOKKOS_INLINE_FUNCTION void operator()(const size_t i, double& val) const { - val += viewQ(hash(i)); + val += static_cast(viewQ(hash(i))); } }; -// A simple bunch_type holding a charge attribute +template +double scatterConservationTolerance(const double reference) { + using tolerance_type = std::conditional_t<(std::numeric_limits::digits + < std::numeric_limits::digits), + ScatterType, FieldType>; + return 100.0 * static_cast(std::numeric_limits::epsilon()) + * std::max(1.0, std::abs(reference)); +} + +// A simple bunch_type holding a charge attribute template struct Bunch : public ippl::ParticleBase { Bunch(PLayout& playout) : ippl::ParticleBase(playout) { this->addAttribute(Q); + this->addAttribute(QFloat); } ~Bunch() = default; - typedef ippl::ParticleAttrib charge_container_type; + typedef ippl::ParticleAttrib + charge_container_type; + typedef ippl::ParticleAttrib + mixed_charge_container_type; charge_container_type Q; + mixed_charge_container_type QFloat; }; template @@ -52,30 +71,30 @@ class GatherScatterTest; template class GatherScatterTest>> : public ::testing::Test { public: - using scalar_type = T; - using exec_space = ExecSpace; + using scalar_type = T; + using exec_space = ExecSpace; static const unsigned dim = Dim; - using flayout_type = ippl::FieldLayout; - using mesh_type = ippl::UniformCartesian; - using playout_type = ippl::ParticleSpatialLayout; - using bunch_type = Bunch; + using flayout_type = ippl::FieldLayout; + using mesh_type = ippl::UniformCartesian; + using playout_type = ippl::ParticleSpatialLayout; + using bunch_type = Bunch; // Domain parameters: use a high resolution grid so that cells are small. std::array nPoints; std::array domain; flayout_type layout; mesh_type mesh; - std::shared_ptr playout; + std::shared_ptr playout; std::shared_ptr bunch; // Particle counts for the tests. - size_t nGather = 10; // for gather test: local particles per rank + size_t nGather = 10; // for gather test: local particles per rank size_t nScatter = static_cast(std::pow(64, Dim)); // for scatter tests // Store cell sizes (hx) for use in generating positions. T hx[Dim]; - GatherScatterTest() { } + GatherScatterTest() {} void SetUp() override { // Use a high-resolution grid (e.g. 512 cells per dimension) @@ -91,13 +110,13 @@ class GatherScatterTest>> : public ::testing: std::array isParallel; isParallel.fill(true); auto owned_tu = std::make_from_tuple>(owned); - layout = flayout_type(MPI_COMM_WORLD, owned_tu, isParallel); + layout = flayout_type(MPI_COMM_WORLD, owned_tu, isParallel); ippl::Vector hx_vec; ippl::Vector origin; for (size_t d = 0; d < Dim; d++) { hx_vec[d] = domain[d] / nPoints[d]; - hx[d] = hx_vec[d]; // store cell size for distribution + hx[d] = hx_vec[d]; // store cell size for distribution origin[d] = 0; } mesh = mesh_type(owned_tu, hx_vec, origin); @@ -136,32 +155,42 @@ class GatherScatterTest>> : public ::testing: Kokkos::deep_copy(bunch->Q.getView(), Q_host); ippl::Comm->barrier(); } + + void fillAttributeQFloat(float value) { + auto Q_host = bunch->QFloat.getHostMirror(); + for (size_t i = 0; i < Q_host.size(); ++i) { + Q_host(i) = value; + } + Kokkos::deep_copy(bunch->QFloat.getView(), Q_host); + ippl::Comm->barrier(); + } }; -using TestTypes = ::testing::Types< - Parameters>, - Parameters>, - Parameters>//, - //Parameters>, - //Parameters>, - //Parameters> ->; +using TestTypes = ::testing::Types>, + Parameters>, + Parameters> //, + // Parameters>, + // Parameters>, + // Parameters> + >; TYPED_TEST_SUITE(GatherScatterTest, TestTypes); // -// GatherTest: +// GatherTest: // First, set each local Q to 10.0. -// Then, call gather with addToAttribute = false so that Q becomes 1.0 (should replace value with 1.0). -// If Q != 0, then the values were 1. not replaced and 2. not correctly gathered from the field. -// Note: for a constant field, there should not be an error during linear interpolation. +// Then, call gather with addToAttribute = false so that Q becomes 1.0 (should replace value +// with 1.0). If Q != 0, then the values were 1. not replaced and 2. not correctly gathered from the +// field. Note: for a constant field, there should not be an error during linear interpolation. // TYPED_TEST(GatherScatterTest, GatherTestReplace) { const size_t n = this->nGather; this->fillRandomPositions(n); this->fillAttributeQ(10.0); - using Mesh_t = typename TestFixture::mesh_type; - using FieldType = ippl::Field; + using Mesh_t = typename TestFixture::mesh_type; + using FieldType = + ippl::Field; FieldType field; field.initialize(this->mesh, this->layout); field = 1.0; @@ -178,7 +207,7 @@ TYPED_TEST(GatherScatterTest, GatherTestReplace) { } // -// GatherTest: +// GatherTest: // First, set each local Q to 1.0. // Then, call gather with addToAttribute = true so that Q becomes 2.0 (should add 1.0 per particle). // @@ -187,8 +216,10 @@ TYPED_TEST(GatherScatterTest, GatherTestIncrement) { this->fillRandomPositions(n); this->fillAttributeQ(1.0); - using Mesh_t = typename TestFixture::mesh_type; - using FieldType = ippl::Field; + using Mesh_t = typename TestFixture::mesh_type; + using FieldType = + ippl::Field; FieldType field; field.initialize(this->mesh, this->layout); field = 1.0; @@ -216,46 +247,80 @@ TYPED_TEST(GatherScatterTest, ScatterSimpleTest) { this->fillAttributeQ(1.0); // Create and initialize a field/mesh. - using Mesh_t = typename TestFixture::mesh_type; - using FieldType = ippl::Field; + using Mesh_t = typename TestFixture::mesh_type; + using FieldType = + ippl::Field; FieldType field; field.initialize(this->mesh, this->layout); - + field = 0.0; // Perform the simple scatter operation (extended functionality is tested below). scatter(this->bunch->Q, field, this->bunch->R); // Compute the total charge in the field and from the particles. - double total_field = field.sum(); + double total_field = field.sum(); double total_particles = this->bunch->Q.sum(); // Check that the scattered field conserves charge. ASSERT_NEAR(total_field, total_particles, 1e-6); } +// +// ScatterMixedValueTypeTest: +// Scatter a lower-precision particle attribute into a higher-precision field. +// This covers ParticleAttrib -> Field, which is the mixed type +// use case enabled by ParticleAttrib::scatter accepting a field value type that +// differs from the attribute value type. +// +TYPED_TEST(GatherScatterTest, ScatterMixedValueTypeTest) { + using ScatterType = float; + using FieldType = double; + + const unsigned int n = this->nScatter; + this->fillRandomPositions(n); + this->fillAttributeQFloat(ScatterType(1.0)); + + using Mesh_t = typename TestFixture::mesh_type; + using Field = ippl::Field; + Field field; + field.initialize(this->mesh, this->layout); + field = FieldType(0.0); + + scatter(this->bunch->QFloat, field, this->bunch->R); + + const double total_field = field.sum(); + const double total_particles = static_cast(this->bunch->QFloat.sum()); + const double tolerance = scatterConservationTolerance(total_particles); + + ASSERT_NEAR(total_field, total_particles, tolerance); +} // -// ScatterCustomRangeTest: +// ScatterCustomRangeTest: // Set Q = 1.0 for all particles and scatter only a subset defined by a custom range policy. // Then compare the total charge in the field to the expected value. // TYPED_TEST(GatherScatterTest, ScatterCustomRangeTest) { const size_t n = this->nScatter; - if(n % ippl::Comm->size() != 0) { + if (n % ippl::Comm->size() != 0) { GTEST_SKIP() << "nScatter not divisible by number of ranks."; } this->fillRandomPositions(n); this->fillAttributeQ(1.0); - using Mesh_t = typename TestFixture::mesh_type; - using FieldType = ippl::Field; + using Mesh_t = typename TestFixture::mesh_type; + using FieldType = + ippl::Field; FieldType field; field.initialize(this->mesh, this->layout); field = 0.0; - size_t rank = ippl::Comm->rank(); - size_t nLoc = this->bunch->getLocalNum(); + size_t rank = ippl::Comm->rank(); + size_t nLoc = this->bunch->getLocalNum(); size_t NScattered = nLoc / 2 + rank; double Q_total = 1.0 * NScattered; @@ -269,7 +334,7 @@ TYPED_TEST(GatherScatterTest, ScatterCustomRangeTest) { } // -// ScatterCustomHashTest: +// ScatterCustomHashTest: // Assign random charges (in [0.5, 1.5]), create and shuffle an index array, // use it as a custom hash, scatter the first NScattered particles accordingly, // and compare the field’s total charge to the expected total. @@ -280,10 +345,10 @@ TYPED_TEST(GatherScatterTest, ScatterCustomHashTest) { GTEST_SKIP() << "nScatter not divisible by number of ranks."; } this->fillRandomPositions(n); - + size_t rank = ippl::Comm->rank(); - size_t nLoc = this->bunch->getLocalNum(); // since update() might change number of particles - size_t NScattered = nLoc / 2 + rank; // can be anything + size_t nLoc = this->bunch->getLocalNum(); // since update() might change number of particles + size_t NScattered = nLoc / 2 + rank; // can be anything // Assign random charges to particles std::mt19937_64 eng(42); @@ -295,8 +360,10 @@ TYPED_TEST(GatherScatterTest, ScatterCustomHashTest) { Kokkos::deep_copy(this->bunch->Q.getView(), Q_host); // Create and initialize a field/mesh. - using Mesh_t = typename TestFixture::mesh_type; - using FieldType = ippl::Field; + using Mesh_t = typename TestFixture::mesh_type; + using FieldType = + ippl::Field; FieldType field; field.initialize(this->mesh, this->layout); field = 0.0; @@ -315,15 +382,16 @@ TYPED_TEST(GatherScatterTest, ScatterCustomHashTest) { } Kokkos::deep_copy(hash, hash_host); - // First compute the total charge of the first NScattered particles as determined by the hash map + // First compute the total charge of the first NScattered particles as determined by the hash + // map double Q_total = 0.0; - auto viewQ = this->bunch->Q.getView(); + auto viewQ = this->bunch->Q.getView(); - ComputeTotalChargeLambda lambda(viewQ, hash); - Kokkos::parallel_reduce("computeTotalCharge", - Kokkos::RangePolicy(0, NScattered), - lambda, Q_total); - /*Kokkos::parallel_reduce("computeTotalCharge", + ComputeTotalChargeLambda lambda(viewQ, hash); + Kokkos::parallel_reduce("computeTotalCharge", + Kokkos::RangePolicy(0, NScattered), + lambda, Q_total); + /*Kokkos::parallel_reduce("computeTotalCharge", Kokkos::RangePolicy(0, NScattered), KOKKOS_LAMBDA(const size_t i, double& val) { val += viewQ(hash(i)); @@ -339,6 +407,69 @@ TYPED_TEST(GatherScatterTest, ScatterCustomHashTest) { ASSERT_NEAR(Q_total, Total_charge_field, 1e-6); } +// +// ScatterMixedValueTypeCustomHashTest: +// Scatter ParticleAttrib into Field through a custom hash array. +// This verifies that both the mixed value type and the hash-dispatched scatter +// path conserve the charge represented by the lower-precision attribute. +// +TYPED_TEST(GatherScatterTest, ScatterMixedValueTypeCustomHashTest) { + using ScatterType = float; + using FieldType = double; + + const size_t n = this->nScatter / ippl::Comm->size(); + if (this->nScatter % ippl::Comm->size() > 0) { + GTEST_SKIP() << "nScatter not divisible by number of ranks."; + } + this->fillRandomPositions(n); + + const size_t rank = ippl::Comm->rank(); + const size_t nLoc = this->bunch->getLocalNum(); + const size_t NScattered = nLoc / 2 + rank; + + std::mt19937_64 eng(42); + std::uniform_real_distribution unif_charge(0.5f, 1.5f); + auto Q_host = this->bunch->QFloat.getHostMirror(); + for (size_t i = 0; i < nLoc; ++i) { + Q_host(i) = unif_charge(eng); + } + Kokkos::deep_copy(this->bunch->QFloat.getView(), Q_host); + + using Mesh_t = typename TestFixture::mesh_type; + using Field = ippl::Field; + Field field; + field.initialize(this->mesh, this->layout); + field = FieldType(0.0); + + using hash_type = typename TestFixture::bunch_type::mixed_charge_container_type::hash_type; + hash_type hash("mixedIndexArray", nLoc); + std::vector host_indices(nLoc); + std::iota(host_indices.begin(), host_indices.end(), 0); + std::shuffle(host_indices.begin(), host_indices.end(), eng); + + auto hash_host = Kokkos::create_mirror_view(hash); + for (size_t i = 0; i < nLoc; ++i) { + hash_host(i) = host_indices[i]; + } + Kokkos::deep_copy(hash, hash_host); + + double Q_total = 0.0; + auto viewQ = this->bunch->QFloat.getView(); + ComputeTotalChargeLambda lambda(viewQ, hash); + Kokkos::parallel_reduce("computeMixedTotalCharge", + Kokkos::RangePolicy(0, NScattered), + lambda, Q_total); + ippl::Comm->allreduce(Q_total, 1, std::plus()); + + Kokkos::RangePolicy policy(0, NScattered); + scatter(this->bunch->QFloat, field, this->bunch->R, policy, hash); + + const double Total_charge_field = field.sum(); + const double tolerance = scatterConservationTolerance(Q_total); + ASSERT_NEAR(Q_total, Total_charge_field, tolerance); +} + int main(int argc, char* argv[]) { ippl::initialize(argc, argv); int result = 1; From 1ec9288f32b8ab382026e64127d93faf526e9669 Mon Sep 17 00:00:00 2001 From: John Biddiscombe Date: Mon, 14 Sep 2026 19:33:18 +0200 Subject: [PATCH 3/6] fix(scatter): remove dead code and broken timer start in ParticleAttrib::scatter The refactor moved the kernel implementation into detail::particleAttribScatterImpl, but ParticleAttrib::scatter was left with the timer start and a full set of unused locals (view, mesh, dx, origin, invdx, layout, lDom, nghost). The leftover IpplTimings::startTimer() was never stopped, so the scatter timer was double-started and only stopped once inside the detail helper, corrupting timing data. Remove the unused setup code and leave only the hash-view validation and dispatch. --- src/Particle/ParticleAttrib.hpp | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/src/Particle/ParticleAttrib.hpp b/src/Particle/ParticleAttrib.hpp index 03e4d56222..dd2d106aed 100644 --- a/src/Particle/ParticleAttrib.hpp +++ b/src/Particle/ParticleAttrib.hpp @@ -199,29 +199,6 @@ namespace ippl { void ParticleAttrib::scatter( Field& f, const ParticleAttrib, Properties...>& pp, policy_type iteration_policy, hash_type hash_array) const { - constexpr unsigned Dim = Field::dim; - using PositionType = typename Field::Mesh_t::value_type; - - static IpplTimings::TimerRef scatterTimer = IpplTimings::getTimer("scatter"); - IpplTimings::startTimer(scatterTimer); - using view_type = typename Field::view_type; - view_type view = f.getView(); - - using mesh_type = typename Field::Mesh_t; - const mesh_type& mesh = f.get_mesh(); - - using vector_type = typename mesh_type::vector_type; - using value_type = typename ParticleAttrib::value_type; - - const vector_type& dx = mesh.getMeshSpacing(); - const vector_type& origin = mesh.getOrigin(); - const vector_type invdx = 1.0 / dx; - - const FieldLayout& layout = f.getLayout(); - const NDIndex& lDom = layout.getLocalNDIndex(); - const int nghost = f.getNghost(); - - // using policy_type = Kokkos::RangePolicy; const bool useHashView = hash_array.extent(0) > 0; if (useHashView && std::cmp_greater(iteration_policy.end(), hash_array.extent(0))) { Inform m("scatter"); From cbae4330fb696209eef743c3aa43e1e46a96783e Mon Sep 17 00:00:00 2001 From: John Biddiscombe Date: Mon, 14 Sep 2026 19:33:23 +0200 Subject: [PATCH 4/6] fix(cic): use Val(1) as default argument for scatterToField With Val != T (e.g. float attribute -> double field), the default Val val = T(1) forces an implicit conversion from the mesh/weight type T to the scattered value type Val. Use Val(1) so the default is constructed in the value type directly, consistent with the mixed-value scatter design. --- src/Interpolation/CIC.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Interpolation/CIC.h b/src/Interpolation/CIC.h index 80e7b3a05f..86f7dbe6f0 100644 --- a/src/Interpolation/CIC.h +++ b/src/Interpolation/CIC.h @@ -80,7 +80,7 @@ namespace ippl { KOKKOS_INLINE_FUNCTION constexpr void scatterToField( const std::index_sequence&, const View& view, const Vector& wlo, const Vector& whi, - const Vector& args, Val val = T(1)); + const Vector& args, Val val = Val(1)); /*! * Gathers from a field at a single point From 26a7c9ff02c54a03406ed7462620c80e08acca0b Mon Sep 17 00:00:00 2001 From: John Biddiscombe Date: Mon, 14 Sep 2026 19:33:25 +0200 Subject: [PATCH 5/6] style(tests): fix naming conventions in TestScatter.cpp Rename new public member variables Q1/Q2 to QFloat_m/QDouble_m and the new local Total_charge_field to totalChargeField, per AGENTS.md naming rules. --- test/particle/TestScatter.cpp | 34 ++++++++++++++++++---------------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/test/particle/TestScatter.cpp b/test/particle/TestScatter.cpp index 8ca2403a13..e9e726b68e 100644 --- a/test/particle/TestScatter.cpp +++ b/test/particle/TestScatter.cpp @@ -6,17 +6,17 @@ template struct Bunch : public ippl::ParticleBase { Bunch(PLayout& playout) : ippl::ParticleBase(playout) { - this->addAttribute(Q1); - this->addAttribute(Q2); + this->addAttribute(QFloat_m); + this->addAttribute(QDouble_m); } ~Bunch() {} typedef ippl::ParticleAttrib charge_container_typeF; - charge_container_typeF Q1; + charge_container_typeF QFloat_m; typedef ippl::ParticleAttrib charge_container_typeD; - charge_container_typeD Q2; + charge_container_typeD QDouble_m; }; int main(int argc, char* argv[]) { @@ -89,41 +89,43 @@ int main(int argc, char* argv[]) { std::cout << "Sum coord: " << global_sum_coord << std::endl; } - bunch.Q1 = 1.0; + bunch.QFloat_m = 1.0; bunch.update(); field = 0.0; - scatter(bunch.Q1, field, bunch.R); + scatter(bunch.QFloat_m, field, bunch.R); // Check charge conservation try { - double Total_charge_field = field.sum(); + double totalChargeField = field.sum(); - std::cout << "Float:: Total charge in the field:" << Total_charge_field << std::endl; - std::cout << "Float:: Total charge of the particles:" << bunch.Q1.sum() << std::endl; - std::cout << "Float:: Error:" << std::fabs(bunch.Q1.sum() - Total_charge_field) + std::cout << "Float:: Total charge in the field:" << totalChargeField << std::endl; + std::cout << "Float:: Total charge of the particles:" << bunch.QFloat_m.sum() + << std::endl; + std::cout << "Float:: Error:" << std::fabs(bunch.QFloat_m.sum() - totalChargeField) << std::endl; } catch (const std::exception& e) { std::cout << e.what() << std::endl; } - bunch.Q2 = 1.0; + bunch.QDouble_m = 1.0; bunch.update(); field = 0.0; - scatter(bunch.Q2, field, bunch.R); + scatter(bunch.QDouble_m, field, bunch.R); // Check charge conservation try { - double Total_charge_field = field.sum(); + double totalChargeField = field.sum(); - std::cout << "Double:: Total charge in the field:" << Total_charge_field << std::endl; - std::cout << "Double:: Total charge of the particles:" << bunch.Q2.sum() << std::endl; - std::cout << "Double:: Error:" << std::fabs(bunch.Q2.sum() - Total_charge_field) + std::cout << "Double:: Total charge in the field:" << totalChargeField << std::endl; + std::cout << "Double:: Total charge of the particles:" << bunch.QDouble_m.sum() + << std::endl; + std::cout << "Double:: Error:" << std::fabs(bunch.QDouble_m.sum() - totalChargeField) << std::endl; } catch (const std::exception& e) { std::cout << e.what() << std::endl; From 85cb26f2c6ef5af42391ceb29e5a4409a77fe73d Mon Sep 17 00:00:00 2001 From: John Biddiscombe Date: Mon, 14 Sep 2026 19:33:28 +0200 Subject: [PATCH 6/6] style(tests): fix naming conventions in GatherScatterTest.cpp Rename the new public member QFloat to QFloat_m and the new local Total_charge_field to totalChargeField, per AGENTS.md naming rules. --- unit_tests/Particle/GatherScatterTest.cpp | 26 +++++++++++------------ 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/unit_tests/Particle/GatherScatterTest.cpp b/unit_tests/Particle/GatherScatterTest.cpp index 387ce03ec0..2f42ddf969 100644 --- a/unit_tests/Particle/GatherScatterTest.cpp +++ b/unit_tests/Particle/GatherScatterTest.cpp @@ -53,7 +53,7 @@ struct Bunch : public ippl::ParticleBase { Bunch(PLayout& playout) : ippl::ParticleBase(playout) { this->addAttribute(Q); - this->addAttribute(QFloat); + this->addAttribute(QFloat_m); } ~Bunch() = default; @@ -62,7 +62,7 @@ struct Bunch : public ippl::ParticleBase { typedef ippl::ParticleAttrib mixed_charge_container_type; charge_container_type Q; - mixed_charge_container_type QFloat; + mixed_charge_container_type QFloat_m; }; template @@ -157,11 +157,11 @@ class GatherScatterTest>> : public ::testing: } void fillAttributeQFloat(float value) { - auto Q_host = bunch->QFloat.getHostMirror(); + auto Q_host = bunch->QFloat_m.getHostMirror(); for (size_t i = 0; i < Q_host.size(); ++i) { Q_host(i) = value; } - Kokkos::deep_copy(bunch->QFloat.getView(), Q_host); + Kokkos::deep_copy(bunch->QFloat_m.getView(), Q_host); ippl::Comm->barrier(); } }; @@ -289,10 +289,10 @@ TYPED_TEST(GatherScatterTest, ScatterMixedValueTypeTest) { field.initialize(this->mesh, this->layout); field = FieldType(0.0); - scatter(this->bunch->QFloat, field, this->bunch->R); + scatter(this->bunch->QFloat_m, field, this->bunch->R); const double total_field = field.sum(); - const double total_particles = static_cast(this->bunch->QFloat.sum()); + const double total_particles = static_cast(this->bunch->QFloat_m.sum()); const double tolerance = scatterConservationTolerance(total_particles); ASSERT_NEAR(total_field, total_particles, tolerance); @@ -429,11 +429,11 @@ TYPED_TEST(GatherScatterTest, ScatterMixedValueTypeCustomHashTest) { std::mt19937_64 eng(42); std::uniform_real_distribution unif_charge(0.5f, 1.5f); - auto Q_host = this->bunch->QFloat.getHostMirror(); + auto Q_host = this->bunch->QFloat_m.getHostMirror(); for (size_t i = 0; i < nLoc; ++i) { Q_host(i) = unif_charge(eng); } - Kokkos::deep_copy(this->bunch->QFloat.getView(), Q_host); + Kokkos::deep_copy(this->bunch->QFloat_m.getView(), Q_host); using Mesh_t = typename TestFixture::mesh_type; using Field = ippl::Fieldbunch->QFloat.getView(); + auto viewQ = this->bunch->QFloat_m.getView(); ComputeTotalChargeLambda lambda(viewQ, hash); Kokkos::parallel_reduce("computeMixedTotalCharge", Kokkos::RangePolicy(0, NScattered), @@ -463,11 +463,11 @@ TYPED_TEST(GatherScatterTest, ScatterMixedValueTypeCustomHashTest) { ippl::Comm->allreduce(Q_total, 1, std::plus()); Kokkos::RangePolicy policy(0, NScattered); - scatter(this->bunch->QFloat, field, this->bunch->R, policy, hash); + scatter(this->bunch->QFloat_m, field, this->bunch->R, policy, hash); - const double Total_charge_field = field.sum(); - const double tolerance = scatterConservationTolerance(Q_total); - ASSERT_NEAR(Q_total, Total_charge_field, tolerance); + const double totalChargeField = field.sum(); + const double tolerance = scatterConservationTolerance(Q_total); + ASSERT_NEAR(Q_total, totalChargeField, tolerance); } int main(int argc, char* argv[]) {