diff --git a/src/Interpolation/CIC.h b/src/Interpolation/CIC.h index 5adac510d3..86f7dbe6f0 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 = Val(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 6d89990081..6dcbe0b713 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()); } diff --git a/src/Particle/ParticleAttrib.hpp b/src/Particle/ParticleAttrib.hpp index 5fa34bc084..dd8056c82e 100644 --- a/src/Particle/ParticleAttrib.hpp +++ b/src/Particle/ParticleAttrib.hpp @@ -16,7 +16,6 @@ #include "Ippl.h" #include - #include #include @@ -32,6 +31,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) { @@ -134,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"); @@ -164,30 +206,12 @@ namespace ippl { << 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; - - 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); + if (useHashView) { + detail::particleAttribScatterImpl(f, dview_m, pp, iteration_policy, hash_array); + } else { + detail::particleAttribScatterImpl(f, dview_m, pp, iteration_policy, hash_array); + } } template @@ -411,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; @@ -453,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..e9e726b68e 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(QFloat_m); + this->addAttribute(QDouble_m); } ~Bunch() {} - typedef ippl::ParticleAttrib charge_container_type; - charge_container_type Q; + typedef ippl::ParticleAttrib charge_container_typeF; + charge_container_typeF QFloat_m; + + typedef ippl::ParticleAttrib charge_container_typeD; + charge_container_typeD QDouble_m; }; int main(int argc, char* argv[]) { @@ -85,26 +89,48 @@ int main(int argc, char* argv[]) { std::cout << "Sum coord: " << global_sum_coord << std::endl; } - bunch.Q = 1.0; + bunch.QFloat_m = 1.0; bunch.update(); field = 0.0; - scatter(bunch.Q, 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 << "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 << "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.QDouble_m = 1.0; + + bunch.update(); + + field = 0.0; + + scatter(bunch.QDouble_m, field, bunch.R); + + // Check charge conservation + try { + double totalChargeField = field.sum(); + + 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; } } ippl::finalize(); - return 0; } diff --git a/unit_tests/Particle/GatherScatterTest.cpp b/unit_tests/Particle/GatherScatterTest.cpp index 5a05f16189..2f42ddf969 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_m); } ~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_m; }; 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_m.getHostMirror(); + for (size_t i = 0; i < Q_host.size(); ++i) { + Q_host(i) = value; + } + Kokkos::deep_copy(bunch->QFloat_m.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_m, field, this->bunch->R); + + const double total_field = field.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); +} // -// 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_m.getHostMirror(); + for (size_t i = 0; i < nLoc; ++i) { + Q_host(i) = unif_charge(eng); + } + Kokkos::deep_copy(this->bunch->QFloat_m.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_m.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_m, field, this->bunch->R, policy, hash); + + const double totalChargeField = field.sum(); + const double tolerance = scatterConservationTolerance(Q_total); + ASSERT_NEAR(Q_total, totalChargeField, tolerance); +} + int main(int argc, char* argv[]) { ippl::initialize(argc, argv); int result = 1;