src/Particle/ParticleAttrib.cpp

Go to the documentation of this file.

// -*- C++ -*-
/***************************************************************************
 *
 * The IPPL Framework
 * 
 * This program was prepared by PSI. 
 * All rights in the program are reserved by PSI.
 * Neither PSI nor the author(s)
 * makes any warranty, express or implied, or assumes any liability or
 * responsibility for the use of this software
 *
 * Visit http://www.acl.lanl.gov/POOMS for more details
 *
 ***************************************************************************/

// -*- C++ -*-
/***************************************************************************
 *
 * The IPPL Framework
 * 
 *
 * Visit http://people.web.psi.ch/adelmann/ for more details
 *
 ***************************************************************************/

// include files
#include "Particle/ParticleAttrib.h"
#include "Field/Field.h"
#include "Field/LField.h"
#include "Message/Message.h"
#include "Utility/IpplInfo.h"
#include "Utility/PAssert.h"
#include "Utility/IpplStats.h"
#include "AppTypes/AppTypeTraits.h"
#include "Profile/Profiler.h"


// Create a ParticleAttribElem to allow the user to access just the Nth
// element of the attribute stored here.
template <class T>
ParticleAttribElem<T,1U>
ParticleAttrib<T>::operator()(unsigned i) {
  PInsist(AppTypeTraits<T>::ElemDim > 0,
          "No operator()(unsigned) for this element type!!");
  return ParticleAttribElem<T,1U>(*this, vec<unsigned,1U>(i));
}


// Same as above, but specifying two indices
template <class T>
ParticleAttribElem<T,2U>
ParticleAttrib<T>::operator()(unsigned i, unsigned j) {
  PInsist(AppTypeTraits<T>::ElemDim > 1,
          "No operator()(unsigned,unsigned) for this element type!!");
  return ParticleAttribElem<T,2U>(*this, vec<unsigned,2U>(i,j));
}


// Same as above, but specifying three indices
template <class T>
ParticleAttribElem<T,3U>
ParticleAttrib<T>::operator()(unsigned i, unsigned j, unsigned k) {
  PInsist(AppTypeTraits<T>::ElemDim > 2,
          "No operator()(unsigned,unsigned,unsigned) for this element type!!");
  return ParticleAttribElem<T,3U>(*this, vec<unsigned,3U>(i,j,k));
}


// Create storage for M particle attributes.  The storage is uninitialized.
// New items are appended to the end of the array.
template<class T>
void ParticleAttrib<T>::create(size_t M) {
  TAU_PROFILE("ParticleAttrib::create()", "void (unsigned)", TAU_PARTICLE);
  // make sure we have storage for M more items
  // and push back M items, using default value
  if (M > 0)
    ParticleList.insert(ParticleList.end(), M, T());
}


// Delete the attribute storage for M particle attributes, starting at
// the position I.  Boolean flag indicates whether to use optimized 
// destroy method.  This function really erases the data, which will
// change local indices of the data.  The optimized method just copies
// data from the end of the storage into the selected block.  Otherwise,
// we use the std::vector erase method, which preserves data ordering.

template<class T>
void ParticleAttrib<T>::destroy(size_t M, size_t I, bool optDestroy) {
  TAU_PROFILE("ParticleAttrib::destroy()", "void (unsigned, unsigned, bool)", TAU_PARTICLE);

  if (M == 0) return;
  if (optDestroy) {
    // get iterators for where the data to be deleted begins, and where
    // the data we copy from the end begins
    typename ParticleList_t::iterator putloc = ParticleList.begin() + I;
    typename ParticleList_t::iterator getloc = ParticleList.end() - M;
    typename ParticleList_t::iterator endloc = ParticleList.end();

    // make sure we do not copy too much
    if ((I + M) > (ParticleList.size() - M))
      getloc = putloc + M;

    // copy over the data
    while (getloc != endloc)
      *putloc++ = *getloc++;
    // delete the last M items
    ParticleList.erase(endloc - M, endloc);
  }
  else {
    // just use the erase method
    typename ParticleList_t::iterator loc = ParticleList.begin() + I;
    ParticleList.erase(loc, loc + M);
  }
  return;
}


