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
|
