src/Particle/ParticleBase.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/ParticleBase.h"
#include "Particle/ParticleLayout.h"
#include "Particle/ParticleAttrib.h"
#include "Message/Message.h"
#include "Message/Communicate.h"
#include "Utility/Inform.h"
#include "Utility/PAssert.h"
#include "Utility/IpplInfo.h"
#include "Utility/IpplStats.h"
#include "Profile/Profiler.h"

#ifdef IPPL_STDSTL
#include <algorithm>
using namespace std;
#else
#include <algo.h>
#endif // IPPL_STDSTL


// For a ParticleBase that was created with the default constructor,
// initialize performs the same actions as are done in the non-default
// constructor.  If this object has already been initialized, it is
// an error.  For initialize, you must supply a layout instance.
template<class PLayout>
void ParticleBase<PLayout>::initialize(PLayout *layout) {
  TAU_TYPE_STRING(taustr, CT(*this) + " void (" + CT(*layout) + ")"); 
  TAU_PROFILE("ParticleBase::initialize()", taustr, TAU_PARTICLE);

  // make sure we have not already been initialized, and that we have
  // been given a good layout
  PAssert(Layout == 0);
  PAssert(layout != 0);

  // save the layout, and perform setup tasks
  Layout = layout;
  setup();
}


// set up this new object:  add attributes and check in to the layout
template<class PLayout>
void ParticleBase<PLayout>::setup() {
  TAU_TYPE_STRING(taustr, CT(*this) + " void ()"); 
  TAU_PROFILE("ParticleBase::setup()", taustr, TAU_PARTICLE);

  TotalNum = 0;
  LocalNum = 0;
  DestroyNum = 0;
  GhostNum = 0;

  // shift ID back so that first ID retrieved is myNode
  NextID = Ippl::Comm->myNode() - Ippl::Comm->getNodes();

  // Create position R and index ID attributes for the particles,
  // and add them to our list of attributes.
  addAttribute(R);
  addAttribute(ID);

  // indicate we have created a new ParticleBase object
  INCIPPLSTAT(incParticleBases);
}


// Return a boolean value indicating if we are on a processor which can
// be used for single-node particle creation and initialization
template<class PLayout>
bool ParticleBase<PLayout>::singleInitNode() const {
  TAU_TYPE_STRING(taustr, CT(*this) + " void ()"); 
  TAU_PROFILE("ParticleBase::singleInitNode()", taustr, TAU_PARTICLE);
  
  return (Ippl::Comm->myNode() == 0);
}


// Return a new unique ID value for use by new particles.
// The ID number = (i * numprocs) + myproc, i = 0, 1, 2, ...
template<class PLayout>
unsigned ParticleBase<PLayout>::getNextID() {
  TAU_TYPE_STRING(taustr, CT(*this) + " void ()"); 
  TAU_PROFILE("ParticleBase::getNextID()", taustr, TAU_PARTICLE);

  return (NextID += Ippl::Comm->getNodes());
}

// Reset the particle ID's to be globally consecutive, 0 thru TotalNum.
template <class PLayout>
void ParticleBase<PLayout>::resetID(void) {
  TAU_TYPE_STRING(taustr, CT(*this) + " void ()");
  TAU_PROFILE("ParticleBase::resetID()", taustr, TAU_PARTICLE);

  unsigned int nodes = Ippl::getNodes();
  unsigned int myNode = Ippl::myNode();
  size_t localNum = this->getLocalNum();
  size_t ip;
  int master = 0;
  if (myNode == master) {
    // Node 0 can immediately set its ID's to 0 thru LocalNum-1
    for (ip=0; ip<localNum; ++ip) 
      this->ID[ip] = ip;
    // if there is only one processor, we are done
    if (nodes == 1) return;
    // parallel case: must find out how many particles each processor has
    // Node 0 gathers this information into an array
    size_t *lp;
    lp = new size_t[nodes];
    lp[0] = localNum;  // enter our own number of particles
    // get next message tag and receive messages
    int tag1 = Ippl::Comm->next_tag(P_RESET_ID_TAG,P_LAYOUT_CYCLE);
    Message* msg1;
    for (ip=1; ip<nodes; ++ip) {
      int rnode = COMM_ANY_NODE;
      msg1 = Ippl::Comm->receive_block(rnode,tag1);
      PAssert(msg1);
      msg1->get(lp[rnode]);
      delete msg1;
    }
    // now we should have all the localnum values.
    // figure out starting ID for each processor and send back
    size_t current, sum = 0;
    for (ip=0; ip<nodes; ++ip) {
      current = lp[ip];
      lp[ip] = sum;
      sum += current;
    }
    // send initial ID values back out
    int tag2 = Ippl::Comm->next_tag(P_RESET_ID_TAG,P_LAYOUT_CYCLE);
    for (ip=1; ip<nodes; ++ip) {
      Message* msg2 = new Message;
      msg2->put(lp[ip]);
      bool success = Ippl::Comm->send(msg2,ip,tag2);
      PAssert(success);
    }
    // we are done
    return;
  }
  else {
    // first send number of local particles to Node 0
    int tag1 = Ippl::Comm->next_tag(P_RESET_ID_TAG,P_LAYOUT_CYCLE);
    Message* msg1 = new Message;
    msg1->put(localNum);
    bool success = Ippl::Comm->send(msg1,master,tag1);
    PAssert(success);
    // now receive back our initial ID number
    size_t initialID;
    int tag2 = Ippl::Comm->next_tag(P_RESET_ID_TAG,P_LAYOUT_CYCLE);
    Message* msg2 = Ippl::Comm->receive_block(master,tag2);
    PAssert(msg2);
    msg2->get(initialID);
    delete msg2;
    // now reset our particle ID's using this initial value
    for (ip=0; ip<localNum; ++ip) 
      this->ID[ip] = ip + initialID;
    // we are done
    return;
  }
}