// Delete the attribute storage for a list of particle destroy events
// This really erases the data, which will change local indices
// of the data.  If we are using the optimized destroy method,
// this just copies data from the end of the storage into the selected
// block.  Otherwise, we use a leading/trailing iterator semantic to
// copy data from below and  preserve data ordering.

template <class T>
void ParticleAttrib<T>::destroy(const vector< pair<size_t,size_t> >& dlist,
                                bool optDestroy)
{
  TAU_PROFILE("ParticleAttrib::destroy()", "void (vector< pair<unsigned,unsigned> >, bool)", TAU_PARTICLE); 

  if (dlist.empty()) return;
  typedef vector< pair<size_t,size_t> > dlist_t;
  if (optDestroy) {
    // process list in reverse order, since we are backfilling
    dlist_t::const_reverse_iterator rbeg, rend = dlist.rend();
    // find point to copy data from
    typename ParticleList_t::iterator putloc, saveloc;
    typename ParticleList_t::iterator getloc = ParticleList.end();
    typename ParticleList_t::iterator endloc = ParticleList.end();
    // loop over destroy list and copy data from end of particle list
    size_t I, M, numParts=0;
    for (rbeg = dlist.rbegin(); rbeg != rend; ++rbeg) {
      I = (*rbeg).first;   // index number to begin destroy
      M = (*rbeg).second;  // number of particles to destroy
      numParts += M;       // running total of number of particles destroyed
      // set iterators for data copy
      putloc = ParticleList.begin() + I;
      // make sure we do not copy too much
      if ((I + M) > ((getloc - ParticleList.begin()) - M)) {
        // we cannot fill all M slots
        saveloc = getloc;  // save endpoint of valid particle data
        getloc = putloc + M;  // move to just past end of section to delete
        // copy over the data
        while (getloc != saveloc) {
          *putloc++ = *getloc++;
        }
        // reset getloc for next copy
        getloc = putloc;  // set to end of last copy
      }
      else {
        // fill all M slots using data from end of particle list
        getloc = getloc - M;
        saveloc = getloc;  // save new endpoint of valid particle data
        // copy over the data
        for (size_t m=0; m<M; ++m)
          *putloc++ = *getloc++;
        // reset getloc for next copy
        getloc = saveloc;  // set to new endpoint of valid particle data
      }
    }
    // delete storage at end of particle list
    ParticleList.erase( endloc - numParts, endloc );
  }
  else {
    // just process destroy list using leading/trailing iterators
    dlist_t::const_iterator dnext = dlist.begin(), dend = dlist.end();
    size_t putIndex, getIndex, endIndex = ParticleList.size();
    putIndex = (*dnext).first;  // first index to delete
    getIndex = putIndex + (*dnext).second;  // move past end of destroy event
    ++dnext;  // move to next destroy event
    // make sure getIndex is not pointing to a deleted particle
    while (dnext != dend && getIndex == (*dnext).first) {
      getIndex += (*dnext).second;  // move past end of destroy event
      ++dnext;                      // move to next destroy event
    }
    while (dnext != dend) {
      // copy into deleted slot
      ParticleList[putIndex++] = ParticleList[getIndex++];
      // make sure getIndex points to next non-deleted particle
      while (dnext != dend && getIndex == (*dnext).first) {
        getIndex += (*dnext).second;  // move past end of destroy event
        ++dnext;                      // move to next destroy event
      }
    }
    // one more loop to do any remaining data copying beyond last destroy
    while (getIndex < endIndex) {
      // copy into deleted slot
      ParticleList[putIndex++] = ParticleList[getIndex++];
    }
    // now erase any data below last copy
    typename ParticleList_t::iterator loc = ParticleList.begin() + putIndex;
    ParticleList.erase(loc, ParticleList.end());
  }
  return;
}


