Refactor ParticleAttrib::scatter()

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 <unsigned long... ScatterPoint, typename View, typename T,
-                  typename IndexType = size_t>
+                  typename IndexType = size_t, typename Val = T>
         KOKKOS_INLINE_FUNCTION constexpr void scatterToField(
             const std::index_sequence<ScatterPoint...>&, const View& view,
             const Vector<T, View::rank>& wlo, const Vector<T, View::rank>& whi,
-            const Vector<IndexType, View::rank>& args, T val = 1);
+            const Vector<IndexType, View::rank>& args, Val val = Val(1));
 
         /*!
          * Gathers from a field at a single point

src/Interpolation/CIC.hpp

@@ -1,5 +1,6 @@
 #include <Kokkos_Core.hpp>
 
+#include <type_traits>
 
 namespace ippl {
     namespace detail {
@@ -33,14 +34,18 @@ namespace ippl {
                                val * (interpolationWeight<ScatterPoint, Index>(wlo, whi) * ...));
         }
 
-        template <unsigned long... ScatterPoint, typename View, typename T, typename IndexType>
+        template <unsigned long... ScatterPoint, typename View, typename T, typename IndexType,
+                  typename Val>
         KOKKOS_INLINE_FUNCTION constexpr void scatterToField(
             const std::index_sequence<ScatterPoint...>&, const View& view,
             const Vector<T, View::rank>& wlo, const Vector<T, View::rank>& whi,
-            const Vector<IndexType, View::rank>& args, T val) {
+            const Vector<IndexType, View::rank>& args, Val val) {
             // The number of indices is equal to the view rank
+            using out_type =
+                std::remove_cv_t<std::remove_reference_t<typename View::non_const_value_type>>;
+            const out_type field_value = static_cast<out_type>(val);
             (scatterToPoint<ScatterPoint>(std::make_index_sequence<View::rank>{}, view, wlo, whi,
-                                          args, val),
+                                          args, field_value),
              ...);
         }
 

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<size_type>(*(this->localNum_mp))));
+            return Kokkos::subview(
+                dview_m,
+                Kokkos::make_pair(size_type(0), static_cast<size_type>(*(this->localNum_mp))));
         }
         const view_type getView() const {
-            return Kokkos::subview(dview_m,
-                                   Kokkos::make_pair(size_type(0),
-                                                     static_cast<size_type>(*(this->localNum_mp))));
+            return Kokkos::subview(
+                dview_m,
+                Kokkos::make_pair(size_type(0), static_cast<size_type>(*(this->localNum_mp))));
         }
 
         host_mirror_type getHostMirror() const { return Kokkos::create_mirror(getView()); }

src/Particle/ParticleAttrib.hpp

@@ -16,7 +16,6 @@
 #include "Ippl.h"
 
 #include <Kokkos_MathematicalConstants.hpp>
-
 #include <cmath>
 #include <utility>
 
@@ -32,6 +31,72 @@
 #include "Particle/SortBuffer.h"
 
 namespace ippl {
+    namespace detail {
+        template <bool UseHashView, typename HashView>
+        KOKKOS_INLINE_FUNCTION size_t scatterMappedIndex(const size_t idx,
+                                                         const HashView& hashView) {
+            if constexpr (UseHashView) {
+                return static_cast<size_t>(hashView(idx));
+            } else {
+                (void)hashView;
+                return idx;
+            }
+        }
+
+        template <bool UseHashView, typename Field, typename ValuesView, typename PositionAttrib,
+                  typename policy_type, typename HashView>
+        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<Dim>& layout = f.getLayout();
+            const NDIndex<Dim>& 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<UseHashView>(idx, hashView);
+
+                    vector_type l                 = (ppview(mapped_idx) - origin) * invdx + 0.5;
+                    Vector<int, Field::dim> index = l;
+                    Vector<PositionType, Field::dim> whi = l - index;
+                    Vector<PositionType, Field::dim> wlo = 1.0 - whi;
+
+                    Vector<size_t, Field::dim> 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 <typename T, class... Properties>
     void ParticleAttrib<T, Properties...>::create(size_type n, bool non_destructive) {
@@ -134,60 +199,19 @@ namespace ippl {
     void ParticleAttrib<T, Properties...>::scatter(
         Field& f, const ParticleAttrib<Vector<PT, Field::dim>, 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<T, Properties...>::value_type;
-
-        const vector_type& dx     = mesh.getMeshSpacing();
-        const vector_type& origin = mesh.getOrigin();
-        const vector_type invdx   = 1.0 / dx;
-
-        const FieldLayout<Dim>& layout = f.getLayout();
-        const NDIndex<Dim>& lDom       = layout.getLocalNDIndex();
-        const int nghost               = f.getNghost();
-
-        // using policy_type = Kokkos::RangePolicy<execution_space>;
         const bool useHashView = hash_array.extent(0) > 0;
         if (useHashView && std::cmp_greater(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;
-
-                vector_type l                        = (ppview(mapped_idx) - origin) * invdx + 0.5;
-                Vector<int, Field::dim> index        = l;
-                Vector<PositionType, Field::dim> whi = l - index;
-                Vector<PositionType, Field::dim> wlo = 1.0 - whi;
 
-                Vector<size_t, Field::dim> 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<true>(f, dview_m, pp, iteration_policy, hash_array);
+        } else {
+            detail::particleAttribScatterImpl<false>(f, dview_m, pp, iteration_policy, hash_array);
+        }
     }
 
     template <typename T, class... Properties>
@@ -411,8 +435,7 @@ namespace ippl {
         Field<FT, Dim, M, C>& f, Field<ST, Dim, M, C>& Sk,
         const ParticleAttrib<Vector<PT, Dim>, Properties...>& pp,
         FFT<NUFFTransform, Field<ST, Dim, M, C>>* nufft, ParticleAttrib<PT, Properties...>& q) {
-        static IpplTimings::TimerRef gatherPIFNUFFTTimer =
-            IpplTimings::getTimer("GatherPIFNUFFT");
+        static IpplTimings::TimerRef gatherPIFNUFFTTimer = IpplTimings::getTimer("GatherPIFNUFFT");
         IpplTimings::startTimer(gatherPIFNUFFTTimer);
 
         typename Field<FT, Dim, M, C>::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<int, 3> iVec = {i, j, k};
                     for (unsigned d = 0; d < Dim; ++d) {

test/particle/TestScatter.cpp

@@ -6,13 +6,17 @@ template <class PLayout>
 struct Bunch : public ippl::ParticleBase<PLayout> {
     Bunch(PLayout& playout)
         : ippl::ParticleBase<PLayout>(playout) {
-        this->addAttribute(Q);
+        this->addAttribute(QFloat_m);
+        this->addAttribute(QDouble_m);
     }
 
     ~Bunch() {}
 
-    typedef ippl::ParticleAttrib<double> charge_container_type;
-    charge_container_type Q;
+    typedef ippl::ParticleAttrib<float> charge_container_typeF;
+    charge_container_typeF QFloat_m;
+
+    typedef ippl::ParticleAttrib<double> 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;
 }