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 252 : void ParticleAttrib<T, Properties...>::create(size_type n) {
26 252 : size_type required = *(this->localNum_mp) + n;
27 252 : if (this->size() < required) {
28 252 : int overalloc = Comm->getDefaultOverallocation();
29 252 : this->realloc(required * overalloc);
30 : }
31 252 : }
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 48 : ParticleAttrib<T, Properties...>& ParticleAttrib<T, Properties...>::operator=(T x) {
82 : using policy_type = Kokkos::RangePolicy<execution_space>;
83 48 : Kokkos::parallel_for(
84 48 : "ParticleAttrib::operator=()", policy_type(0, *(this->localNum_mp)),
85 3936 : KOKKOS_CLASS_LAMBDA(const size_t i) { dview_m(i) = x; });
86 48 : 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 24 : 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 24 : constexpr unsigned Dim = Field::dim;
110 : using PositionType = typename Field::Mesh_t::value_type;
111 :
112 24 : static IpplTimings::TimerRef scatterTimer = IpplTimings::getTimer("scatter");
113 24 : IpplTimings::startTimer(scatterTimer);
114 : using view_type = typename Field::view_type;
115 24 : view_type view = f.getView();
116 :
117 : using mesh_type = typename Field::Mesh_t;
118 24 : 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 24 : const vector_type& dx = mesh.getMeshSpacing();
124 24 : const vector_type& origin = mesh.getOrigin();
125 24 : const vector_type invdx = 1.0 / dx;
126 :
127 24 : const FieldLayout<Dim>& layout = f.getLayout();
128 24 : const NDIndex<Dim>& lDom = layout.getLocalNDIndex();
129 24 : const int nghost = f.getNghost();
130 :
131 : //using policy_type = Kokkos::RangePolicy<execution_space>;
132 24 : const bool useHashView = hash_array.extent(0) > 0;
133 24 : 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 24 : Kokkos::parallel_for(
139 : "ParticleAttrib::scatter", iteration_policy,
140 1968 : KOKKOS_CLASS_LAMBDA(const size_t idx) {
141 : // map index to possible hash_map
142 1920 : size_t mapped_idx = useHashView ? hash_array(idx) : idx;
143 :
144 : // find nearest grid point
145 1920 : vector_type l = (pp(mapped_idx) - origin) * invdx + 0.5;
146 1920 : Vector<int, Field::dim> index = l;
147 1920 : Vector<PositionType, Field::dim> whi = l - index;
148 1920 : Vector<PositionType, Field::dim> wlo = 1.0 - whi;
149 :
150 1920 : Vector<size_t, Field::dim> args = index - lDom.first() + nghost;
151 :
152 : // scatter
153 1920 : const value_type& val = dview_m(mapped_idx);
154 1920 : detail::scatterToField(std::make_index_sequence<1 << Field::dim>{}, view, wlo, whi,
155 : args, val);
156 1920 : });
157 24 : IpplTimings::stopTimer(scatterTimer);
158 :
159 24 : static IpplTimings::TimerRef accumulateHaloTimer = IpplTimings::getTimer("accumulateHalo");
160 24 : IpplTimings::startTimer(accumulateHaloTimer);
161 24 : f.accumulateHalo();
162 24 : IpplTimings::stopTimer(accumulateHaloTimer);
163 24 : }
164 :
165 : template <typename T, class... Properties>
166 : template <typename Field, typename P2>
167 12 : void ParticleAttrib<T, Properties...>::gather(
168 : Field& f, const ParticleAttrib<Vector<P2, Field::dim>, Properties...>& pp,
169 : const bool addToAttribute) {
170 12 : constexpr unsigned Dim = Field::dim;
171 : using PositionType = typename Field::Mesh_t::value_type;
172 :
173 12 : static IpplTimings::TimerRef fillHaloTimer = IpplTimings::getTimer("fillHalo");
174 12 : IpplTimings::startTimer(fillHaloTimer);
175 12 : f.fillHalo();
176 12 : IpplTimings::stopTimer(fillHaloTimer);
177 :
178 12 : static IpplTimings::TimerRef gatherTimer = IpplTimings::getTimer("gather");
179 12 : IpplTimings::startTimer(gatherTimer);
180 12 : const typename Field::view_type view = f.getView();
181 :
182 : using mesh_type = typename Field::Mesh_t;
183 12 : const mesh_type& mesh = f.get_mesh();
184 :
185 : using vector_type = typename mesh_type::vector_type;
186 :
187 12 : const vector_type& dx = mesh.getMeshSpacing();
188 12 : const vector_type& origin = mesh.getOrigin();
189 12 : const vector_type invdx = 1.0 / dx;
190 :
191 12 : const FieldLayout<Dim>& layout = f.getLayout();
192 12 : const NDIndex<Dim>& lDom = layout.getLocalNDIndex();
193 12 : const int nghost = f.getNghost();
194 :
195 : using policy_type = Kokkos::RangePolicy<execution_space>;
196 12 : Kokkos::parallel_for(
197 12 : "ParticleAttrib::gather", policy_type(0, *(this->localNum_mp)),
198 408 : KOKKOS_CLASS_LAMBDA(const size_t idx) {
199 : // find nearest grid point
200 384 : vector_type l = (pp(idx) - origin) * invdx + 0.5;
201 384 : Vector<int, Field::dim> index = l;
202 384 : Vector<PositionType, Field::dim> whi = l - index;
203 384 : Vector<PositionType, Field::dim> wlo = 1.0 - whi;
204 :
205 384 : 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 384 : view, wlo, whi, args);
210 384 : if (addToAttribute) {
211 0 : dview_m(idx) += gathered;
212 : } else {
213 768 : dview_m(idx) = gathered;
214 : }
215 384 : });
216 12 : IpplTimings::stopTimer(gatherTimer);
217 12 : }
218 :
219 : /*
220 : * Non-class function
221 : *
222 : */
223 :
224 : /**
225 : * @brief Non-class interface for scattering particle attribute data onto a field.
226 : *
227 : * This overload preserves legacy functionality by providing a default iteration policy.
228 : * It calls the member scatter() with a default Kokkos::RangePolicy.
229 : *
230 : * @note The default behaviour is to scatter all particles without any custom index mapping.
231 : *
232 : * @tparam Attrib1 The type of the particle attribute.
233 : * @tparam Field The type of the field.
234 : * @tparam Attrib2 The type of the particle position attribute.
235 : * @tparam policy_type (Default: `Kokkos::RangePolicy<typename Field::execution_space>`)
236 : * @param attrib The particle attribute to scatter.
237 : * @param f The field onto which the data is scattered.
238 : * @param pp The ParticleAttrib representing particle positions.
239 : */
240 : template <typename Attrib1, typename Field, typename Attrib2,
241 : typename policy_type = Kokkos::RangePolicy<typename Field::execution_space>>
242 24 : inline void scatter(const Attrib1& attrib, Field& f, const Attrib2& pp) {
243 24 : attrib.scatter(f, pp, policy_type(0, attrib.getParticleCount()));
244 24 : }
245 :
246 : /**
247 : * @brief Non-class interface for scattering with a custom iteration policy and optional index array.
248 : *
249 : * This overload allows the caller to specify a custom `Kokkos::range_policy` and an optional
250 : * `ippl::hash_type` array. It forwards the parameters to the member scatter() function.
251 : *
252 : * @note See ParticleAttrib::scatter() for more information on the custom iteration functionality.
253 : *
254 : * @tparam Attrib1 The type of the particle attribute.
255 : * @tparam Field The type of the field.
256 : * @tparam Attrib2 The type of the particle position attribute.
257 : * @tparam policy_type (Default: `Kokkos::RangePolicy<typename Field::execution_space>`)
258 : * @param attrib The particle attribute to scatter.
259 : * @param f The field onto which the data is scattered.
260 : * @param pp The ParticleAttrib representing particle positions.
261 : * @param iteration_policy A custom `Kokkos::range_policy` defining the iteration range.
262 : * @param hash_array An optional `ippl::hash_type` array for index mapping.
263 : */
264 : template <typename Attrib1, typename Field, typename Attrib2,
265 : typename policy_type = Kokkos::RangePolicy<typename Field::execution_space>>
266 : inline void scatter(const Attrib1& attrib, Field& f, const Attrib2& pp,
267 : policy_type iteration_policy, typename Attrib1::hash_type hash_array = {}) {
268 : attrib.scatter(f, pp, iteration_policy, hash_array);
269 : }
270 :
271 : /**
272 : * @brief Non-class interface for gathering field data into a particle attribute.
273 : *
274 : * This interface calls the member ParticleAttrib::gather() function with the provided
275 : * parameters and preserving legacy behavior by assigning `addToAttribute` a default value.
276 : *
277 : * @note See ParticleAttrib::gather() for more information on the behavior of `addToAttribute`.
278 : *
279 : * @tparam Attrib1 The type of the particle attribute.
280 : * @tparam Field The type of the field.
281 : * @tparam Attrib2 The type of the particle position attribute.
282 : * @param attrib The particle attribute to gather data into.
283 : * @param f The field from which data is gathered.
284 : * @param pp The ParticleAttrib representing particle positions.
285 : * @param addToAttribute If true, the gathered field value is added to the current attribute value;
286 : * otherwise, the attribute value is overwritten.
287 : */
288 : template <typename Attrib1, typename Field, typename Attrib2>
289 12 : inline void gather(Attrib1& attrib, Field& f, const Attrib2& pp,
290 : const bool addToAttribute = false) {
291 12 : attrib.gather(f, pp, addToAttribute);
292 12 : }
293 :
294 : #define DefineParticleReduction(fun, name, op, MPI_Op) \
295 : template <typename T, class... Properties> \
296 : T ParticleAttrib<T, Properties...>::name() { \
297 : T temp = 0.0; \
298 : using policy_type = Kokkos::RangePolicy<execution_space>; \
299 : Kokkos::parallel_reduce( \
300 : "fun", policy_type(0, *(this->localNum_mp)), \
301 : KOKKOS_CLASS_LAMBDA(const size_t i, T& valL) { \
302 : T myVal = dview_m(i); \
303 : op; \
304 : }, \
305 : Kokkos::fun<T>(temp)); \
306 : T globaltemp = 0.0; \
307 : Comm->allreduce(temp, globaltemp, 1, MPI_Op<T>()); \
308 : return globaltemp; \
309 : }
310 :
311 432 : DefineParticleReduction(Sum, sum, valL += myVal, std::plus)
312 : DefineParticleReduction(Max, max, if (myVal > valL) valL = myVal, std::greater)
313 : DefineParticleReduction(Min, min, if (myVal < valL) valL = myVal, std::less)
314 : DefineParticleReduction(Prod, prod, valL *= myVal, std::multiplies)
315 : } // namespace ippl
|