// put the data for M particles into a message, starting from index I.
// This will either put in local or ghost particle data, but not both.
template<class T>
size_t
ParticleAttrib<T>::putMessage(Message& msg, size_t M, size_t I)
{
  TAU_PROFILE("ParticleAttrib::putMessage()", 
              "unsigned (Message, unsigned, unsigned)", TAU_PARTICLE);

  if (M > 0) {
    if (isTemporary()) {
      ::putMessage(msg, M);
    }
    else {
      typename ParticleList_t::iterator currp = ParticleList.begin() + I;
      typename ParticleList_t::iterator endp = currp + M;
      ::putMessage(msg, currp, endp);
    }
  }

  return M;
}


// put the data for particles in a list into a Message
// This will only work for local particle data right now.
template<class T>
size_t
ParticleAttrib<T>::putMessage(Message& msg,
                              const vector<size_t>& putList)
{  
  TAU_PROFILE("ParticleAttrib::putMessage()", 
    "unsigned (Message, vector<unsigned>)", TAU_PARTICLE); 

  vector<size_t>::size_type M = putList.size();

  if (M > 0) {
    if (isTemporary()) {
      ::putMessage(msg, M);
    }
    else {
      ::putMessage(msg, putList, ParticleList.begin());
    }
  }

  return M;
}


// Get data out of a Message containing M particle's attribute data,
// and store it here.  Data is appended to the end of the list.  Return
// the number of particles retrieved.
template<class T>
size_t
ParticleAttrib<T>::getMessage(Message& msg, size_t M)
{
  TAU_PROFILE("ParticleAttrib::getMessage()", 
    "unsigned (Message, unsigned)", TAU_PARTICLE); 

  if (M > 0) {
    if (isTemporary()) {
      size_t checksize;
      ::getMessage(msg, checksize);
      PAssert(checksize == M);
      create(M);
    }
    else {
      size_t currsize = size();
      create(M);
      ::getMessage_iter(msg, ParticleList.begin() + currsize);
    }
  }

  return M;
}


// Print out information for debugging purposes.  This version just
// prints out static information, so it is static
template<class T>
void ParticleAttrib<T>::printDebug(Inform& o)
{
  TAU_PROFILE("ParticleAttrib::printDebug()", 
    "void (Inform)", TAU_PARTICLE | TAU_IO); 

  o << "PAttr: size = " << ParticleList.size()
    << ", capacity = " << ParticleList.capacity()
    << ", temporary = " << isTemporary();
}


template<class T>
struct PASortCompare
{
  static bool compare(const T &, const T &, bool)
  {
    // by default, just return false indicating "no change"
    return false;
  }
};

#define PA_SORT_COMPARE_SCALAR(SCALAR)                                  \
template<>                                                              \
struct PASortCompare<SCALAR>                                            \
{                                                                       \
  static bool compare(const SCALAR &a, const SCALAR &b, bool ascending) \
  {                                                                     \
    return (ascending ? (a < b) : (a > b));                             \
  }                                                                     \
};

PA_SORT_COMPARE_SCALAR(char)
PA_SORT_COMPARE_SCALAR(unsigned char)
PA_SORT_COMPARE_SCALAR(short)
PA_SORT_COMPARE_SCALAR(unsigned short)
PA_SORT_COMPARE_SCALAR(int)
PA_SORT_COMPARE_SCALAR(unsigned int)
PA_SORT_COMPARE_SCALAR(long)
PA_SORT_COMPARE_SCALAR(unsigned long)
PA_SORT_COMPARE_SCALAR(float)
PA_SORT_COMPARE_SCALAR(double)


