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
|
