Branch data Line data Source code
1 : : //
2 : : // Class ParticleBase
3 : : // Base class for all user-defined particle classes.
4 : : //
5 : : // ParticleBase is a container and manager for a set of particles.
6 : : // The user must define a class derived from ParticleBase which describes
7 : : // what specific data attributes the particle has (e.g., mass or charge).
8 : : // Each attribute is an instance of a ParticleAttribute<T> class; ParticleBase
9 : : // keeps a list of pointers to these attributes, and performs particle creation
10 : : // and destruction.
11 : : //
12 : : // ParticleBase is templated on the ParticleLayout mechanism for the particles.
13 : : // This template parameter should be a class derived from ParticleLayout.
14 : : // ParticleLayout-derived classes maintain the info on which particles are
15 : : // located on which processor, and performs the specific communication
16 : : // required between processors for the particles. The ParticleLayout is
17 : : // templated on the type and dimension of the atom position attribute, and
18 : : // ParticleBase uses the same types for these items as the given
19 : : // ParticleLayout.
20 : : //
21 : : // ParticleBase and all derived classes have the following common
22 : : // characteristics:
23 : : // - The spatial positions of the N particles are stored in the
24 : : // particle_position_type variable R
25 : : // - The global index of the N particles are stored in the
26 : : // particle_index_type variable ID
27 : : // - A pointer to an allocated layout class. When you construct a
28 : : // ParticleBase, you must provide a layout instance, and ParticleBase
29 : : // will delete this instance when it (the ParticleBase) is deleted.
30 : : //
31 : : // To use this class, the user defines a derived class with the same
32 : : // structure as in this example:
33 : : //
34 : : // class UserParticles :
35 : : // public ParticleBase< ParticleSpatialLayout<double,3> > {
36 : : // public:
37 : : // // attributes for this class
38 : : // ParticleAttribute<double> rad; // radius
39 : : // particle_position_type vel; // velocity, same storage type as R
40 : : //
41 : : // // constructor: add attributes to base class
42 : : // UserParticles(ParticleSpatialLayout<double,2>* L) : ParticleBase(L) {
43 : : // addAttribute(rad);
44 : : // addAttribute(vel);
45 : : // }
46 : : // };
47 : : //
48 : : // This example defines a user class with 3D position and two extra
49 : : // attributes: a radius rad (double), and a velocity vel (a 3D Vector).
50 : : //
51 : :
52 : : namespace ippl {
53 : :
54 : : template <class PLayout, typename... IP>
55 : 204 : ParticleBase<PLayout, IP...>::ParticleBase()
56 : 204 : : layout_m(nullptr)
57 : 204 : , localNum_m(0)
58 : 204 : , totalNum_m(0)
59 : 204 : , nextID_m(Comm->rank())
60 [ + - + - : 408 : , numNodes_m(Comm->size()) {
+ - + - ]
61 : : if constexpr (EnableIDs) {
62 [ + - ]: 48 : addAttribute(ID);
63 : : }
64 [ + - ]: 204 : addAttribute(R);
65 : 204 : }
66 : :
67 : : template <class PLayout, typename... IP>
68 : 192 : ParticleBase<PLayout, IP...>::ParticleBase(PLayout& layout)
69 : 192 : : ParticleBase() {
70 : 192 : initialize(layout);
71 : 192 : }
72 : :
73 : : template <class PLayout, typename... IP>
74 : : template <typename MemorySpace>
75 : 336 : void ParticleBase<PLayout, IP...>::addAttribute(detail::ParticleAttribBase<MemorySpace>& pa) {
76 [ + - ]: 336 : attributes_m.template get<MemorySpace>().push_back(&pa);
77 : 336 : pa.setParticleCount(localNum_m);
78 : 336 : }
79 : :
80 : : template <class PLayout, typename... IP>
81 : 204 : void ParticleBase<PLayout, IP...>::initialize(PLayout& layout) {
82 : : // PAssert(layout_m == nullptr);
83 : :
84 : : // save the layout, and perform setup tasks
85 : 204 : layout_m = &layout;
86 : 204 : }
87 : :
88 : : template <class PLayout, typename... IP>
89 : 168 : void ParticleBase<PLayout, IP...>::create(size_type nLocal) {
90 [ - + ]: 168 : PAssert(layout_m != nullptr);
91 : :
92 [ + - ]: 168 : if (nLocal > 0) {
93 [ + - ]: 672 : forAllAttributes([&]<typename Attribute>(Attribute& attribute) {
94 : 252 : attribute->create(nLocal);
95 : : });
96 : :
97 : : if constexpr (EnableIDs) {
98 : : // set the unique ID value for these new particles
99 : : using policy_type =
100 : : Kokkos::RangePolicy<size_type, typename particle_index_type::execution_space>;
101 [ + - ]: 12 : auto pIDs = ID.getView();
102 : 12 : auto nextID = this->nextID_m;
103 : 12 : auto numNodes = this->numNodes_m;
104 [ + - + - ]: 12 : Kokkos::parallel_for(
105 : 12 : "ParticleBase<...>::create(size_t)", policy_type(localNum_m, nLocal),
106 [ + - + - : 12024 : KOKKOS_LAMBDA(const std::int64_t i) { pIDs(i) = nextID + numNodes * i; });
- - ]
107 : : // nextID_m += numNodes_m * (nLocal - localNum_m);
108 : 12 : nextID_m += numNodes_m * nLocal;
109 : 12 : }
110 : :
111 : : // remember that we're creating these new particles
112 : 168 : localNum_m += nLocal;
113 : : }
114 : :
115 : 168 : Comm->allreduce(localNum_m, totalNum_m, 1, std::plus<size_type>());
116 : 168 : }
117 : :
118 : : template <class PLayout, typename... IP>
119 : : void ParticleBase<PLayout, IP...>::createWithID(index_type id) {
120 : : PAssert(layout_m != nullptr);
121 : :
122 : : size_type n = (id > -1) ? 1 : 0;
123 : :
124 : : // temporary change
125 : : index_type tmpNextID = nextID_m;
126 : : nextID_m = id;
127 : : numNodes_m = 0;
128 : :
129 : : create(n);
130 : :
131 : : nextID_m = tmpNextID;
132 : : numNodes_m = Comm->size();
133 : : }
134 : :
135 : : template <class PLayout, typename... IP>
136 : : void ParticleBase<PLayout, IP...>::globalCreate(size_type nTotal) {
137 : : PAssert(layout_m != nullptr);
138 : :
139 : : // Compute the number of particles local to each processor
140 : : size_type nLocal = nTotal / numNodes_m;
141 : :
142 : : const size_t rank = Comm->rank();
143 : :
144 : : size_type rest = nTotal - nLocal * rank;
145 : : if (rank < rest) {
146 : : ++nLocal;
147 : : }
148 : :
149 : : create(nLocal);
150 : : }
151 : :
152 : : template <class PLayout, typename... IP>
153 : : template <typename... Properties>
154 : 12 : void ParticleBase<PLayout, IP...>::destroy(const Kokkos::View<bool*, Properties...>& invalid,
155 : : const size_type destroyNum) {
156 : 12 : this->internalDestroy(invalid, destroyNum);
157 : :
158 : 12 : Comm->allreduce(localNum_m, totalNum_m, 1, std::plus<size_type>());
159 : 12 : }
160 : :
161 : : template <class PLayout, typename... IP>
162 : : template <typename... Properties>
163 : 12 : void ParticleBase<PLayout, IP...>::internalDestroy(
164 : : const Kokkos::View<bool*, Properties...>& invalid, const size_type destroyNum) {
165 [ - + - - ]: 12 : PAssert(destroyNum <= localNum_m);
166 : :
167 : : // If there aren't any particles to delete, do nothing
168 [ - + ]: 12 : if (destroyNum == 0) {
169 : 0 : return;
170 : : }
171 : :
172 : : // If we're deleting all the particles, there's no point in doing
173 : : // anything because the valid region will be empty; we only need to
174 : : // update the particle count
175 [ - + ]: 12 : if (destroyNum == localNum_m) {
176 : 0 : localNum_m = 0;
177 : 0 : return;
178 : : }
179 : :
180 : : using view_type = Kokkos::View<bool*, Properties...>;
181 : : using memory_space = typename view_type::memory_space;
182 : : using execution_space = typename view_type::execution_space;
183 : : using policy_type = Kokkos::RangePolicy<execution_space>;
184 [ + - ]: 12 : auto& locDeleteIndex = deleteIndex_m.get<memory_space>();
185 [ + - ]: 12 : auto& locKeepIndex = keepIndex_m.get<memory_space>();
186 : :
187 : : // Resize buffers, if necessary
188 [ + - ]: 36 : detail::runForAllSpaces([&]<typename MemorySpace>() {
189 [ # # ][ + - : 12 : if (attributes_m.template get<memory_space>().size() > 0) {
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - ]
190 : 12 : int overalloc = Comm->getDefaultOverallocation();
191 : 12 : auto& del = deleteIndex_m.get<memory_space>();
192 : 12 : auto& keep = keepIndex_m.get<memory_space>();
193 [ # # ][ + - : 12 : if (del.size() < destroyNum) {
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - ]
194 : 12 : Kokkos::realloc(del, destroyNum * overalloc);
195 : 12 : Kokkos::realloc(keep, destroyNum * overalloc);
196 : : }
197 : : }
198 : : });
199 : :
200 : : // Reset index buffer
201 [ + - ]: 12 : Kokkos::deep_copy(locDeleteIndex, -1);
202 : :
203 : : // Find the indices of the invalid particles in the valid region
204 [ + - + - ]: 12 : Kokkos::parallel_scan(
205 : 12 : "Scan in ParticleBase::destroy()", policy_type(0, localNum_m - destroyNum),
206 [ + - + - : 12024 : KOKKOS_LAMBDA(const size_t i, int& idx, const bool final) {
+ - - - ]
207 [ # # # # : 12000 : if (final && invalid(i)) {
# # ][ + -
+ + + + +
- + + + +
+ - + + +
+ + - + +
+ + + - +
+ + + + -
+ + + + +
- + + + +
+ - + + +
+ + - + +
+ + + - +
+ + + + -
+ + + + +
- + + +
+ ]
208 : 6000 : locDeleteIndex(idx) = i;
209 : : }
210 [ # # ][ + + : 12000 : if (invalid(i)) {
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + ]
211 : 3000 : idx += 1;
212 : : }
213 : : });
214 [ + - + - ]: 12 : Kokkos::fence();
215 : :
216 : : // Determine the total number of invalid particles in the valid region
217 : 12 : size_type maxDeleteIndex = 0;
218 [ + - + - ]: 12 : Kokkos::parallel_reduce(
219 [ + - ]: 24 : "Reduce in ParticleBase::destroy()", policy_type(0, destroyNum),
220 [ + - ]: 12024 : KOKKOS_LAMBDA(const size_t i, size_t& maxIdx) {
221 [ # # # # : 12000 : if (locDeleteIndex(i) >= 0 && i > maxIdx) {
# # ][ + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + +
+ ]
222 : 2988 : maxIdx = i;
223 : : }
224 : : },
225 [ + - ]: 24 : Kokkos::Max<size_type>(maxDeleteIndex));
226 : :
227 : : // Find the indices of the valid particles in the invalid region
228 [ + - + - ]: 12 : Kokkos::parallel_scan(
229 : : "Second scan in ParticleBase::destroy()",
230 : 12 : Kokkos::RangePolicy<size_type, execution_space>(localNum_m - destroyNum, localNum_m),
231 [ + - + - : 12024 : KOKKOS_LAMBDA(const size_t i, int& idx, const bool final) {
+ - - - ]
232 [ # # # # : 12000 : if (final && !invalid(i)) {
# # ][ + -
+ + + + +
- + + + +
+ - + + +
+ + - + +
+ + + - +
+ + + + -
+ + + + +
- + + + +
+ - + + +
+ + - + +
+ + + - +
+ + + + -
+ + + + +
- + + +
+ ]
233 : 6000 : locKeepIndex(idx) = i;
234 : : }
235 [ # # ][ + + : 12000 : if (!invalid(i)) {
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + ]
236 : 3000 : idx += 1;
237 : : }
238 : : });
239 : :
240 [ + - + - ]: 12 : Kokkos::fence();
241 : :
242 : 12 : localNum_m -= destroyNum;
243 : :
244 : : // We need to delete particles in all memory spaces. If there are any attributes not stored
245 : : // in the memory space we've already been using, we need to copy the index views to the
246 : : // other spaces.
247 : 12 : auto filter = [&]<typename MemorySpace>() {
248 : : return attributes_m.template get<MemorySpace>().size() > 0;
249 : : };
250 [ + - ]: 12 : deleteIndex_m.copyToOtherSpaces<memory_space>(filter);
251 [ + - ]: 12 : keepIndex_m.copyToOtherSpaces<memory_space>(filter);
252 : :
253 : : // Partition the attributes into valid and invalid regions
254 : : // NOTE: The vector elements are pointers, but we want to extract
255 : : // the memory space from the class type, so we explicitly
256 : : // make the lambda argument a pointer to the template parameter
257 [ + - ]: 60 : forAllAttributes([&]<typename Attribute>(Attribute*& attribute) {
258 : : using att_memory_space = typename Attribute::memory_space;
259 : 24 : auto& del = deleteIndex_m.get<att_memory_space>();
260 : 24 : auto& keep = keepIndex_m.get<att_memory_space>();
261 : 24 : attribute->destroy(del, keep, maxDeleteIndex + 1);
262 : : });
263 : : }
264 : :
265 : : template <class PLayout, typename... IP>
266 : : template <typename HashType>
267 : 0 : void ParticleBase<PLayout, IP...>::sendToRank(int rank, int tag,
268 : : std::vector<MPI_Request>& requests,
269 : : const HashType& hash) {
270 [ # # ]: 0 : size_type nSends = hash.size();
271 [ # # ]: 0 : requests.resize(requests.size() + 1);
272 : :
273 [ # # ]: 0 : auto hashes = hash_container_type(hash, [&]<typename MemorySpace>() {
274 : : return attributes_m.template get<MemorySpace>().size() > 0;
275 : : });
276 [ # # ]: 0 : pack(hashes);
277 [ # # ]: 0 : detail::runForAllSpaces([&]<typename MemorySpace>() {
278 [ # # ][ # # : 0 : size_type bufSize = packedSize<MemorySpace>(nSends);
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
279 [ # # ][ # # : 0 : if (bufSize == 0) {
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
280 : 0 : return;
281 : : }
282 : :
283 [ # # ][ # # : 0 : auto buf = Comm->getBuffer<MemorySpace>(bufSize);
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
284 : :
285 [ # # ][ # # : 0 : Comm->isend(rank, tag++, *this, *buf, requests.back(), nSends);
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
286 : 0 : buf->resetWritePos();
287 : 0 : });
288 : 0 : }
289 : :
290 : : template <class PLayout, typename... IP>
291 : 0 : void ParticleBase<PLayout, IP...>::recvFromRank(int rank, int tag, size_type nRecvs) {
292 [ # # ]: 0 : detail::runForAllSpaces([&]<typename MemorySpace>() {
293 [ # # ][ # # : 0 : size_type bufSize = packedSize<MemorySpace>(nRecvs);
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
294 [ # # ][ # # : 0 : if (bufSize == 0) {
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
295 : 0 : return;
296 : : }
297 : :
298 [ # # ][ # # : 0 : auto buf = Comm->getBuffer<MemorySpace>(bufSize);
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
299 : :
300 [ # # ][ # # : 0 : Comm->recv(rank, tag++, *this, *buf, bufSize, nRecvs);
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
301 : 0 : buf->resetReadPos();
302 : 0 : });
303 : 0 : unpack(nRecvs);
304 : 0 : }
305 : :
306 : : template <class PLayout, typename... IP>
307 : : template <typename Archive>
308 : 0 : void ParticleBase<PLayout, IP...>::serialize(Archive& ar, size_type nsends) {
309 : : using memory_space = typename Archive::buffer_type::memory_space;
310 [ # # ]: 0 : forAllAttributes<memory_space>([&]<typename Attribute>(Attribute& att) {
311 : 0 : att->serialize(ar, nsends);
312 : : });
313 : 0 : }
314 : :
315 : : template <class PLayout, typename... IP>
316 : : template <typename Archive>
317 : 0 : void ParticleBase<PLayout, IP...>::deserialize(Archive& ar, size_type nrecvs) {
318 : : using memory_space = typename Archive::buffer_type::memory_space;
319 [ # # ]: 0 : forAllAttributes<memory_space>([&]<typename Attribute>(Attribute& att) {
320 : 0 : att->deserialize(ar, nrecvs);
321 : : });
322 : 0 : }
323 : :
324 : : template <class PLayout, typename... IP>
325 : : template <typename MemorySpace>
326 : 0 : detail::size_type ParticleBase<PLayout, IP...>::packedSize(const size_type count) const {
327 : 0 : size_type total = 0;
328 [ # # ]: 0 : forAllAttributes<MemorySpace>([&]<typename Attribute>(const Attribute& att) {
329 : 0 : total += att->packedSize(count);
330 : : });
331 : 0 : return total;
332 : : }
333 : :
334 : : template <class PLayout, typename... IP>
335 : 0 : void ParticleBase<PLayout, IP...>::pack(const hash_container_type& hash) {
336 [ # # ]: 0 : detail::runForAllSpaces([&]<typename MemorySpace>() {
337 : 0 : auto& att = attributes_m.template get<MemorySpace>();
338 [ # # ][ # # : 0 : for (unsigned j = 0; j < att.size(); j++) {
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
339 : 0 : att[j]->pack(hash.template get<MemorySpace>());
340 : : }
341 : : });
342 : 0 : }
343 : :
344 : : template <class PLayout, typename... IP>
345 : 0 : void ParticleBase<PLayout, IP...>::unpack(size_type nrecvs) {
346 [ # # ]: 0 : detail::runForAllSpaces([&]<typename MemorySpace>() {
347 : 0 : auto& att = attributes_m.template get<MemorySpace>();
348 [ # # ][ # # : 0 : for (unsigned j = 0; j < att.size(); j++) {
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
349 : 0 : att[j]->unpack(nrecvs);
350 : : }
351 : : });
352 : 0 : localNum_m += nrecvs;
353 : 0 : }
354 : : } // namespace ippl
|