// Calculate a "sort list", which is an array of data of the same
// length as this attribute, with each element indicating the
// (local) index wherethe ith particle shoulkd go.  For example, 
// if there are four particles, and the sort-list is {3,1,0,2}, that
// means the particle currently with index=0 should be moved to the third
// position, the one with index=1 should stay where it is, etc.
// The optional second argument indicates if the sort should be ascending
// (true, the default) or descending (false).
template<class T>
void ParticleAttrib<T>::calcSortList(SortList_t &slist, bool ascending)
{
  int i, j;

  //Inform dbgmsg("PA<T>::calcSortList");

  // Resize the sort list, if necessary, to our own size
  SortList_t::size_type slsize = slist.size();
  size_t mysize = size();
  if (slsize < mysize) {
    // dbgmsg << "Resizing provided sort-list: new size = ";
    slist.insert(slist.end(), mysize - slsize, (SortList_t::value_type) 0);
    // dbgmsg << slist.size() << ", attrib size = " << mysize << endl;
  }

  // Initialize the sort-list with a negative value, since we check
  // it later when determing what items to consider in the sort.  This
  // is done to avoid changing the attribute values.
  for (i=0; i < mysize; ++i)
    slist[i] = (-1);

  // OK, this is a VERY simple sort routine, O(N^2): Find min or max
  // of all elems, store where it goes, then for N-1 elems, etc.  I
  // am sure there is a better way.
  int firstindx = 0;
  int lastindx = (mysize - 1);
  for (i=0; i < mysize; ++i) {
    int currindx = firstindx;
    T currval = ParticleList[currindx];

    for (j=(firstindx + 1); j <= lastindx; ++j) {
      // skip looking at this item if we already know where it goes
      if (slist[j] < 0) {
        // compare current to jth item, if the new one is different
        // in the right way, save that index
        if (PASortCompare<T>::compare(ParticleList[j], currval, ascending)) {
          currindx = j;
          currval = ParticleList[currindx];
        }
      }
    }

    // We've found the min or max element, it has index = currindx.
    // So the currindx's item in the sort-list should say "this will be
    // the ith item".
    slist[currindx] = i;
    // dbgmsg << "Found min/max value " << i << " at position " << currindx;
    // dbgmsg << " with value = " << currval << endl;

    // Adjust the min/max index range to look at next time, if necessary
    while (slist[firstindx] >= 0 && firstindx < lastindx)
      firstindx++;
    while (slist[lastindx] >= 0 && firstindx < lastindx)
      lastindx--;
    // dbgmsg << " firstindx = " << firstindx << ", lastindx = " << lastindx;
    // dbgmsg << endl;
  }
}


// Process a sort-list, as described for "calcSortList", to reorder
// the elements in this attribute.  All indices in the sort list are
// considered "local", so they should be in the range 0 ... localnum-1.
// The sort-list does not have to have been calcualted by calcSortList,
// it could be calculated by some other means, but it does have to
// be in the same format.  Note that the routine may need to modify
// the sort-list temporarily, but it will return it in the same state.
template<class T>
void ParticleAttrib<T>::sort(SortList_t &slist)
{
  // Make sure the sort-list has the proper length.
  PAssert(slist.size() >= size());

  // Inform dbgmsg("PA<T>::sort");
  // dbgmsg << "Sorting " << size() << " items." << endl;

  // Go through the sort-list instructions, and move items around.
  int i, mysize = size();
  for (i=0; i < mysize; ++i) {
    PAssert(slist[i] < mysize);

    // skip this swap if the swap-list value is negative.  This
    // happens when we've already put the item in the proper place.
    if (slist[i] >= 0 && slist[i] != i) {
      // We should not have a negative slist value for the destination
      PAssert(slist[slist[i]] >= 0);

      // OK, swap the items
      T tmpval = ParticleList[i];
      ParticleList[i] = ParticleList[slist[i]];
      ParticleList[slist[i]] = tmpval;
      //dbgmsg << "Swapping item " << i << " to position " << slist[i] << endl;

      // then indicate that we've put this
      // item in the proper location.
      slist[slist[i]] -= mysize;
      // } else {
      //  dbgmsg << "Skipping item " << i << " in sort: slist[" << i << "] = ";
      //  dbgmsg << slist[i] << endl;
    }
  }

  // Restore the sort-list
  for (i=0; i < mysize; ++i) {
    if (slist[i] < 0)
      slist[i] += mysize;
    // dbgmsg << "At end of sort: restored slist[" << i << "] = " << slist[i];
    // dbgmsg << ", data[" << i << "] = " << ParticleList[i] << endl;
  }
}


