LCOV - code coverage report
Current view: top level - src/Particle - ParticleAttrib.hpp (source / functions) Coverage Total Hit
Test: final_report.info Lines: 57.1 % 140 80
Test Date: 2025-07-18 09:50:00 Functions: 19.9 % 317 63

            Line data    Source code
       1              : //
       2              : // Class ParticleAttrib
       3              : //   Templated class for all particle attribute classes.
       4              : //
       5              : //   This templated class is used to represent a single particle attribute.
       6              : //   An attribute is one data element within a particle object, and is
       7              : //   stored as a Kokkos::View. This class stores the type information for the
       8              : //   attribute, and provides methods to create and destroy new items, and
       9              : //   to perform operations involving this attribute with others.
      10              : //
      11              : //   ParticleAttrib is the primary element involved in expressions for
      12              : //   particles (just as BareField is the primary element there).  This file
      13              : //   defines the necessary templated classes and functions to make
      14              : //   ParticleAttrib a capable expression-template participant.
      15              : //
      16              : #include "Ippl.h"
      17              : 
      18              : #include "Communicate/DataTypes.h"
      19              : 
      20              : #include "Utility/IpplTimings.h"
      21              : 
      22              : namespace ippl {
      23              : 
      24              :     template <typename T, class... Properties>
      25          142 :     void ParticleAttrib<T, Properties...>::create(size_type n) {
      26          142 :         size_type required = *(this->localNum_mp) + n;
      27          142 :         if (this->size() < required) {
      28          142 :             int overalloc = Comm->getDefaultOverallocation();
      29          142 :             this->realloc(required * overalloc);
      30              :         }
      31          142 :     }
      32              : 
      33              :     template <typename T, class... Properties>
      34           24 :     void ParticleAttrib<T, Properties...>::destroy(const hash_type& deleteIndex,
      35              :                                                    const hash_type& keepIndex,
      36              :                                                    size_type invalidCount) {
      37              :         // Replace all invalid particles in the valid region with valid
      38              :         // particles in the invalid region
      39              :         using policy_type = Kokkos::RangePolicy<execution_space>;
      40           24 :         Kokkos::parallel_for(
      41           24 :             "ParticleAttrib::destroy()", policy_type(0, invalidCount),
      42        12048 :             KOKKOS_CLASS_LAMBDA(const size_t i) {
      43            0 :                 dview_m(deleteIndex(i)) = dview_m(keepIndex(i));
      44              :             });
      45           24 :     }
      46              : 
      47              :     template <typename T, class... Properties>
      48            0 :     void ParticleAttrib<T, Properties...>::pack(const hash_type& hash) {
      49            0 :         auto size = hash.extent(0);
      50            0 :         if (buf_m.extent(0) < size) {
      51            0 :             int overalloc = Comm->getDefaultOverallocation();
      52            0 :             Kokkos::realloc(buf_m, size * overalloc);
      53              :         }
      54              : 
      55              :         using policy_type = Kokkos::RangePolicy<execution_space>;
      56            0 :         Kokkos::parallel_for(
      57            0 :             "ParticleAttrib::pack()", policy_type(0, size),
      58            0 :             KOKKOS_CLASS_LAMBDA(const size_t i) { buf_m(i) = dview_m(hash(i)); });
      59            0 :         Kokkos::fence();
      60            0 :     }
      61              : 
      62              :     template <typename T, class... Properties>
      63            0 :     void ParticleAttrib<T, Properties...>::unpack(size_type nrecvs) {
      64            0 :         auto size          = dview_m.extent(0);
      65            0 :         size_type required = *(this->localNum_mp) + nrecvs;
      66            0 :         if (size < required) {
      67            0 :             int overalloc = Comm->getDefaultOverallocation();
      68            0 :             this->resize(required * overalloc);
      69              :         }
      70              : 
      71            0 :         size_type count   = *(this->localNum_mp);
      72              :         using policy_type = Kokkos::RangePolicy<execution_space>;
      73            0 :         Kokkos::parallel_for(
      74            0 :             "ParticleAttrib::unpack()", policy_type(0, nrecvs),
      75            0 :             KOKKOS_CLASS_LAMBDA(const size_t i) { dview_m(count + i) = buf_m(i); });
      76            0 :         Kokkos::fence();
      77            0 :     }
      78              : 
      79              :     template <typename T, class... Properties>
      80              :     // KOKKOS_INLINE_FUNCTION
      81            7 :     ParticleAttrib<T, Properties...>& ParticleAttrib<T, Properties...>::operator=(T x) {
      82              :         using policy_type = Kokkos::RangePolicy<execution_space>;
      83            7 :         Kokkos::parallel_for(
      84            7 :             "ParticleAttrib::operator=()", policy_type(0, *(this->localNum_mp)),
      85          718 :             KOKKOS_CLASS_LAMBDA(const size_t i) { dview_m(i) = x; });
      86            7 :         return *this;
      87              :     }
      88              : 
      89              :     template <typename T, class... Properties>
      90              :     template <typename E, size_t N>
      91              :     // KOKKOS_INLINE_FUNCTION
      92            0 :     ParticleAttrib<T, Properties...>& ParticleAttrib<T, Properties...>::operator=(
      93              :         detail::Expression<E, N> const& expr) {
      94              :         using capture_type = detail::CapturedExpression<E, N>;
      95            0 :         capture_type expr_ = reinterpret_cast<const capture_type&>(expr);
      96              : 
      97              :         using policy_type = Kokkos::RangePolicy<execution_space>;
      98            0 :         Kokkos::parallel_for(
      99            0 :             "ParticleAttrib::operator=()", policy_type(0, *(this->localNum_mp)),
     100            0 :             KOKKOS_CLASS_LAMBDA(const size_t i) { dview_m(i) = expr_(i); });
     101            0 :         return *this;
     102              :     }
     103              : 
     104              :     template <typename T, class... Properties>
     105              :     template <typename Field, class PT, typename policy_type>
     106            4 :     void ParticleAttrib<T, Properties...>::scatter(
     107              :         Field& f, const ParticleAttrib<Vector<PT, Field::dim>, Properties...>& pp,
     108              :         policy_type iteration_policy, hash_type hash_array) const {
     109            4 :         constexpr unsigned Dim = Field::dim;
     110              :         using PositionType     = typename Field::Mesh_t::value_type;
     111              : 
     112            4 :         static IpplTimings::TimerRef scatterTimer = IpplTimings::getTimer("scatter");
     113            4 :         IpplTimings::startTimer(scatterTimer);
     114              :         using view_type = typename Field::view_type;
     115            4 :         view_type view  = f.getView();
     116              : 
     117              :         using mesh_type       = typename Field::Mesh_t;
     118            4 :         const mesh_type& mesh = f.get_mesh();
     119              : 
     120              :         using vector_type = typename mesh_type::vector_type;
     121              :         using value_type  = typename ParticleAttrib<T, Properties...>::value_type;
     122              : 
     123            4 :         const vector_type& dx     = mesh.getMeshSpacing();
     124            4 :         const vector_type& origin = mesh.getOrigin();
     125            4 :         const vector_type invdx   = 1.0 / dx;
     126              : 
     127            4 :         const FieldLayout<Dim>& layout = f.getLayout();
     128            4 :         const NDIndex<Dim>& lDom       = layout.getLocalNDIndex();
     129            4 :         const int nghost               = f.getNghost();
     130              : 
     131              :         //using policy_type = Kokkos::RangePolicy<execution_space>;
     132            4 :         const bool useHashView = hash_array.extent(0) > 0;
     133            4 :         if (useHashView && (iteration_policy.end() > hash_array.extent(0))) {
     134            0 :             Inform m("scatter");
     135            0 :             m << "Hash array was passed to scatter, but size does not match iteration policy." << endl;
     136            0 :             ippl::Comm->abort();
     137            0 :         }
     138            4 :         Kokkos::parallel_for(
     139              :             "ParticleAttrib::scatter", iteration_policy,
     140          424 :             KOKKOS_CLASS_LAMBDA(const size_t idx) {
     141              :                 // map index to possible hash_map
     142          416 :                 size_t mapped_idx = useHashView ? hash_array(idx) : idx;
     143              : 
     144              :                 // find nearest grid point
     145          416 :                 vector_type l                        = (pp(mapped_idx) - origin) * invdx + 0.5;
     146          416 :                 Vector<int, Field::dim> index        = l;
     147          416 :                 Vector<PositionType, Field::dim> whi = l - index;
     148          416 :                 Vector<PositionType, Field::dim> wlo = 1.0 - whi;
     149              : 
     150          416 :                 Vector<size_t, Field::dim> args = index - lDom.first() + nghost;
     151              : 
     152              :                 // scatter
     153            0 :                 const value_type& val = dview_m(mapped_idx);
     154          416 :                 detail::scatterToField(std::make_index_sequence<1 << Field::dim>{}, view, wlo, whi,
     155              :                                        args, val);
     156            0 :             });
     157            4 :         IpplTimings::stopTimer(scatterTimer);
     158              : 
     159            4 :         static IpplTimings::TimerRef accumulateHaloTimer = IpplTimings::getTimer("accumulateHalo");
     160            4 :         IpplTimings::startTimer(accumulateHaloTimer);
     161            4 :         f.accumulateHalo();
     162            4 :         IpplTimings::stopTimer(accumulateHaloTimer);
     163            4 :     }
     164              : 
     165              :     template <typename T, class... Properties>
     166              :     template <typename Field, typename P2>
     167            1 :     void ParticleAttrib<T, Properties...>::gather(
     168              :         Field& f, const ParticleAttrib<Vector<P2, Field::dim>, Properties...>& pp,
     169              :         const bool addToAttribute) {
     170            1 :         constexpr unsigned Dim = Field::dim;
     171              :         using PositionType     = typename Field::Mesh_t::value_type;
     172              : 
     173            1 :         static IpplTimings::TimerRef fillHaloTimer = IpplTimings::getTimer("fillHalo");
     174            1 :         IpplTimings::startTimer(fillHaloTimer);
     175            1 :         f.fillHalo();
     176            1 :         IpplTimings::stopTimer(fillHaloTimer);
     177              : 
     178            1 :         static IpplTimings::TimerRef gatherTimer = IpplTimings::getTimer("gather");
     179            1 :         IpplTimings::startTimer(gatherTimer);
     180            1 :         const typename Field::view_type view = f.getView();
     181              : 
     182              :         using mesh_type       = typename Field::Mesh_t;
     183            1 :         const mesh_type& mesh = f.get_mesh();
     184              : 
     185              :         using vector_type = typename mesh_type::vector_type;
     186              : 
     187            1 :         const vector_type& dx     = mesh.getMeshSpacing();
     188            1 :         const vector_type& origin = mesh.getOrigin();
     189            1 :         const vector_type invdx   = 1.0 / dx;
     190              : 
     191            1 :         const FieldLayout<Dim>& layout = f.getLayout();
     192            1 :         const NDIndex<Dim>& lDom       = layout.getLocalNDIndex();
     193            1 :         const int nghost               = f.getNghost();
     194              : 
     195              :         using policy_type = Kokkos::RangePolicy<execution_space>;
     196            1 :         Kokkos::parallel_for(
     197            1 :             "ParticleAttrib::gather", policy_type(0, *(this->localNum_mp)),
     198           34 :             KOKKOS_CLASS_LAMBDA(const size_t idx) {
     199              :                 // find nearest grid point
     200           32 :                 vector_type l                        = (pp(idx) - origin) * invdx + 0.5;
     201           32 :                 Vector<int, Field::dim> index        = l;
     202           32 :                 Vector<PositionType, Field::dim> whi = l - index;
     203           32 :                 Vector<PositionType, Field::dim> wlo = 1.0 - whi;
     204              : 
     205           32 :                 Vector<size_t, Field::dim> args = index - lDom.first() + nghost;
     206              : 
     207              :                 // gather
     208            0 :                 value_type gathered = detail::gatherFromField(std::make_index_sequence<1 << Field::dim>{},
     209           32 :                                                               view, wlo, whi, args);
     210           32 :                 if (addToAttribute) {
     211            0 :                     dview_m(idx) += gathered;
     212              :                 } else {
     213            0 :                     dview_m(idx)  = gathered;
     214              :                 }
     215            0 :             });
     216            1 :         IpplTimings::stopTimer(gatherTimer);
     217            1 :     }
     218              : 
     219              :     template <typename T, class... Properties>
     220            0 :     void ParticleAttrib<T, Properties...>::applyPermutation(
     221              :         const hash_type& permutation) {
     222              : 
     223            0 :         const auto view = this->getView();
     224            0 :         const auto size = this->getParticleCount();
     225              : 
     226            0 :         view_type temp("copy", size);
     227              : 
     228              :         using policy_type = Kokkos::RangePolicy<execution_space>;
     229            0 :         Kokkos::parallel_for(
     230            0 :             "Copy to temp", policy_type(0, size),
     231            0 :             KOKKOS_LAMBDA(const size_type& i) { temp(permutation(i)) = view(i); });
     232              : 
     233            0 :         Kokkos::fence();
     234              : 
     235            0 :         Kokkos::deep_copy(Kokkos::subview(view, Kokkos::make_pair<size_type, size_type>(0, size)), temp);
     236            0 :     }
     237              : 
     238              :     template<typename T, class... Properties>
     239            0 :     void ParticleAttrib<T, Properties...>::internalCopy(
     240              :         const hash_type &indices) {
     241            0 :         auto copySize = indices.size();
     242            0 :         create(copySize);
     243              : 
     244            0 :         auto view = this->getView();
     245            0 :         const auto size = this->getParticleCount();
     246              : 
     247              :         using policy_type = Kokkos::RangePolicy<execution_space>;
     248            0 :         Kokkos::parallel_for(
     249            0 :             "Copy to temp", policy_type(0, copySize),
     250            0 :             KOKKOS_LAMBDA(const size_type &i) {
     251            0 :             view(size + i) = view(i);
     252              :         });
     253              : 
     254            0 :         Kokkos::fence();
     255            0 :     }
     256              : 
     257              :     /*
     258              :      * Non-class function
     259              :      *
     260              :      */
     261              : 
     262              :     /**
     263              :      * @brief Non-class interface for scattering particle attribute data onto a field.
     264              :      *
     265              :      * This overload preserves legacy functionality by providing a default iteration policy.
     266              :      * It calls the member scatter() with a default Kokkos::RangePolicy.
     267              :      * 
     268              :      * @note The default behaviour is to scatter all particles without any custom index mapping.
     269              :      *
     270              :      * @tparam Attrib1 The type of the particle attribute.
     271              :      * @tparam Field The type of the field.
     272              :      * @tparam Attrib2 The type of the particle position attribute.
     273              :      * @tparam policy_type (Default: `Kokkos::RangePolicy<typename Field::execution_space>`)
     274              :      * @param attrib The particle attribute to scatter.
     275              :      * @param f The field onto which the data is scattered.
     276              :      * @param pp The ParticleAttrib representing particle positions.
     277              :      */
     278              :     template <typename Attrib1, typename Field, typename Attrib2, 
     279              :                 typename policy_type = Kokkos::RangePolicy<typename Field::execution_space>>
     280            4 :     inline void scatter(const Attrib1& attrib, Field& f, const Attrib2& pp) {
     281            4 :         attrib.scatter(f, pp, policy_type(0, attrib.getParticleCount())); 
     282            4 :     }
     283              : 
     284              :     /**
     285              :      * @brief Non-class interface for scattering with a custom iteration policy and optional index array.
     286              :      *
     287              :      * This overload allows the caller to specify a custom `Kokkos::range_policy` and an optional
     288              :      * `ippl::hash_type` array. It forwards the parameters to the member scatter() function.
     289              :      * 
     290              :      * @note See ParticleAttrib::scatter() for more information on the custom iteration functionality.
     291              :      *
     292              :      * @tparam Attrib1 The type of the particle attribute.
     293              :      * @tparam Field The type of the field.
     294              :      * @tparam Attrib2 The type of the particle position attribute.
     295              :      * @tparam policy_type (Default: `Kokkos::RangePolicy<typename Field::execution_space>`)
     296              :      * @param attrib The particle attribute to scatter.
     297              :      * @param f The field onto which the data is scattered.
     298              :      * @param pp The ParticleAttrib representing particle positions.
     299              :      * @param iteration_policy A custom `Kokkos::range_policy` defining the iteration range.
     300              :      * @param hash_array An optional `ippl::hash_type` array for index mapping.
     301              :      */
     302              :     template <typename Attrib1, typename Field, typename Attrib2, 
     303              :                 typename policy_type = Kokkos::RangePolicy<typename Field::execution_space>>
     304              :     inline void scatter(const Attrib1& attrib, Field& f, const Attrib2& pp, 
     305              :                         policy_type iteration_policy, typename Attrib1::hash_type hash_array = {}) {
     306              :         attrib.scatter(f, pp, iteration_policy, hash_array);
     307              :     }
     308              : 
     309              :     /**
     310              :      * @brief Non-class interface for gathering field data into a particle attribute.
     311              :      * 
     312              :      * This interface calls the member ParticleAttrib::gather() function with the provided
     313              :      * parameters and preserving legacy behavior by assigning `addToAttribute` a default value.
     314              :      * 
     315              :      * @note See ParticleAttrib::gather() for more information on the behavior of `addToAttribute`.
     316              :      * 
     317              :      * @tparam Attrib1 The type of the particle attribute.
     318              :      * @tparam Field The type of the field.
     319              :      * @tparam Attrib2 The type of the particle position attribute.
     320              :      * @param attrib The particle attribute to gather data into.
     321              :      * @param f The field from which data is gathered.
     322              :      * @param pp The ParticleAttrib representing particle positions.
     323              :      * @param addToAttribute If true, the gathered field value is added to the current attribute value;
     324              :      *                       otherwise, the attribute value is overwritten.
     325              :      */
     326              :     template <typename Attrib1, typename Field, typename Attrib2>
     327            1 :     inline void gather(Attrib1& attrib, Field& f, const Attrib2& pp, 
     328              :                         const bool addToAttribute = false) {
     329            1 :         attrib.gather(f, pp, addToAttribute);
     330            1 :     }
     331              : 
     332              : #define DefineParticleReduction(fun, name, op, MPI_Op)            \
     333              :     template <typename T, class... Properties>                    \
     334              :     T ParticleAttrib<T, Properties...>::name() {                  \
     335              :         T temp            = 0.0;                                  \
     336              :         using policy_type = Kokkos::RangePolicy<execution_space>; \
     337              :         Kokkos::parallel_reduce(                                  \
     338              :             "fun", policy_type(0, *(this->localNum_mp)),          \
     339              :             KOKKOS_CLASS_LAMBDA(const size_t i, T& valL) {        \
     340              :                 T myVal = dview_m(i);                             \
     341              :                 op;                                               \
     342              :             },                                                    \
     343              :             Kokkos::fun<T>(temp));                                \
     344              :         T globaltemp = 0.0;                                       \
     345              :         Comm->allreduce(temp, globaltemp, 1, MPI_Op<T>());        \
     346              :         return globaltemp;                                        \
     347              :     }
     348              : 
     349           36 :     DefineParticleReduction(Sum, sum, valL += myVal, std::plus)
     350              :     DefineParticleReduction(Max, max, if (myVal > valL) valL = myVal, std::greater)
     351              :     DefineParticleReduction(Min, min, if (myVal < valL) valL = myVal, std::less)
     352              :     DefineParticleReduction(Prod, prod, valL *= myVal, std::multiplies)
     353              : }  // namespace ippl
        

Generated by: LCOV version 2.0-1