// put the data for M particles starting from local index I in a Message
template<class PLayout>
size_t
ParticleBase<PLayout>::putMessage(Message& msg, size_t M, size_t I) {
  TAU_PROFILE("ParticleBase::putMessage()", 
    "unsigned (Message, unsigned, unsigned)", TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  // put into message the number of items in the message
  msg.put(M);

  // go through all the attributes and put their data in the message
  if (M > 0) {
    // this routine should only be called for local particles; call
    // ghostPutMessage to put in particles which might be ghost particles
    PAssert(I < R.size() && (I+M) <= R.size());

    attrib_container_t::iterator abeg = AttribList.begin();
    attrib_container_t::iterator aend = AttribList.end();
    for ( ; abeg != aend; abeg++ )
      (*abeg)->putMessage(msg, M, I);
  }

  return M;
}


// put the data for a list of particles in a Message
template<class PLayout>
size_t
ParticleBase<PLayout>::putMessage(Message& msg,
                                  const vector<size_t>& putList)
{
  TAU_PROFILE("ParticleBase::putMessage()", 
    "unsigned (Message, vector<unsigned>)", TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

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

  // go through all the attributes and put their data in the message
  if (M > 0) {
    attrib_container_t::iterator abeg = AttribList.begin();
    attrib_container_t::iterator aend = AttribList.end();
    for ( ; abeg != aend; ++abeg )
      (*abeg)->putMessage(msg, putList);
  }

  return M;
}


// retrieve particles from the given message and store them
template<class PLayout>
size_t ParticleBase<PLayout>::getMessage(Message& msg) {
  TAU_PROFILE("ParticleBase::getMessage()", "unsigned (Message)",
              TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  // get the number of items in the message
  size_t numitems;
  msg.get(numitems);

  // go through all the attributes and get their data from the message
  if (numitems > 0) {
    attrib_container_t::iterator abeg = AttribList.begin();
    attrib_container_t::iterator aend = AttribList.end();
    for ( ; abeg != aend; abeg++ )
      (*abeg)->getMessage(msg, numitems);
  }

  return numitems;
}


// retrieve particles from the given message and store them, also 
// signaling we are creating the given number of particles.  Return the
// number of particles created.
template<class PLayout>
size_t ParticleBase<PLayout>::getMessageAndCreate(Message& msg) {
  TAU_PROFILE("ParticleBase::getMessageAndCreate()", "unsigned (Message)",
              TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  // call the regular create, and add in the particles to our LocalNum
  size_t numcreate = getMessage(msg);
  LocalNum += numcreate;
  ADDIPPLSTAT(incParticlesCreated,numcreate);
  return numcreate;
}


// create M new particles on this processor
template<class PLayout>
void ParticleBase<PLayout>::create(size_t M) {
  TAU_PROFILE("ParticleBase::create()", "void (unsigned)", TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  // go through all the attributes, and allocate space for M new particles
  attrib_container_t::iterator abeg = AttribList.begin();
  attrib_container_t::iterator aend = AttribList.end();
  for ( ; abeg != aend; abeg++ )
    (*abeg)->create(M);

  // set the unique ID value for these new particles
  size_t i1 = LocalNum;
  size_t i2 = i1 + M;
  while (i1 < i2)
    ID[i1++] = getNextID();

  // remember that we're creating these new particles
  LocalNum += M;
  ADDIPPLSTAT(incParticlesCreated,M);
}


// create np new particles globally, equally distributed among all processors
template<class PLayout>
void ParticleBase<PLayout>::globalCreate(size_t np) {
  TAU_PROFILE("ParticleBase::globalCreate()", "void (unsigned)", TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  // Compute the number of particles local to each processor:
  unsigned nPE = IpplInfo::getNodes();
  unsigned npLocal = np/nPE;

  // Handle the case where np/nPE is not an integer:
  unsigned myPE = IpplInfo::myNode();
  unsigned rem = np - npLocal * nPE;
  if (myPE < rem) ++npLocal;

  // Now each PE calls the local ParticleBase::create() function to create it's
  // local number of particles:
  create(npLocal);
}


// delete M particles, starting with the Ith particle.  If the last argument
// is true, the destroy will be done immediately, otherwise the request
// will be cached.
template<class PLayout>
void ParticleBase<PLayout>::destroy(size_t M, size_t I, bool doNow) {
  TAU_PROFILE("ParticleBase::destroy()", 
    "void (unsigned, unsigned, bool)", TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  if (M > 0) {
    if (doNow) {
      // find out if we are using optimized destroy method
      bool optDestroy = getUpdateFlag(PLayout::OPTDESTROY);
      // loop over attributes and carry out the destroy request
      attrib_container_t::iterator abeg, aend = AttribList.end();
      for (abeg = AttribList.begin(); abeg != aend; ++abeg)
        (*abeg)->destroy(M,I,optDestroy);
      LocalNum -= M;
    }
    else {
      // add this group of particle indices to our list of items to destroy
      pair<size_t,size_t> destroyEvent(I,M);
      DestroyList.push_back(destroyEvent);
      DestroyNum += M;
    }

    // remember we have this many more items to destroy (or have destroyed)
    ADDIPPLSTAT(incParticlesDestroyed,M);
  }
}


// Update the particle object after a timestep.  This routine will change
// our local, total, create particle counts properly.
template<class PLayout>
void ParticleBase<PLayout>::update() {
  TAU_TYPE_STRING(taustr, CT(*this) + " void ()"); 
  TAU_PROFILE("ParticleBase::update()", taustr, TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  // ask the layout manager to update our atoms, etc.
  Layout->update(*this);
  INCIPPLSTAT(incParticleUpdates);
}


// Update the particle object after a timestep.  This routine will change
// our local, total, create particle counts properly.
template<class PLayout>
void ParticleBase<PLayout>::update(const ParticleAttrib<char>& canSwap) {
  TAU_TYPE_STRING(taustr, CT(*this) + " void (ParticleAttrib<char>)"); 
  TAU_PROFILE("ParticleBase::update()", taustr, TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  // ask the layout manager to update our atoms, etc.
  Layout->update(*this, &canSwap);
  INCIPPLSTAT(incParticleUpdates);
}


// Actually perform the delete atoms action for all the attributes; the
// calls to destroy() only stored a list of what to do.  This actually
// does it.  This should in most cases only be called by the layout manager.
template<class PLayout>
void ParticleBase<PLayout>::performDestroy() {
  TAU_TYPE_STRING(taustr, CT(*this) + " void ()"); 
  TAU_PROFILE("ParticleBase::performDestroy()", taustr, TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  // nothing to do if destroy list is empty
  if (DestroyList.empty()) return;

  // before processing the list, we should make sure it is sorted
  bool isSorted = true;
  typedef vector< pair<size_t,size_t> > dlist_t;
  dlist_t::const_iterator curr = DestroyList.begin(),
                          last = DestroyList.end();
  dlist_t::const_iterator next = curr + 1;
  while (next != last && isSorted) {
    if (*next++ < *curr++) isSorted = false;
  }
  if (!isSorted)
     ::sort(DestroyList.begin(),DestroyList.end());
 
  // find out if we are using optimized destroy method
  bool optDestroy = getUpdateFlag(PLayout::OPTDESTROY);

  // loop over attributes and process destroy list
  attrib_container_t::iterator abeg, aend = AttribList.end();
  for (abeg = AttribList.begin(); abeg != aend; ++abeg)
    (*abeg)->destroy(DestroyList,optDestroy);

  // clear destroy list and update destroy num counter
  DestroyList.erase(DestroyList.begin(),DestroyList.end());
  DestroyNum = 0;
}


// delete M ghost particles, starting with the Ith particle.
// This is done immediately.
template<class PLayout>
void ParticleBase<PLayout>::ghostDestroy(size_t M, size_t I) {
  TAU_TYPE_STRING(taustr, CT(*this) + " void (unsigned, unsigned)"); 
  TAU_PROFILE("ParticleBase::ghostDestroy()", taustr, TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  if (M > 0) {
    // delete the data from the attribute containers
    size_t dnum = 0;
    attrib_container_t::iterator abeg = AttribList.begin();
    attrib_container_t::iterator aend = AttribList.end();
    for ( ; abeg != aend; ++abeg )
      dnum = (*abeg)->ghostDestroy(M, I);
    GhostNum -= dnum;
  }
}


// Put the data for M particles starting from local index I in a Message.
// Return the number of particles put in the Message.  This is for building
// ghost particle interaction lists.
template<class PLayout>
size_t
ParticleBase<PLayout>::ghostPutMessage(Message &msg, size_t M, size_t I) {
  TAU_PROFILE("ParticleBase::ghostPutMessage()", 
              "unsigned (Message, unsigned, unsigned)", TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  // put into message the number of items in the message
  if (I >= R.size()) {
    // we're putting in ghost particles ...
    if ((I + M) > (R.size() + GhostNum))
      M = (R.size() + GhostNum) - I;
  } else {
    // we're putting in local particles ...
    if ((I + M) > R.size())
      M = R.size() - I;
  }
  msg.put(M);

  // go through all the attributes and put their data in the message
  if (M > 0) {
    attrib_container_t::iterator abeg = AttribList.begin();
    attrib_container_t::iterator aend = AttribList.end();
    for ( ; abeg != aend; abeg++ )
      (*abeg)->ghostPutMessage(msg, M, I);
  }

  return M;
}


// put the data for particles on a list into a Message, given list of indices
// Return the number of particles put in the Message.  This is for building
// ghost particle interaction lists.
template<class PLayout>
size_t
ParticleBase<PLayout>::ghostPutMessage(Message &msg,
                                       const vector<size_t>& pl) {
  TAU_PROFILE("ParticleBase::ghostPutMessage()", 
              "unsigned (Message, vector<unsigned>)", TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  vector<size_t>::size_type M = pl.size();  
  msg.put(M);

  // go through all the attributes and put their data in the message
  if (M > 0) {
    attrib_container_t::iterator abeg = AttribList.begin();
    attrib_container_t::iterator aend = AttribList.end();
    for ( ; abeg != aend; ++abeg )
      (*abeg)->ghostPutMessage(msg, pl);
  }

  return M;
}


// retrieve particles from the given message and sending node and store them
template<class PLayout>
size_t
ParticleBase<PLayout>::ghostGetMessage(Message& msg, int node) {
  TAU_TYPE_STRING(taustr, CT(*this) + " unsigned (Message, int)"); 
  TAU_PROFILE("ParticleBase::ghostGetMessage()", taustr, TAU_PARTICLE);

  // make sure we've been initialized
  PAssert(Layout != 0);

  // get the number of items in the message
  size_t numitems;
  msg.get(numitems);
  GhostNum += numitems;

  // go through all the attributes and get their data from the message
  if (numitems > 0) {
    attrib_container_t::iterator abeg = AttribList.begin();
    attrib_container_t::iterator aend = AttribList.end();
    for ( ; abeg != aend; abeg++ )
      (*abeg)->ghostGetMessage(msg, numitems);
  }

  return numitems;
}


// Apply the given sort-list to all the attributes.  The sort-list
// may be temporarily modified, thus it must be passed by non-const ref.
template<class PLayout>
void ParticleBase<PLayout>::sort(SortList_t &sortlist) {
  attrib_container_t::iterator abeg = AttribList.begin();
  attrib_container_t::iterator aend = AttribList.end();
  for ( ; abeg != aend; ++abeg )
    (*abeg)->sort(sortlist);
}


// print it out
template<class PLayout>
ostream& operator<<(ostream& out, const ParticleBase<PLayout>& P) {
  TAU_TYPE_STRING(taustr, "ostream (ostream, " + CT(P) + " )" ); 
  TAU_PROFILE("operator<<()", taustr, TAU_PARTICLE | TAU_IO);
  out << "Particle object contents:";
  out << "\n  Total particles: " << P.getTotalNum();
  out << "\n  Local particles: " << P.getLocalNum();
  out << "\n  Attributes (including R and ID): " << P.numAttributes();
  out << "\n  Layout = " << P.getLayout();
  return out;
}


// print out debugging information
template<class PLayout>
void ParticleBase<PLayout>::printDebug(Inform& o) {
  TAU_PROFILE("printDebug()", "void (Inform)", TAU_PARTICLE | TAU_IO);
  o << "PBase: total = " << getTotalNum() << ", local = " << getLocalNum();
  o << ", attributes = " << AttribList.size() << endl;
  for (attrib_container_t::size_type i=0; i < AttribList.size(); ++i) {
    o << "    ";
    AttribList[i]->printDebug(o);
    o << endl;
  }
  o << "    ";
  Layout->printDebug(o);
}


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

---