// scatter functions
//mwerks Moved into class definition (.h file).


// gather functions
//mwerks Moved into class definition (.h file).

// Global function templates


// scatter functions for number density

template <class FT, unsigned Dim, class M, class C, class PT, class IntOp>
void
scatter(Field<FT,Dim,M,C>& f, const ParticleAttrib< Vektor<PT,Dim> >& pp,
        const IntOp& intop, FT val) {
  TAU_TYPE_STRING(taustr, "void (" + CT(f) + ", " + CT(pp) + ", " 
    + CT(intop)  + ", " + CT(val)  + " )" );
  TAU_PROFILE("scatter()", taustr, TAU_PARTICLE);

  // make sure field is uncompressed and guard cells are zeroed
  f.Uncompress();
  FT zero = 0;
  f.setGuardCells(zero);

  const M& mesh = f.get_mesh();
  // iterate through particles and call scatter operation
  typename ParticleAttrib< Vektor<PT,Dim> >::iterator ppiter, ppend=pp.end();
  for (ppiter = pp.begin(); ppiter != ppend; ++ppiter)
    IntOp::scatter(val,f,*ppiter,mesh);

  // accumulate values in guard cells
  f.accumGuardCells();

  INCIPPLSTAT(incParticleScatters);
  return;
}

template <class FT, unsigned Dim, class M, class C, class PT,
          class IntOp, class CacheData>
void
scatter(Field<FT,Dim,M,C>& f, const ParticleAttrib< Vektor<PT,Dim> >& pp,
  const IntOp& intop, ParticleAttrib<CacheData>& cache, FT val) {

  TAU_TYPE_STRING(taustr, "void (" + CT(f) + ", " + CT(pp) + ", " 
    + CT(intop)  + ", " + CT(cache)  + ", " + CT(val)  + " )" );
  TAU_PROFILE("scatter()", taustr, TAU_PARTICLE);

  // make sure field is uncompressed and guard cells are zeroed
  f.Uncompress();
  FT zero = 0;
  f.setGuardCells(zero);

  const M& mesh = f.get_mesh();
  // iterate through particles and call scatter operation
  typename ParticleAttrib< Vektor<PT,Dim> >::iterator ppiter, ppend=pp.end();
  typename ParticleAttrib<CacheData>::iterator citer=cache.begin();
  for (ppiter = pp.begin(); ppiter != ppend; ++ppiter, ++citer)
    IntOp::scatter(val,f,*ppiter,mesh,*citer);

  // accumulate values in guard cells
  f.accumGuardCells();

  INCIPPLSTAT(incParticleScatters);
  return;
}

template <class FT, unsigned Dim, class M, class C,
          class CacheData, class IntOp>
void
scatter(Field<FT,Dim,M,C>& f, const IntOp& intop,
        const ParticleAttrib<CacheData>& cache, FT val) {

  TAU_TYPE_STRING(taustr, "void (" + CT(f) + ", " + CT(intop) + ", " 
    + CT(cache)  + ", " + CT(val)  + " )" );
  TAU_PROFILE("scatter()", taustr, TAU_PARTICLE);

  // make sure field is uncompressed and guard cells are zeroed
  f.Uncompress();
  FT zero = 0;
  f.setGuardCells(zero);

  // iterate through particles and call scatter operation
  typename ParticleAttrib<CacheData>::iterator citer, cend=cache.end();
  for (citer = cache.begin(); citer != cend; ++citer)
    IntOp::scatter(val,f,*citer);

  // accumulate values in guard cells
  f.accumGuardCells();

  INCIPPLSTAT(incParticleScatters);
  return;
}


/***************************************************************************
 * $RCSfile: ParticleAttrib.cpp,v $   $Author: adelmann $
 * $Revision: 1.1.1.1 $   $Date: 2003/01/23 07:40:28 $
 * IPPL_VERSION_ID: $Id: ParticleAttrib.cpp,v 1.1.1.1 2003/01/23 07:40:28 adelmann Exp $ 
 ***************************************************************************/

---