Branch data 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
|
