OPAL (Object Oriented Parallel Accelerator Library)  2021.1.99
OPAL
ParticleAttrib.hpp
Go to the documentation of this file.
1 // -*- C++ -*-
2 /***************************************************************************
3  *
4  * The IPPL Framework
5  *
6  * This program was prepared by PSI.
7  * All rights in the program are reserved by PSI.
8  * Neither PSI nor the author(s)
9  * makes any warranty, express or implied, or assumes any liability or
10  * responsibility for the use of this software
11  *
12  * Visit www.amas.web.psi for more details
13  *
14  ***************************************************************************/
15 
16 // -*- C++ -*-
17 /***************************************************************************
18  *
19  * The IPPL Framework
20  *
21  *
22  * Visit http://people.web.psi.ch/adelmann/ for more details
23  *
24  ***************************************************************************/
25 
26 // include files
28 #include "Field/Field.h"
29 #include "Field/LField.h"
30 #include "Message/Message.h"
31 #include "Utility/IpplInfo.h"
32 #include "Utility/PAssert.h"
33 #include "Utility/IpplStats.h"
34 #include "AppTypes/AppTypeTraits.h"
35 
37 // Create a ParticleAttribElem to allow the user to access just the Nth
38 // element of the attribute stored here.
39 template <class T>
43  "No operator()(unsigned) for this element type!!");
45 }
46 
47 
49 // Same as above, but specifying two indices
50 template <class T>
52 ParticleAttrib<T>::operator()(unsigned i, unsigned j) {
54  "No operator()(unsigned,unsigned) for this element type!!");
55  return ParticleAttribElem<T,2U>(*this, vec<unsigned,2U>(i,j));
56 }
57 
58 
60 // Same as above, but specifying three indices
61 template <class T>
63 ParticleAttrib<T>::operator()(unsigned i, unsigned j, unsigned k) {
65  "No operator()(unsigned,unsigned,unsigned) for this element type!!");
66  return ParticleAttribElem<T,3U>(*this, vec<unsigned,3U>(i,j,k));
67 }
68 
69 
71 // Create storage for M particle attributes. The storage is uninitialized.
72 // New items are appended to the end of the array.
73 template<class T>
74 void ParticleAttrib<T>::create(size_t M) {
75 
76  // make sure we have storage for M more items
77  // and push back M items, using default value
78  if (M > 0)
79  {
80  //ParticleList.insert(ParticleList.end(), M, T());
81  ParticleList.insert(ParticleList.begin()+LocalSize, M, T());
82  LocalSize+=M;
83  attributeIsDirty_ = true;
84  }
85 
86 }
87 
88 
90 // Delete the attribute storage for M particle attributes, starting at
91 // the position I. Boolean flag indicates whether to use optimized
92 // destroy method. This function really erases the data, which will
93 // change local indices of the data. The optimized method just copies
94 // data from the end of the storage into the selected block. Otherwise,
95 // we use the std::vector erase method, which preserves data ordering.
96 
97 template<class T>
98 void ParticleAttrib<T>::destroy(size_t M, size_t I, bool optDestroy) {
99 
100 
101  if (M == 0) return;
102  if (optDestroy) {
103  // get iterators for where the data to be deleted begins, and where
104  // the data we copy from the end begins
105  typename ParticleList_t::iterator putloc = ParticleList.begin() + I;
106  typename ParticleList_t::iterator getloc = ParticleList.begin()+LocalSize - M;
107  typename ParticleList_t::iterator endloc = ParticleList.begin()+LocalSize;
108 
109  // make sure we do not copy too much
110  //if ((I + M) > (ParticleList.size() - M))
111  if ((I + M) > (LocalSize - M))
112  getloc = putloc + M;
113 
114  // copy over the data
115  while (getloc != endloc)
116  *putloc++ = *getloc++;
117  // delete the last M items
118  ParticleList.erase(endloc - M, endloc);
119  }
120  else {
121  // just use the erase method
122  typename ParticleList_t::iterator loc = ParticleList.begin() + I;
123  ParticleList.erase(loc, loc + M);
124  }
125  LocalSize-=M;
126  attributeIsDirty_ = true;
127  return;
128 }
129 
130 
132 // Delete the attribute storage for a list of particle destroy events
133 // This really erases the data, which will change local indices
134 // of the data. If we are using the optimized destroy method,
135 // this just copies data from the end of the storage into the selected
136 // block. Otherwise, we use a leading/trailing iterator semantic to
137 // copy data from below and preserve data ordering.
138 
139 template <class T>
140 void ParticleAttrib<T>::destroy(const std::vector< std::pair<size_t,size_t> >& dlist,
141  bool optDestroy)
142 {
143 
144 
145  if (dlist.empty()) return;
146  typedef std::vector< std::pair<size_t,size_t> > dlist_t;
147  if (optDestroy) {
148  // process list in reverse order, since we are backfilling
149  dlist_t::const_reverse_iterator rbeg, rend = dlist.rend();
150  // find point to copy data from
151  typename ParticleList_t::iterator putloc, saveloc;
152  typename ParticleList_t::iterator getloc = ParticleList.begin()+LocalSize;
153  typename ParticleList_t::iterator endloc = ParticleList.begin()+LocalSize;
154  // loop over destroy list and copy data from end of particle list
155  size_t I, M, numParts=0;
156  for (rbeg = dlist.rbegin(); rbeg != rend; ++rbeg) {
157  I = (*rbeg).first; // index number to begin destroy
158  M = (*rbeg).second; // number of particles to destroy
159  numParts += M; // running total of number of particles destroyed
160  // set iterators for data copy
161  putloc = ParticleList.begin() + I;
162  // make sure we do not copy too much
163  if ((I + M) > ((getloc - ParticleList.begin()) - M)) {
164  // we cannot fill all M slots
165  saveloc = getloc; // save endpoint of valid particle data
166  getloc = putloc + M; // move to just past end of section to delete
167  // copy over the data
168  while (getloc != saveloc) {
169  *putloc++ = *getloc++;
170  }
171  // reset getloc for next copy
172  getloc = putloc; // set to end of last copy
173  }
174  else {
175  // fill all M slots using data from end of particle list
176  getloc = getloc - M;
177  saveloc = getloc; // save new endpoint of valid particle data
178  // copy over the data
179  for (size_t m=0; m<M; ++m)
180  *putloc++ = *getloc++;
181  // reset getloc for next copy
182  getloc = saveloc; // set to new endpoint of valid particle data
183  }
184  LocalSize-=M;
185  }
186  // delete storage at end of particle list
187  ParticleList.erase( endloc - numParts, endloc );
188  }
189  else {
190  // just process destroy list using leading/trailing iterators
191  dlist_t::const_iterator dnext = dlist.begin(), dend = dlist.end();
192  size_t putIndex, getIndex, endIndex = LocalSize;
193  putIndex = (*dnext).first; // first index to delete
194  getIndex = putIndex + (*dnext).second; // move past end of destroy event
195  ++dnext; // move to next destroy event
196  // make sure getIndex is not pointing to a deleted particle
197  while (dnext != dend && getIndex == (*dnext).first) {
198  getIndex += (*dnext).second; // move past end of destroy event
199  ++dnext; // move to next destroy event
200  }
201  while (dnext != dend) {
202  // copy into deleted slot
203  ParticleList[putIndex++] = ParticleList[getIndex++];
204  // make sure getIndex points to next non-deleted particle
205  while (dnext != dend && getIndex == (*dnext).first) {
206  getIndex += (*dnext).second; // move past end of destroy event
207  ++dnext; // move to next destroy event
208  }
209  }
210  // one more loop to do any remaining data copying beyond last destroy
211  while (getIndex < endIndex) {
212  // copy into deleted slot
213  ParticleList[putIndex++] = ParticleList[getIndex++];
214  }
215  // now erase any data below last copy
216  typename ParticleList_t::iterator loc = ParticleList.begin() + putIndex;
217  ParticleList.erase(loc, ParticleList.begin()+LocalSize);
218  LocalSize -= ParticleList.begin()+LocalSize - loc;
219  }
220 
221  attributeIsDirty_ = true;
222  return;
223 }
224 
225 
227 // put the data for M particles into a message, starting from index I.
228 // This will either put in local or ghost particle data, but not both.
229 template<class T>
230 size_t
231 ParticleAttrib<T>::putMessage(Message& msg, size_t M, size_t I)
232 {
233 
234  if (M > 0) {
235  if (isTemporary()) {
236  ::putMessage(msg, M);
237  }
238  else {
239  typename ParticleList_t::iterator currp = ParticleList.begin() + I;
240  typename ParticleList_t::iterator endp = currp + M;
241  ::putMessage(msg, currp, endp);
242  }
243  }
244 
245  return M;
246 }
247 
248 
250 // put the data for particles in a list into a Message
251 // This will only work for local particle data right now.
252 template<class T>
253 size_t
255  const std::vector<size_t>& putList)
256 {
257 
258  std::vector<size_t>::size_type M = putList.size();
259 
260  if (M > 0) {
261  if (isTemporary()) {
262  ::putMessage(msg, M);
263  }
264  else {
265  ::putMessage(msg, putList, ParticleList.begin());
266  }
267  }
268 
269  return M;
270 }
271 
272 
274 // Get data out of a Message containing M particle's attribute data,
275 // and store it here. Data is appended to the end of the list. Return
276 // the number of particles retrieved.
277 template<class T>
278 size_t
280 {
281 
282  if (M > 0) {
283  if (isTemporary()) {
284  size_t checksize;
285  ::getMessage(msg, checksize);
286  PAssert_EQ(checksize, M);
287  create(M);
288  }
289  else {
290  size_t currsize = size();
291  create(M);
292  ::getMessage_iter(msg, ParticleList.begin() + currsize);
293  }
294  }
295 
296  return M;
297 }
298 
299 //~ virtual size_t ghostDestroy(size_t, size_t) {
300  //~ return 0;
301  //~ }
302  //~
303  //~ virtual void ghostCreate(size_t)
304  //~ {
305  //~
306  //~ }
307  //~ // puts M particle's data starting from index I into a Message.
308  //~ // Return the number of particles put into the message. This is for
309  //~ // when particles are being swapped to build ghost particle interaction
310  //~ // lists.
311  template<class T>
312  size_t ParticleAttrib<T>::ghostPutMessage(Message&, size_t, size_t) {
313  return 0;
314  }
315  // puts data for a list of particles into a Message, for interaction lists.
316  // Return the number of particles put into the message.
317  template<class T>
318  size_t ParticleAttrib<T>::ghostPutMessage(Message&, const std::vector<size_t>&) {
319  return 0;
320  }
321 //~
322  //~ // Get ghost particle data from a message.
323  //~ virtual size_t ghostGetMessage(Message&, size_t) {
324  //~ return 0;
325  //~ }
326 
327 template<class T>
329 
330  // make sure we have storage for M more items
331  // and push back M items, using default value
332  if (M > 0)
333  {
334  //ParticleList.insert(ParticleList.end(), M, T());
335  ParticleList.insert(ParticleList.end(), M, T());
336  }
337 
338 }
339 
340 template<class T>
341 size_t ParticleAttrib<T>::ghostDestroy(size_t M, size_t I) {
342 
343 
344  if (M > 0)
345  {
346  //ParticleList.insert(ParticleList.end(), M, T());
347  ParticleList.erase(ParticleList.begin() + LocalSize + I, ParticleList.begin() + LocalSize + I + M);
348  }
349  return M;
350 }
351 
352 template<class T>
353 size_t
355 {
356 
357  if (M > 0) {
358  size_t currsize = ParticleList.size();
359  ghostCreate(M);
360  ::getMessage_iter(msg, ParticleList.begin() + currsize);
361  }
362 
363 
364  return M;
365 }
366 
368 // Print out information for debugging purposes. This version just
369 // prints out static information, so it is static
370 template<class T>
372 {
373 
374  o << "PAttr: size = " << ParticleList.size()
375  << ", capacity = " << ParticleList.capacity()
376  << ", temporary = " << isTemporary();
377 }
378 
379 
380 template<class T>
382 {
383  static bool compare(const T &, const T &, bool)
384  {
385  // by default, just return false indicating "no change"
386  return false;
387  }
388 };
389 
390 #define PA_SORT_COMPARE_SCALAR(SCALAR) \
391 template<> \
392 struct PASortCompare<SCALAR> \
393 { \
394  static bool compare(const SCALAR &a, const SCALAR &b, bool ascending) \
395  { \
396  return (ascending ? (a < b) : (a > b)); \
397  } \
398 };
399 
403 PA_SORT_COMPARE_SCALAR(unsigned short)
410 
411 
413 // Calculate a "sort list", which is an array of data of the same
414 // length as this attribute, with each element indicating the
415 // (local) index wherethe ith particle shoulkd go. For example,
416 // if there are four particles, and the sort-list is {3,1,0,2}, that
417 // means the particle currently with index=0 should be moved to the third
418 // position, the one with index=1 should stay where it is, etc.
419 // The optional second argument indicates if the sort should be ascending
420 // (true, the default) or descending (false).
421 template<class T>
422 void ParticleAttrib<T>::calcSortList(SortList_t &slist, bool ascending)
423 {
424  unsigned int i;
425  int j;
426 
427  //Inform dbgmsg("PA<T>::calcSortList");
428 
429  // Resize the sort list, if necessary, to our own size
430  SortList_t::size_type slsize = slist.size();
431  size_t mysize = size();
432  if (slsize < mysize) {
433  // dbgmsg << "Resizing provided sort-list: new size = ";
434  slist.insert(slist.end(), mysize - slsize, (SortList_t::value_type) 0);
435  // dbgmsg << slist.size() << ", attrib size = " << mysize << endl;
436  }
437 
438  // Initialize the sort-list with a negative value, since we check
439  // it later when determing what items to consider in the sort. This
440  // is done to avoid changing the attribute values.
441  for (i=0; i < mysize; ++i)
442  slist[i] = (-1);
443 
444  // OK, this is a VERY simple sort routine, O(N^2): Find min or max
445  // of all elems, store where it goes, then for N-1 elems, etc. I
446  // am sure there is a better way.
447  int firstindx = 0;
448  int lastindx = (mysize - 1);
449  for (i=0; i < mysize; ++i) {
450  int currindx = firstindx;
451  T currval = ParticleList[currindx];
452 
453  for (j=(firstindx + 1); j <= lastindx; ++j) {
454  // skip looking at this item if we already know where it goes
455  if (slist[j] < 0) {
456  // compare current to jth item, if the new one is different
457  // in the right way, save that index
458  if (PASortCompare<T>::compare(ParticleList[j], currval, ascending)) {
459  currindx = j;
460  currval = ParticleList[currindx];
461  }
462  }
463  }
464 
465  // We've found the min or max element, it has index = currindx.
466  // So the currindx's item in the sort-list should say "this will be
467  // the ith item".
468  slist[currindx] = i;
469  // dbgmsg << "Found min/max value " << i << " at position " << currindx;
470  // dbgmsg << " with value = " << currval << endl;
471 
472  // Adjust the min/max index range to look at next time, if necessary
473  while (slist[firstindx] >= 0 && firstindx < lastindx)
474  firstindx++;
475  while (slist[lastindx] >= 0 && firstindx < lastindx)
476  lastindx--;
477  // dbgmsg << " firstindx = " << firstindx << ", lastindx = " << lastindx;
478  // dbgmsg << endl;
479  }
480 }
481 
482 
484 // Process a sort-list, as described for "calcSortList", to reorder
485 // the elements in this attribute. All indices in the sort list are
486 // considered "local", so they should be in the range 0 ... localnum-1.
487 // The sort-list does not have to have been calculated by calcSortList,
488 // it could be calculated by some other means, but it does have to
489 // be in the same format. Note that the routine may need to modify
490 // the sort-list temporarily, but it will return it in the same state.
491 template<class T>
493 {
494  // Make sure the sort-list has the proper length.
495  PAssert_GE(slist.size(), size());
496 
497  // Inform dbgmsg("PA<T>::sort");
498  // dbgmsg << "Sorting " << size() << " items." << endl;
499 
500  // Go through the sort-list instructions, and move items around.
501  int i = 0, j = 0, k = -1, mysize = size();
502  while ( i < mysize ) {
503  PAssert_LT(slist[i], mysize);
504 
505  // skip this swap if the swap-list value is negative. This
506  // happens when we've already put the item in the proper place.
507  if ( i == k || slist[i] < 0 ) {
508  ++i; k = -1;
509  // dbgmsg << "Skipping item " << i << " in sort: slist[" << i << "] = ";
510  // dbgmsg << slist[i] << endl;
511  continue;
512  }
513 
514  j = ( k > 0 ) ? k : slist[i];
515  k = slist[j];
516 
517  // We should not have a negative slist value for the destination
518  PAssert_GE(k, 0);
519 
520  // OK, swap the items
521  std::iter_swap(ParticleList.begin() + i, ParticleList.begin() + j);
522  // dbgmsg << "Swapping item " << i << " to position " << slist[i] << endl;
523 
524 
525  // then indicate that we've put this
526  // item in the proper location.
527  slist[j] -= mysize;
528  }
529 
530 
531  // Restore the sort-list
532  for (i=0; i < mysize; ++i) {
533  if (slist[i] < 0)
534  slist[i] += mysize;
535  // dbgmsg << "At end of sort: restored slist[" << i << "] = " << slist[i];
536  // dbgmsg << ", data[" << i << "] = " << ParticleList[i] << endl;
537  }
538 }
539 
541 // scatter functions
543 //mwerks Moved into class definition (.h file).
544 
545 
547 // gather functions
549 //mwerks Moved into class definition (.h file).
550 
552 // Global function templates
554 
555 
557 // scatter functions for number density
559 
560 template <class FT, unsigned Dim, class M, class C, class PT, class IntOp>
561 void
563  const IntOp& /*intop*/, FT val) {
564 
565  // make sure field is uncompressed and guard cells are zeroed
566  f.Uncompress();
567  FT zero = 0;
568  f.setGuardCells(zero);
569 
570  const M& mesh = f.get_mesh();
571  // iterate through particles and call scatter operation
572  typename ParticleAttrib< Vektor<PT,Dim> >::const_iterator ppiter;
573  size_t i = 0;
574  for (ppiter = pp.cbegin(); i < pp.size(); ++i,++ppiter)
575  IntOp::scatter(val,f,*ppiter,mesh);
576 
577  // accumulate values in guard cells
578  f.accumGuardCells();
579 
580  INCIPPLSTAT(incParticleScatters);
581  return;
582 }
583 
584 template <class FT, unsigned Dim, class M, class C, class PT,
585  class IntOp, class CacheData>
586 void
588  const IntOp& /*intop*/, ParticleAttrib<CacheData>& cache, FT val) {
589 
590  // make sure field is uncompressed and guard cells are zeroed
591  f.Uncompress();
592  FT zero = 0;
593  f.setGuardCells(zero);
594 
595  const M& mesh = f.get_mesh();
596  // iterate through particles and call scatter operation
598  typename ParticleAttrib<CacheData>::iterator citer=cache.begin();
599  size_t i = 0;
600  for (ppiter = pp.begin(); i < pp.size(); ++i, ++ppiter, ++citer)
601  IntOp::scatter(val,f,*ppiter,mesh,*citer);
602 
603  // accumulate values in guard cells
604  f.accumGuardCells();
605 
606  INCIPPLSTAT(incParticleScatters);
607  return;
608 }
609 
610 template <class FT, unsigned Dim, class M, class C,
611  class IntOp, class CacheData>
612 void
613 scatter(Field<FT,Dim,M,C>& f, const IntOp& /*intop*/,
614  const ParticleAttrib<CacheData>& cache, FT val) {
615 
616  // make sure field is uncompressed and guard cells are zeroed
617  f.Uncompress();
618  FT zero = 0;
619  f.setGuardCells(zero);
620 
621  // iterate through particles and call scatter operation
622  typename ParticleAttrib<CacheData>::iterator citer, cend=cache.begin()+cache.size();//not sure jp
623  for (citer = cache.begin(); citer != cend; ++citer)
624  IntOp::scatter(val,f,*citer);
625 
626  // accumulate values in guard cells
627  f.accumGuardCells();
628 
629  INCIPPLSTAT(incParticleScatters);
630  return;
631 }
632 
633 
634 /***************************************************************************
635  * $RCSfile: ParticleAttrib.cpp,v $ $Author: adelmann $
636  * $Revision: 1.1.1.1 $ $Date: 2003/01/23 07:40:28 $
637  * IPPL_VERSION_ID: $Id: ParticleAttrib.cpp,v 1.1.1.1 2003/01/23 07:40:28 adelmann Exp $
638  ***************************************************************************/
T * value_type(const SliceIterator< T > &)
const unsigned Dim
void putMessage(Message &m, const T &t)
Definition: Message.h:549
void getMessage_iter(Message &m, OutputIterator o)
Definition: Message.h:595
void getMessage(Message &m, T &t)
Definition: Message.h:572
void scatter(Field< FT, Dim, M, C > &f, const ParticleAttrib< Vektor< PT, Dim > > &pp, const IntOp &, FT val)
#define PA_SORT_COMPARE_SCALAR(SCALAR)
#define INCIPPLSTAT(stat)
Definition: IpplStats.h:236
#define PAssert_LT(a, b)
Definition: PAssert.h:106
#define PInsist(c, m)
Definition: PAssert.h:120
#define PAssert_EQ(a, b)
Definition: PAssert.h:104
#define PAssert_GE(a, b)
Definition: PAssert.h:109
std::string::iterator iterator
Definition: MSLang.h:16
Definition: Vektor.h:32
Definition: Field.h:33
Mesh_t & get_mesh() const
Definition: Field.h:110
ParticleAttribElem< T, 1U > operator()(unsigned)
virtual void printDebug(Inform &)
size_t size(void) const
ParticleAttribBase::SortList_t SortList_t
virtual void ghostCreate(size_t)
virtual void sort(SortList_t &slist)
virtual size_t ghostGetMessage(Message &, size_t)
iterator begin()
virtual size_t putMessage(Message &, size_t, size_t)
virtual void destroy(size_t M, size_t I, bool optDestroy=true)
virtual size_t getMessage(Message &, size_t)
const_iterator cbegin() const
virtual size_t ghostDestroy(size_t, size_t)
virtual size_t ghostPutMessage(Message &, size_t, size_t)
virtual void create(size_t)
void accumGuardCells()
Definition: BareField.hpp:698
void Uncompress() const
Definition: BareField.hpp:1005
void setGuardCells(const T &) const
Definition: BareField.hpp:616
static bool compare(const T &, const T &, bool)
Definition: Inform.h:42
Definition: Vec.h:22