dataLL.tar

dataLL.h0000755000175000017500000003114111610330336012273 0ustar  mrobinsmrobins
#ifndef DATALL_H
#define DATALL_H

#include <vector>
#include <list>
#include <blitz/array.h>
#include <mpi.h>
#include <sys/time.h>
#include "customConstants.h"
#include "particle.h"
#include "sphConstants.h"
#include "globalVars.h"


class CpInfoLink;

/*
 * CpInfo contains information on a particle (Cparticle) that is used by
 * CdataLL. This information is: 
 *            - key        used to spatially sort the particles
 *            - neighbrs   list of neighbours     
 *            - ppInfoLink link to entry for particle in pInfoLink.
 */
class CpInfo {
public:
      Cparticle *p;
      CpInfoLink *ppInfoLink;
      unsigned long key;
      vector<Cparticle *> neighbrs;

      CpInfo() {}

      CpInfo(const CpInfo &pInfo) {
         p = pInfo.p;
         ppInfoLink = pInfo.ppInfoLink;
         key = pInfo.key;
      }

      CpInfo& operator=(CpInfo &pInfo) {
         p = pInfo.p;
         ppInfoLink = pInfo.ppInfoLink;
         key = pInfo.key;
         return *this;
      }
};

/*
 * contains back-link to particle information. Used to implement vector of
 * pointers pInfoLinks that exists alongside ps, the vector of particles. pInfoLinks points to particle
 * information in pInfos, a list that is spatially sorted. This way the CdataLL
 * class can loop through a sorted list of particle (pInfos) or an unsorted
 * list (ps or pInfoLinks). 
 */
class CpInfoLink {
public:
      CpInfo *ppInfo;
};

/*
* helper classes for the functions
*     CdataLL::getParticleFromTag(int tag);
*     CdataLL::getPinfoFromTag(int tag);
*/
class CtagPP {
public:
     int tag;
     Cparticle *p;
     CtagPP() {
        tag = -1;
     }
};
class CtagPI {
public:
      int tag;
      CpInfo *pInfo;
      CtagPI() {
         tag = -1;
      }
};


/*
 * the main particle data class. This handles:
 *    -storing the particle data in ps
 *    -applying functions to particles (traverse())
 *    -applying functions to particles and their neighbours (neighbours() and neighboursGroups())
 *    -providing a sorted (pInfos) and unsorted (pInfoLinks and ps) list of
 *        particles and particle information (neighbours etC). 
 *        Sorting done spatially to improve cache performance
 *    -mpi communication of particle and ghost particles to neighbouring cpus
 *    -get particle by tag
 *    -misc others....
 */
class CdataLL {
public:

      /*
       * constructor takes a vector of particles (in_ps), global variables
       * (in_globals) and a flag that determines whether spatial ordering is to be performed
       */
      CdataLL(particleContainer &in_ps, CglobalVars &in_globals,bool opt);
      
      /*
       * reset - Whenever the particle positions are changed the neighbour
       * lists and the ghost particles from neighbouring cpus are no longer valid.
       * This function resets these, and all other data contained in dataLL
       * that depends on the particle positions.
       * The user would normally call this after the particles have been moved.
       */
      void reset();

      /*
       * traverse - This applies a function (Thefunct) to all the particles.
       *            [optional] -can supply a function (Iffunct) that only
       *                        applies Thefunct if Iffunct returns true.
       *                       -can supply input constants or state variables.
       *                        see below for details....
       *
       * big TODO: all these different options have evolved over the life of
       *           this code. They can, and should, be combined into a smaller
       *           number (optimally just one) of functions. 
       */
      template <void Thefunct(Cparticle &,CglobalVars &)>
      void traverse();
      template <void Thefunct(Cparticle &,CglobalVars &,double)>
      void traverse(double);
      template <class T, void Thefunct(Cparticle &,CglobalVars &,T &)>
      void traverse(T &);
      template <void Thefunct(Cparticle &,CglobalVars &), bool Iffunct(Cparticle &)>
      void traverse();
      template <void Thefunct(Cparticle &,CglobalVars &,double), bool Iffunct(Cparticle &)>
      void traverse(double);
      template <class T, void Thefunct(Cparticle &,CglobalVars &,T &), bool Iffunct(Cparticle &)>
      void traverse(T &);
      template <void Thefunct(Cparticle &,CglobalVars &,double,double)>
      void traverse(double,double);
      template <void Thefunct(Cparticle &,CglobalVars &,double,double), bool Iffunct(Cparticle &)>
      void traverse(double,double);
      template <void Thefunct(Cparticle &,CglobalVars &,double,double,double), bool Iffunct(Cparticle &)>
      void traverse(double,double,double);


      /*
       * neighbours - This applies a function (Thefunct) to all the particle
       *              pairs that are neighbours with each other
       *              [optional] -can supply a function (Iffunct) that only
       *                        applies Thefunct if Iffunct returns true.
       */
      template <void Thefunct(Cparticle &,Cparticle &,CglobalVars &)>
      void neighbours();
      template <void Thefunct(Cparticle &,Cparticle &,CglobalVars &), bool Iffunct(Cparticle &)>
      void neighbours();

      /*
       * neighboursGroup - This applies a function (Thefunct) to all the neighbours
       *              of each particle. The neighbours are passed to the
       *              function using a vector of particle pointers.
       *              [optional] -can supply a function (Iffunct) that only
       *                          applies Thefunct if Iffunct returns true.
       *                         -can supply a state variable (with an
       *                          arbitrary class)
       */
      template <void Thefunct(Cparticle &,vector<Cparticle *> &,CglobalVars &)>
      void neighboursGroup();
      template <void Thefunct(Cparticle &,vector<Cparticle *> &,CglobalVars &), bool Iffunct(Cparticle &)>
      void neighboursGroup();
      template <class T, void Thefunct(Cparticle &,vector<Cparticle *> &,CglobalVars &,T &), bool Iffunct(Cparticle &)>
      void neighboursGroup(T &);

      /*
       * neighboursUsing - This applies a function (Thefunct) to all the particles
       *              that are neighbours with each particle in _ps.
       *              [optional] -Thefunct can either take a particle pair or a
       *                          particle and all its neighbours (as a vector of particle pointers).
       */
      template <void Thefunct(Cparticle &,Cparticle &,CglobalVars &)>
      void neighboursUsing(vector<Cparticle> &_ps);
      template <void Thefunct(Cparticle &,vector<Cparticle *> &,CglobalVars &)>
      void neighboursUsing(vector<Cparticle> &_ps);
      
      /*
       * syncParticlesBetweenProcs - syncs data between neighbouring cpus. The
       * data to be synced is read/writen using the templated functions (to be
       * applied to each particles).
       */
      template <class T, T readFunct(Cparticle &), void writeFunct(Cparticle &,T)>
      void syncParticlesBetweenProcs();

      /*
       * functionOverGrid - apply theFunct over a set of grid points and the
       * particles that neighbour these grid points. theFunct
       * takes as input a particle pair, the first particle is a grid point, 
       * the second is a neighbour of that grid point. The final grid points
       * are output in outPs
       */
      template <void theFunct(Cparticle &,Cparticle &,CglobalVars &)>
      void functionOverGrid(vector<Cparticle> &outPs,vectInt &gridDimN);

      /*
       * find neighbours of particle _p. neighbours returned in _neighbrs
       */
      void calcNeighbours(vector<Cparticle *> &_neighbrs,Cparticle &_p);

      /*
       * set a timestep (sent to all cpus, the minimum timestep given by all
       * cpus is used)
       */
      void setGlobalTimestep(double dt);

      /*
       * print info on dataLL class
       */
      void printMe();

      /*
       * update the spatial ordering of the particles. This improves cache
       * performance. Only use once particles are sufficiently disordered.
       * (eg. every output timestep)
       */
      void calcTraverseOrder();

      int numParticles() { return ps.size(); };
      particleContainer *getParticles() { return &ps; }

      /*
       * useful if you have just made massive changes to the particle array 
       * in general or particle positions in general. Reconstructs the spatial
       * ordering from scratch (warning: very slow!)
       */
      void updateParticles();

      /*
       * insertion/deletion functions
       */
      void insertNewParticle(Cparticle &newP);
      void markForDeletion(Cparticle &p);
      void deleteParticles();

      /*
       * get particle or particle info from a (unique) tag
       */
      Cparticle *getParticleFromTag(int tag);
      CpInfo *getPinfoFromTag(int tag);

      /*
       * sum a value over all cpus
       */
      void sumOverProcs(vect *data,int num);
      void sumOverProcs(double *data,int num);

      /*
       * global variables. TODO: this should be private
       */
      CglobalVars &globals;

private:

      /*
       * functions to allocate mpi buffers
       */
      int calcBufferSize(vectInt coords);
      int calcPBufferSize(vectInt coords);
      void allocateBuffers();

      /*
       * delete particle from ps, pInfos and pInfoLink
       */
      void deleteParticle(vector<CpInfoLink>::iterator &ppInfoLink,CpInfo *ppInfo);

      /*
       * insert a particle into the bucket array
       */
      void insert(Cparticle &);

      void calcDomainLimits();

      /*
       * add _toAdd to list of neighbours of _p, only if _toAdd is a neighbour
       * of _p
       */
      void addIfNeighbr(Cparticle &_p,vector<Cparticle *> &_neighbrs,vector<Cparticle *> &_toAdd);

      /*
       * maps from position to bucket cell index
       */
      inline vectInt getCellI(vect r) {
         return static_cast<vectInt>((r-rmin)/(KERNAL_RADIUS*hmax)+2);
      }

      /*
       * functions to sort particles spatially
       *    calcKey. calculates a key for each particle using morton ordering
       *    initPinfos - clears current pInfos (vector of particle information,
       *    like neighbours, keys etc. and recreates it
       *    checkPinfos - checks that pInfos is consistent with pInfoLinks
       *    sortPInfos - sorts pInfos according to key
       *    initTraverseOrder - runs previous three functions.
       */
      unsigned long calcKey(vect r);
      void initPinfos();
      void checkPinfos();
      void sortPInfos();
      void initTraverseOrder();
 
      /*
       * functions that implement mpi communication
       */
      void updateDomain();
      void sendRecvGhosts(const vectInt coords,const int num,MPI_Request *request,int *ghostSizesSend,int *ghostSizesRecv);
      void sendRecvParticles(const vectInt coords,const int num,MPI_Request *request,int *particleSizesSend,int *particleSizesRecv);
      void sendRecvDataSync(const vectInt coords,const int num,MPI_Request *request,int sendSize, void *sendBuffer, int recvSize, void *recvBuffer);
      
      /*
       * vectors that hold pointers to particles or pInfos. These pointers are
       * sorted by tag, so that the particle or pInfo with a given tag can easily be found
       */
      vector<CtagPP> tagSortedParticlePointers;
      vector<CtagPI> tagSortedPinfoPointers;


      /*
       * domain extents and max smoothing length
       */
      vect rmax;
      vect rmin;
      double hmax;

      Array<vector<Cparticle *>,NDIM> cells; //bucket list
      int n;                                 //number of particles
      particleContainer &ps;                 //holds the set of particles
      Array<Cparticle **,NDIM> ghostedParticles; //points to the particles ghosted to neighbouring cpus
      vector<Cparticle> procGhostParticles;  //ghost particles received from neighbouring cpus
      Array<double,NDIM> procGhostMax2H;     //max 2h from neighbouring cpus


      /*
       * buffers and buffer sizes for mpi communication
       */
      Array<CghostData*,NDIM> ghostBuffersSend;
      Array<CghostData*,NDIM> ghostBuffersRecv;
      Array<Cparticle*,NDIM> particleBuffersSend;
      Array<Cparticle*,NDIM> particleBuffersRecv;
      Array<void*,NDIM> syncBuffersSend;
      Array<void*,NDIM> syncBuffersRecv;

      Array<int,NDIM> pBufferSizes;
      Array<int,NDIM> gBufferSizes;
      Array<int,NDIM> sendSizesGhosts;
      Array<int,NDIM> recvSizesGhosts;

      /*
       * pInfos holds particle information that is used by this class
       * (neighbours, sorting key etc.). 
       * pInfoLinks is a vector where the i'th entry contains a pointer to the
       * pInfo for particle i in the main particle vector (ps)
       */
      vector<CpInfoLink> pInfoLinks;
      vector<CpInfo> pInfos;
};

#endif
dataLL.impl.h0000755000175000017500000003271211610061124013234 0ustar  mrobinsmrobins#ifndef DATALL_IMPL_H
#define DATALL_IMPL_H

#include "dataLL.h"
#include "misc.h"


/*
 * basic traverse function. Use this to apply a per-particle function to the
 * set of particles (eg. a drift or kick operation)
 * loops through the ps vector (vector of Cparticle) and runs T_thefunct on
 * each particle
 * also loops through all the ghost particles (procGhostParticles) and does
 * the same.
 */
template <void T_thefunct(Cparticle &,CglobalVars &)>
void CdataLL::traverse() {
   for (particleContainer::iterator pp = ps.begin();pp!=ps.end();pp++) {
      T_thefunct(*pp,globals);
   }
   for (vector<Cparticle>::iterator pp = procGhostParticles.begin();pp!=procGhostParticles.end();pp++) {
      T_thefunct(*pp,globals);
   }
}

/*
 * same as above but adds an input constant to T_thefunct
 */
template <void T_thefunct(Cparticle &,CglobalVars &, double)>
void CdataLL::traverse(double param) {
   for (particleContainer::iterator pp = ps.begin();pp!=ps.end();pp++) {
      T_thefunct(*pp,globals,param);
   }
   for (vector<Cparticle>::iterator pp = procGhostParticles.begin();pp!=procGhostParticles.end();pp++) {
      T_thefunct(*pp,globals,param);
   }
}

/*
 * same as above but only applies the function if the boolean Iffunct is true
 */
template <void T_thefunct(Cparticle &,CglobalVars &, double),bool Iffunct(Cparticle &)>
void CdataLL::traverse(double param) {
   for (particleContainer::iterator pp = ps.begin();pp!=ps.end();pp++) {
      if (Iffunct(*pp)) {
         T_thefunct(*pp,globals,param);
      }
   }
   //TODO: ghost or not should be indicated by another flag, but this would be a significant change
   if (!procGhostParticles.empty()) {
      enum iamTypes store = procGhostParticles.begin()->iam;
      procGhostParticles.begin()->iam = ghost;
      if (Iffunct(*(procGhostParticles.begin()))) {
         procGhostParticles.begin()->iam = store;
         for (vector<Cparticle>::iterator pp = procGhostParticles.begin();pp!=procGhostParticles.end();pp++) {
            if (Iffunct(*pp)) {
               T_thefunct(*pp,globals,param);
            }
         }
      }
      procGhostParticles.begin()->iam = store;
   }
}

/*
 * and other variants. Look at header file for info...
 */
template <class T, void Thefunct(Cparticle &,CglobalVars &,T &)>
void CdataLL::traverse(T &param) {
   for (particleContainer::iterator pp = ps.begin();pp!=ps.end();pp++) {
      Thefunct(*pp,globals,param);
   }
   for (vector<Cparticle>::iterator pp = procGhostParticles.begin();pp!=procGhostParticles.end();pp++) {
      Thefunct(*pp,globals,param);
   }
}

template <class T, void Thefunct(Cparticle &,CglobalVars &,T &), bool Iffunct(Cparticle &)>
void CdataLL::traverse(T &param) {
   for (particleContainer::iterator pp = ps.begin();pp!=ps.end();pp++) {
      if (Iffunct(*pp)) {
         Thefunct(*pp,globals,param);
      }
   }
   if (!procGhostParticles.empty()) {
      enum iamTypes store = procGhostParticles.begin()->iam;
      procGhostParticles.begin()->iam = ghost;
      if (Iffunct(*(procGhostParticles.begin()))) {
         procGhostParticles.begin()->iam = store;
         for (vector<Cparticle>::iterator pp = procGhostParticles.begin();pp!=procGhostParticles.end();pp++) {
            if (Iffunct(*pp)) {
               Thefunct(*pp,globals,param);
            }
         }
      }
      procGhostParticles.begin()->iam = store;
   }
}

template <void T_thefunct(Cparticle &,CglobalVars &, double, double)>
void CdataLL::traverse(double param1,double param2) {
   for (particleContainer::iterator pp = ps.begin();pp!=ps.end();pp++) {
      T_thefunct(*pp,globals,param1,param2);
   }
   for (vector<Cparticle>::iterator pp = procGhostParticles.begin();pp!=procGhostParticles.end();pp++) {
      T_thefunct(*pp,globals,param1,param2);
   }
}

template <void T_thefunct(Cparticle &,CglobalVars &, double, double),bool Iffunct(Cparticle &)>
void CdataLL::traverse(double param1, double param2) {
   for (particleContainer::iterator pp = ps.begin();pp!=ps.end();pp++) {
      if (Iffunct(*pp)) {
         T_thefunct(*pp,globals,param1,param2);
      }
   }
   if (!procGhostParticles.empty()) {
      enum iamTypes store = procGhostParticles.begin()->iam;
      procGhostParticles.begin()->iam = ghost;
      if (Iffunct(*(procGhostParticles.begin()))) {
         procGhostParticles.begin()->iam = store;
         for (vector<Cparticle>::iterator pp = procGhostParticles.begin();pp!=procGhostParticles.end();pp++) {
            if (Iffunct(*pp)) {
               T_thefunct(*pp,globals,param1,param2);
            }
         }
      }
      procGhostParticles.begin()->iam = store;
   }
}

template <void T_thefunct(Cparticle &,CglobalVars &, double, double, double),bool Iffunct(Cparticle &)>
void CdataLL::traverse(double param1, double param2, double param3) {
   for (particleContainer::iterator pp = ps.begin();pp!=ps.end();pp++) {
      if (Iffunct(*pp)) {
         T_thefunct(*pp,globals,param1,param2,param3);
      }
   }
   if (!procGhostParticles.empty()) {
      enum iamTypes store = procGhostParticles.begin()->iam;
      procGhostParticles.begin()->iam = ghost;
      if (Iffunct(*(procGhostParticles.begin()))) {
         procGhostParticles.begin()->iam = store;
         for (vector<Cparticle>::iterator pp = procGhostParticles.begin();pp!=procGhostParticles.end();pp++) {
            if (Iffunct(*pp)) {
               T_thefunct(*pp,globals,param1,param2,param3);
            }
         }
      }
      procGhostParticles.begin()->iam = store;
   }
}

template <void T_thefunct(Cparticle &,CglobalVars &),bool Iffunct(Cparticle &)>
void CdataLL::traverse() {
   for (particleContainer::iterator pp = ps.begin();pp!=ps.end();pp++) {
      if (Iffunct(*pp)) {
         T_thefunct(*pp,globals);
      }
   }
   if (!procGhostParticles.empty()) {
      enum iamTypes store = procGhostParticles.begin()->iam;
      procGhostParticles.begin()->iam = ghost;
      if (Iffunct(*(procGhostParticles.begin()))) {
         procGhostParticles.begin()->iam = store;
         for (vector<Cparticle>::iterator pp = procGhostParticles.begin();pp!=procGhostParticles.end();pp++) {
            if (Iffunct(*pp)) {
               T_thefunct(*pp,globals);
            }
         }
      }
      procGhostParticles.begin()->iam = store;
   }
}

/*
 * basic neighbours function. This time it loops over pInfos, which is sorted
 * spatially. As Thefunct is calculated over neighbouring particles pairs, this
 * will improve cache performance, as close particles have already been
 * accessed.
 */
template <void Thefunct(Cparticle &,Cparticle &,CglobalVars &)>
void CdataLL::neighbours() {
   for (vector<CpInfo>::iterator ppInfo = pInfos.begin();ppInfo!=pInfos.end();ppInfo++) {
      Cparticle *pp = ppInfo->p;
      if (ppInfo->neighbrs.empty()) {
         calcNeighbours(ppInfo->neighbrs,*pp);
         //cout << "Found "<<neighbrs[i].size()<<" neighbours for particle "<<i<<endl;
      }
      //for each neighbours run function
      for (vector<Cparticle *>::iterator pNeighbr = ppInfo->neighbrs.begin();pNeighbr!=ppInfo->neighbrs.end();pNeighbr++) {
         Thefunct(*pp,**pNeighbr,globals);
      }
   }
}

/*
 * same as above but adds an if function...
 */
template <void Thefunct(Cparticle &,Cparticle &,CglobalVars &),bool Iffunct(Cparticle &)>
void CdataLL::neighbours() {
   for (vector<CpInfo>::iterator ppInfo = pInfos.begin();ppInfo!=pInfos.end();ppInfo++) {
      Cparticle *pp = ppInfo->p;
      if (Iffunct(*pp)) {
         if (ppInfo->neighbrs.empty()) {
            calcNeighbours(ppInfo->neighbrs,*pp);
         }
         //for each neighbours run function
         for (vector<Cparticle *>::iterator pNeighbr = ppInfo->neighbrs.begin();pNeighbr!=ppInfo->neighbrs.end();pNeighbr++) {
            Thefunct(*pp,**pNeighbr,globals);
         }
      }
   }
}

template <void Thefunct(Cparticle &,vector<Cparticle *> &,CglobalVars &)>
void CdataLL::neighboursGroup() {
   for (vector<CpInfo>::iterator ppInfo = pInfos.begin();ppInfo!=pInfos.end();ppInfo++) {
      particleContainer::iterator pp = ppInfo->p;
      if (ppInfo->neighbrs.empty()) {
         calcNeighbours(ppInfo->neighbrs,*pp);
      }
      Thefunct(*pp,ppInfo->neighbrs,globals);
   }
}

template <void Thefunct(Cparticle &,vector<Cparticle *> &,CglobalVars &), bool Iffunct(Cparticle &)>
void CdataLL::neighboursGroup() {
   for (vector<CpInfo>::iterator ppInfo = pInfos.begin();ppInfo!=pInfos.end();ppInfo++) {
      Cparticle *pp = ppInfo->p;
      if (Iffunct(*pp)) {
         if (ppInfo->neighbrs.empty()) {
            calcNeighbours(ppInfo->neighbrs,*pp);
         }
         Thefunct(*pp,ppInfo->neighbrs,globals);
      }
   }
}


template <class T, void Thefunct(Cparticle &,vector<Cparticle *> &,CglobalVars &,T &), bool Iffunct(Cparticle &)>
void CdataLL::neighboursGroup(T &param) {
   for (vector<CpInfo>::iterator ppInfo = pInfos.begin();ppInfo!=pInfos.end();ppInfo++) {
      Cparticle *pp = ppInfo->p;
      if (Iffunct(*pp)) {
         if (ppInfo->neighbrs.empty()) {
            calcNeighbours(ppInfo->neighbrs,*pp);
         }
         Thefunct(*pp,ppInfo->neighbrs,globals,param);
      }
   }
}

template <void Thefunct(Cparticle &,Cparticle &,CglobalVars &)>
void CdataLL::neighboursUsing(vector<Cparticle> &_ps) {
   vector<Cparticle *> _neighbrs;
   for (int i=0;i<_ps.size();i++) {
      _neighbrs.clear();
      calcNeighbours(_neighbrs,_ps[i]);
      //for each neighbours run function
      for (int j=0;j<_neighbrs.size();j++) {
         Cparticle *_pthep = _neighbrs[j];
         Thefunct(_ps[i],*_pthep,globals);
      }
   }
}

template <void Thefunct(Cparticle &,vector<Cparticle *> &,CglobalVars &)>
void CdataLL::neighboursUsing(vector<Cparticle> &_ps) {
   vector<Cparticle *> _neighbrs;
   int n = _ps.size();
   for (int i=0;i<n;i++) {
      _neighbrs.clear();
      calcNeighbours(_neighbrs,_ps[i]);
      Thefunct(_ps[i],_neighbrs,globals);
   }
}


template <void theFunct(Cparticle &,Cparticle &,CglobalVars &)>
void CdataLL::functionOverGrid(vector<Cparticle> &outPs,vectInt &gridDimN) {
   //TODO: I've broken this so that it won't work in parallel anymore. To fix just replace newRmax and min with rmax and rmin
   vect newRmax;
   vect newRmin;
   for (int i=0;i<NDIM;i++) {
      newRmax = RMAX[i];
      newRmin = RMIN[i];
   }
   gridDimN = (newRmax-newRmin)/GRIDSEP;
   vect newGridSep = (newRmax-newRmin)/gridDimN;
   gridDimN += 1;
   int gridN = product(gridDimN); 
   outPs.resize(gridN);
  
   for (int i=0;i<gridDimN[0];i++) {
      if (NDIM > 1) {
         for (int j=0;j<gridDimN[1];j++) {
            int jOffset = j*gridDimN[0];
            if (NDIM > 2) {
               for (int k=0;k<gridDimN[2];k++) {
                  int kOffset = k*gridDimN[0]*gridDimN[1];
                  outPs[kOffset+jOffset+i].r = i*newGridSep[0],j*newGridSep[1],k*newGridSep[2];
                  outPs[kOffset+jOffset+i].r += newRmin; 
               }
            } else {
               outPs[jOffset+i].r = i*newGridSep[0],j*newGridSep[1];
               outPs[jOffset+i].r += newRmin; 
            }
         }
      } else {
         outPs[i].r = i*newGridSep[0];
         outPs[i].r += newRmin; 
      } 
   }
   neighboursUsing<theFunct>(outPs);
}





template <class T, T readFunct(Cparticle &), void writeFunct(Cparticle &,T)>
void CdataLL::syncParticlesBetweenProcs() {
//NOTE: if the size sent per particle is greater than sizeof(CghostData) the send and recv buffers are too small!

   //fill up send buffers
   for (Array<int,NDIM>::iterator ap = globals.procNeighbrs.begin();ap != globals.procNeighbrs.end();ap++) {
      vectInt coords = ap.position();
      if (globals.procNeighbrs(coords)>=0) {
         T *p = (T *)syncBuffersSend(coords);
         for (int j=0;j<sendSizesGhosts(coords);j++) {
            *p = readFunct(*(ghostedParticles(coords)[j]));
            p++;
         }
      }
   }
  
   //send and recv everything (non-blocking)
   int sendRecvArraySize = 2*int(pow(3.0,NDIM));
   MPI_Request requestSendRecv[sendRecvArraySize];
   MPI_Status statusSendRecv[sendRecvArraySize];
   if (globals.mpiRank%2==0) {
      int upper = int(pow(3.0,NDIM));
      for (int i=0;i<upper;i++) {
         vectInt split = 3;
         vectInt coords = Nmisc::numToCoords(i,split);
         if (globals.procNeighbrs(coords)>=0) {
            sendRecvDataSync(coords,i,requestSendRecv+i*2,sizeof(T)*sendSizesGhosts(coords),syncBuffersSend(coords),sizeof(T)*recvSizesGhosts(coords),syncBuffersRecv(coords));
         }
         else {
            for (int j=0;j<2;j++) {
               requestSendRecv[i*2+j] = MPI_REQUEST_NULL;
            }
         }
      }
   } else {
      for (int i=int(pow(3.0,NDIM))-1;i>=0;i--) {
         vectInt split = 3;
         vectInt coords = Nmisc::numToCoords(i,split);
         if (globals.procNeighbrs(coords)>=0) {
            sendRecvDataSync(coords,i,requestSendRecv+i*2,sizeof(T)*sendSizesGhosts(coords),syncBuffersSend(coords),sizeof(T)*recvSizesGhosts(coords),syncBuffersRecv(coords));
         }
         else {
            for (int j=0;j<2;j++) {
               requestSendRecv[i*2+j] = MPI_REQUEST_NULL;
            }
         }
      } 
   }
   MPI_Waitall(sendRecvArraySize, requestSendRecv, statusSendRecv);   

   //sync new data with ghost particles
   vector<Cparticle>::iterator gp = procGhostParticles.begin();
   for (Array<int,NDIM>::iterator ap = globals.procNeighbrs.begin();ap != globals.procNeighbrs.end();ap++) {
      vectInt coords = ap.position();
      if (globals.procNeighbrs(coords)>= 0) {
         T *buffPtrT = (T *)syncBuffersRecv(coords);
         for (int z=0;z<recvSizesGhosts(coords);z++) {
            writeFunct(*gp,*buffPtrT);  
            buffPtrT++;
            gp++;
         }
      }
   }
}

#endif
 
dataLL.cpp0000755000175000017500000010214411610331612012626 0ustar  mrobinsmrobins#include "dataLL.impl.h"
#include "sort.h"

const unsigned int KEY_SIZE = 30;  // must be divisable by NDIM
const unsigned int KEY_DIM_BITS = KEY_SIZE/NDIM;
const unsigned int KEY_DIM_MAX = (unsigned int)floor(pow(2.0,(int)KEY_DIM_BITS)-1);

CdataLL::CdataLL(particleContainer &in_ps, CglobalVars &in_globals,bool opt): ps(in_ps),globals(in_globals) {
   //initialise hmax, the maximum smoothing length
   hmax = H;

   // allocate buffers used for mpi communication
   allocateBuffers();
         
   //keep track of hmax for each neighbouring cpu
   procGhostMax2H.resize(3);
   procGhostMax2H = 1.1*KERNAL_RADIUS*H;

   cout << "Processor "<<globals.mpiRank<<" of "<<globals.mpiSize<<" has domain ";
   for (int i=0;i<NDIM*2;i++) {
      cout<<globals.procDomain[i]<<' ';
   }

   // if prompted, init the sorted traverse order
   if (opt) {
      initTraverseOrder();
   }

   // reset all data (neighbours, ghosts, mpi comms etc)
   reset();
}

/*
 * allocate buffers used for mpi communication. 
 * TODO: these buffers take up quite a bit of memory. Some could be combined,
 *       for example the ones used to sync info and the normal ghost buffers
 */
void CdataLL::allocateBuffers() {
   ghostBuffersSend.resize(3);
   pBufferSizes.resize(3);
   gBufferSizes.resize(3);

   /*
    * sending and receiving ghost buffers 
    */
   for (Array<CghostData*,NDIM>::iterator i=ghostBuffersSend.begin();i!=ghostBuffersSend.end();i++) {
      vectInt coords = i.position()-1;
      if (globals.procNeighbrs(i.position()) >= 0) {
         *i = new CghostData[calcBufferSize(coords)];
         cout << calcBufferSize(coords)<<endl;
         gBufferSizes(i.position()) = calcBufferSize(coords);
      }
   }

   ghostBuffersRecv.resize(3);
   for (Array<CghostData*,NDIM>::iterator i=ghostBuffersRecv.begin();i!=ghostBuffersRecv.end();i++) {
      vectInt coords = i.position()-1;
      if (globals.procNeighbrs(i.position()) >= 0) {
         *i = new CghostData[calcBufferSize(coords)];
         cout << calcBufferSize(coords)<<endl;
      }
   }

   /*
    * ghostedParticles contains links to particles that have been ghosted
    * (for later syncing purposes)
    */
   ghostedParticles.resize(3);
   for (Array<Cparticle**,NDIM>::iterator i=ghostedParticles.begin();i!=ghostedParticles.end();i++) {
      vectInt coords = i.position()-1;
      if (globals.procNeighbrs(i.position()) >= 0) {
         *i = new Cparticle*[calcBufferSize(coords)];
         cout << calcBufferSize(coords)<<endl;
      }
   }

   /*
    * sending and receiving particle buffers 
    */
   particleBuffersSend.resize(3);
   for (Array<Cparticle*,NDIM>::iterator i=particleBuffersSend.begin();i!=particleBuffersSend.end();i++) {
      vectInt coords = i.position()-1;
      if (globals.procNeighbrs(i.position()) >= 0) {
         *i = new Cparticle[calcPBufferSize(coords)];
         cout << calcPBufferSize(coords)<<endl;
         pBufferSizes(i.position()) = calcPBufferSize(coords);
      }
   }

   particleBuffersRecv.resize(3);
   for (Array<Cparticle*,NDIM>::iterator i=particleBuffersRecv.begin();i!=particleBuffersRecv.end();i++) {
      vectInt coords = i.position()-1;
      if (globals.procNeighbrs(i.position()) >= 0) {
         *i = new Cparticle[calcPBufferSize(coords)];
         cout << calcPBufferSize(coords)<<endl;
      }
   }
   
   /*
    * sending and receiving buffers for syncing
    */
   syncBuffersSend.resize(3);
   for (Array<void*,NDIM>::iterator i=syncBuffersSend.begin();i!=syncBuffersSend.end();i++) {
      vectInt coords = i.position()-1;
      if (globals.procNeighbrs(i.position()) >= 0) {
         *i = new unsigned char[sizeof(CghostData)*calcBufferSize(coords)];
         cout << calcBufferSize(coords)<<endl;
      }
   }
   syncBuffersRecv.resize(3);
   for (Array<void*,NDIM>::iterator i=syncBuffersRecv.begin();i!=syncBuffersRecv.end();i++) {
      vectInt coords = i.position()-1;
      if (globals.procNeighbrs(i.position()) >= 0) {
         *i = new unsigned char[sizeof(CghostData)*calcBufferSize(coords)];
         cout << calcBufferSize(coords)<<endl;
      }
   }
   sendSizesGhosts.resize(3);
   recvSizesGhosts.resize(3);
}


/*
 * reset particle positions, find neighbours, mpi communication to neighbor
 * cpus etc.
 * 
 * user will normally call this whenever particle positions are updated 
 */
void CdataLL::reset() {
   
   tagSortedParticlePointers.clear();

   timeval t1,t2;
   gettimeofday(&t1,NULL);
   int psCapacity = ps.capacity();

   // update domain and sync with neighbouring processes
   updateDomain();

   // check that the sizes of ps, pInfos and pInfoLinks are identical
   if ((ps.size() != pInfos.size())||(ps.size() != pInfoLinks.size())) {
      if (globals.sphStep != 0) {
         cout <<"Error: ps has "<<ps.size()<<" particles, pInfos has "<<pInfos.size()<<" particles, pInfoLinks has "<<pInfoLinks.size()<<" partcles"<<endl;
         exit(-1);
      }
   }

   // ps should not reallocate its memory!
   if (ps.size() > psCapacity) {
      cerr << "Error: dataLL:reset(): ps has grown beyond its capacity! Bad things can happen......I'm outta here!"<<endl;
      exit(-1);
   }
   n = ps.size();

   // recalculate limits
   calcDomainLimits();
   
   // clear bucket list and reinsert particles and ghost particles
   cells.resize(static_cast<vectInt>((rmax-rmin)/(KERNAL_RADIUS*hmax))+6);
   for (Array<vector<Cparticle *>,NDIM>::iterator i=cells.begin();i!=cells.end();i++) {
      i->clear();
   }
   for (particleContainer::iterator p = ps.begin();p != ps.end();p++) {
      insert(*p);
   }
   for (vector<Cparticle>::iterator p=procGhostParticles.begin();p!= procGhostParticles.end();p++) {
      insert(*p);
   }

   gettimeofday(&t2,NULL);
   globals.wtTotalReset.tv_sec += t2.tv_sec-t1.tv_sec;
   globals.wtTotalReset.tv_usec += t2.tv_usec-t1.tv_usec;
}

void CdataLL::insertNewParticle(Cparticle &newP) {
   ps.push_back(newP);
   n++;
   if (n>MAX_NUM_PARTICLES_PER_CPU) {
      cout<<"insertNewParticle(): number of particles greater than  max limit!!!"<<endl;
      exit(-1);
   }

   /*
    * This is the tricky bit. Make sure that the inserted particle info is also
    *placed in pInfos and pInfoLinks
    */
   Cparticle *p = &(ps.back());
   CpInfo newpInfo;
   newpInfo.p = p;
   newpInfo.key = calcKey(newpInfo.p->r);
   pInfos.push_back(newpInfo);
   CpInfoLink newpInfoLink;
   newpInfoLink.ppInfo = &(pInfos.back());
   pInfoLinks.push_back(newpInfoLink);
   pInfos.back().ppInfoLink = &(pInfoLinks.back());
   insert(*p);
}


void CdataLL::markForDeletion(Cparticle &p) {
   // using a custom tag for deletion. Could be better....
   p.tag = -111;
}

/*
 * deletes all particles with a tag = -111. This could be done better...
 * (delete flag maybe)
 */
void CdataLL::deleteParticles() {
   for (vector<CpInfoLink>::iterator ppInfoLink = pInfoLinks.begin();ppInfoLink!=pInfoLinks.end();) {
      CpInfo *ppInfo = ppInfoLink->ppInfo;
      Cparticle *thisP = ppInfo->p;
      if (thisP->tag == -111) {
         deleteParticle(ppInfoLink,ppInfo);
         continue;
      }
      ppInfoLink++;
   }
}

/*
 * deletes particle from particle list (ps). Makes sure that the particle
 * infomation is also deleted from pInfoLinks and pInfos
 */
void CdataLL::deleteParticle(vector<CpInfoLink>::iterator &ppInfoLink,CpInfo *ppInfo) {
   Cparticle *thisP = ppInfo->p;
   if (ppInfoLink != pInfoLinks.end()-1) {
      *thisP = ps.back();
      *ppInfoLink = pInfoLinks.back();

      ppInfoLink->ppInfo->p = thisP;
      ppInfoLink->ppInfo->ppInfoLink = &(*ppInfoLink);

      ps.pop_back();
      pInfoLinks.pop_back();
   } else {
      ps.pop_back();
      pInfoLinks.pop_back();
      ppInfoLink = pInfoLinks.end();
   }

   if (ppInfo != &(pInfos.back())) {
      *ppInfo = pInfos.back();
      ppInfo->ppInfoLink->ppInfo = ppInfo;
      pInfos.pop_back();
   } else {
      pInfos.pop_back();
   }
}


/*
 * BIG TODO: This function handles ALL the mpi processing for each timestep.
 * Need to break this up into sub-functions and tidy it all up
 */
void CdataLL::updateDomain() {
   timeval t11,t22;
   gettimeofday(&t11,NULL);

   Array<int,NDIM> ghostIndicies(3);
   Array<int,NDIM> particleIndicies(3);
   ghostIndicies = 0;
   particleIndicies = 0;

   procGhostParticles.clear();

   /*
    * update relevant domain limits (for ghost particles etc) using
    * globals.procDomain.
    */
   vect dmin,dspace,dgmin,dgspace;
   for (int i=0;i<NDIM;i++) {
      dmin[i] = globals.procDomain[i*2];
      dspace[i] = globals.procDomain[i*2+1] - dmin[i];
      vectInt coords = 1;
      coords[i] = 0;
      dgmin[i] = dmin[i] + procGhostMax2H(coords);
      dgspace[i] = dspace[i] - procGhostMax2H(coords);
      coords[i] = 2;
      dgspace[i] -= procGhostMax2H(coords);
   }

   /*
    * loop thru particles and add all particle to be sent to
    * particleBuffersSend. Delete these particles from this CPUs list
    */
   for (vector<CpInfoLink>::iterator ppInfoLink = pInfoLinks.begin();ppInfoLink!=pInfoLinks.end();) {
      CpInfo *ppInfo = ppInfoLink->ppInfo;

      ppInfo->neighbrs.clear();
      Cparticle *thisP = ppInfo->p;
      if (thisP->tag == -111) {
         deleteParticle(ppInfoLink,ppInfo);
         continue;
      }
      vectInt outCoords = static_cast<vectInt>((thisP->r-dmin)/dspace+1);
      if (any(outCoords>2)||any(outCoords<0)) {
         cout << "Particle out of range! Tag = "<<thisP->tag<<" iam = "<<thisP->iam<<" position = "<<thisP->r<<" velocityhat = "<<thisP->vhat<<" velocity = "<<thisP->v<<" density = "<<thisP->dens<<" total force = "<<thisP->f<<" pressure force = "<<thisP->fp<<endl;
         exit(-1);
      }
      if (globals.procNeighbrs(outCoords) >= 0) {
         //send particle to neighbour
         //add to real particle to send
         if (particleIndicies(outCoords)>=pBufferSizes(outCoords)) {
            cerr << "particleBuffers are full, exiting...."<<endl;
            exit(-1);
         }
         particleBuffersSend(outCoords)[particleIndicies(outCoords)] = *thisP;
         particleIndicies(outCoords)++;
         deleteParticle(ppInfoLink,ppInfo);
         continue;
      }
      ppInfoLink++;
   }

   /*
    * Send these particles to the correct neighbouring cpus
    */
   Array<int,NDIM> particleSizesSend(3);
   Array<int,NDIM> particleSizesRecv(3);
   particleSizesSend = particleIndicies*sizeof(Cparticle);
   int sendRecvArraySize = 4*int(pow(3.0,NDIM));
   MPI_Request requestSendRecv[sendRecvArraySize];
   MPI_Status statusSendRecv[sendRecvArraySize];

   timeval t1,t2;
   gettimeofday(&t1,NULL);
   if (globals.mpiRank%2==0) {
      int upper = int(pow(3.0,NDIM));
      for (int i=0;i<upper;i++) {
         vectInt split = 3;
         vectInt coords = Nmisc::numToCoords(i,split);
         if (globals.procNeighbrs(coords)>=0) {
            sendRecvParticles(coords,i,requestSendRecv+i*4,&(particleSizesSend(coords)),&(particleSizesRecv(coords))); 
         }
         else {
            for (int j=0;j<4;j++) {
               requestSendRecv[i*4+j] = MPI_REQUEST_NULL;
            }
         }
      }
   } else {
      for (int i=int(pow(3.0,NDIM))-1;i>=0;i--) {
         vectInt split = 3;
         vectInt coords = Nmisc::numToCoords(i,split);
         if (globals.procNeighbrs(coords)>=0) {
            sendRecvParticles(coords,i,requestSendRecv+i*4,&(particleSizesSend(coords)),&(particleSizesRecv(coords))); 
         }
         else {
            for (int j=0;j<4;j++) {
               requestSendRecv[i*4+j] = MPI_REQUEST_NULL;
            }
         }
      } 
   }
   MPI_Waitall(sendRecvArraySize, requestSendRecv, statusSendRecv);   

   gettimeofday(&t2,NULL);
   globals.wtTotalMPI.tv_sec += t2.tv_sec-t1.tv_sec;
   globals.wtTotalMPI.tv_usec += t2.tv_usec-t1.tv_usec;
     
   /*
    * All particles have arrived. Now loop thru them all and add them to
    * my list. Also, process any special boundary conditions (periodic, ghost etc)
    */
   particleIndicies = particleSizesRecv/sizeof(Cparticle);

   for (Array<int,NDIM>::iterator ap = globals.procNeighbrs.begin();ap != globals.procNeighbrs.end();ap++) {
      vectInt coords = ap.position();
      int neighbr = globals.procNeighbrs(coords);
      if (neighbr >= 0) {
         double maxh = hmax;

         for (int z=0;z<particleIndicies(coords);z++) {
            Cparticle *recvP = &(particleBuffersRecv(coords)[z]);
            ps.push_back(*recvP);
            Cparticle *p = &(ps.back());
            for (int i=0;i<NDIM;i++) {
               if ((PERIODIC[i])&&(coords[i]==0)&&(RMIN[i] >= globals.procDomain[i*2]-PSEP)&&(RMIN[i] <= globals.procDomain[i*2]+PSEP)) {
                  p->r[i] = RMIN[i]-RMAX[i]+p->r[i]; 
                  if (p->iam==sph)
                     p->dens += DENS_DROP[i];
                  #ifdef SLK
                  p->currR[i] = RMIN[i]-RMAX[i]+p->currR[i];
                  #endif
               }
               if ((PERIODIC[i])&&(coords[i]==2)&&(RMAX[i] >= globals.procDomain[i*2+1]-PSEP)&&(RMAX[i] <= globals.procDomain[i*2+1]+PSEP)) {
                  p->r[i] = RMAX[i]-RMIN[i]+p->r[i]; 
                  if (p->iam==sph)
                     p->dens -= DENS_DROP[i];
                  #ifdef SLK
                  p->currR[i] = RMAX[i]-RMIN[i]+p->currR[i];
                  #endif
               }
               if ((GHOST[2*i])&&(coords[i]==0)&&(RMIN[i] >= globals.procDomain[i*2]-PSEP)&&(RMIN[i] <= globals.procDomain[i*2]+PSEP)) {
                  p->r[i] = 2.0*RMIN[i]-p->r[i];
                  if (GHOST[2*i]==1) {
                     p->v[i] = -p->v[i];
                     p->vhat[i] = -p->vhat[i];
                  } else {
                     p->v = -p->v;
                     p->vhat = -p->vhat;
                  }
               }
               if ((GHOST[2*i+1])&&(coords[i]==2)&&(RMAX[i] >= globals.procDomain[i*2+1]-PSEP)&&(RMAX[i] <= globals.procDomain[i*2+1]+PSEP)) {
                  p->r[i] = 2.0*RMAX[i]-p->r[i];
                  if (GHOST[2*i+1]==1) {
                     p->v[i] = -p->v[i];
                     p->vhat[i] = -p->vhat[i];
                  } else {
                     p->v = -p->v;
                     p->vhat = -p->vhat;
                  }
               }
            }
            CpInfo newpInfo;
            newpInfo.p = p;
            newpInfo.key = calcKey(newpInfo.p->r);
            pInfos.push_back(newpInfo);
            CpInfoLink newpInfoLink;
            newpInfoLink.ppInfo = &(pInfos.back());
            pInfoLinks.push_back(newpInfoLink);
            pInfos.back().ppInfoLink = &(pInfoLinks.back());
         }
      }
   }

   /*
    * now loop thru all the particles again and find all the ghost particles
    * that need to be sent to neighbouring cpus
    */
   for (vector<CpInfoLink>::iterator ppInfoLink = pInfoLinks.begin();ppInfoLink!=pInfoLinks.end();ppInfoLink++) {
      CpInfo *ppInfo = ppInfoLink->ppInfo;
      ppInfo->neighbrs.clear();
      Cparticle *thisP = ppInfo->p;
      vectInt inCoords = static_cast<vectInt>((thisP->r-dgmin)/dgspace+1);
      for (int i=0;i<NDIM;i++) {
         if (inCoords[i]>2) inCoords[i]=2;
         if (inCoords[i]<0) inCoords[i]=0;
      }
      if (any(inCoords!=1)) {
         vectInt lowerBounds;
         vectInt upperBounds;
         for (int i=0;i<NDIM;i++) {
            if (inCoords[i]<=1) {
               upperBounds[i] = 1;
               lowerBounds[i] = inCoords[i];
            } else {
               upperBounds[i] = inCoords[i];
               lowerBounds[i] = 1;
            }
         }
         vectInt extent = upperBounds-lowerBounds+1;
         Array<int,NDIM> dummy(lowerBounds,extent);

         for (Array<int,NDIM>::iterator pD = dummy.begin();pD !=dummy.end();pD++) {
            vectInt dCoords = pD.position();
            if (globals.procNeighbrs(dCoords) >= 0) {

               //add to ghost particle array
               if (ghostIndicies(dCoords)>=gBufferSizes(dCoords)) {
                  cerr << "ghostBuffers are full, exiting...."<<endl;
                  exit(-1);
               }
               ghostBuffersSend(dCoords)[ghostIndicies(dCoords)] = *thisP;
               ghostedParticles(dCoords)[ghostIndicies(dCoords)] = thisP;
               ghostIndicies(dCoords)++;
            }
         }
      }
   }
   sendSizesGhosts = ghostIndicies;

   /*
    * ok. We have found all the ghost and put them in ghostBuffersSend. Now
    * send the data to the correct neighbouring cpus
    */
   Array<int,NDIM> ghostSizesSend(3);
   Array<int,NDIM> ghostSizesRecv(3);
   ghostSizesSend = ghostIndicies*sizeof(CghostData);

   gettimeofday(&t1,NULL);
   if (globals.mpiRank%2==0) {
      int upper = int(pow(3.0,NDIM));
      for (int i=0;i<upper;i++) {
         vectInt split = 3;
         vectInt coords = Nmisc::numToCoords(i,split);
         if (globals.procNeighbrs(coords)>=0) {
            sendRecvGhosts(coords,i,requestSendRecv+i*4,&(ghostSizesSend(coords)),&(ghostSizesRecv(coords))); 
         }
         else {
            for (int j=0;j<4;j++) {
               requestSendRecv[i*4+j] = MPI_REQUEST_NULL;
            }
         }
      }
   } else {
      for (int i=int(pow(3.0,NDIM))-1;i>=0;i--) {
         vectInt split = 3;
         vectInt coords = Nmisc::numToCoords(i,split);
         if (globals.procNeighbrs(coords)>=0) {
            sendRecvGhosts(coords,i,requestSendRecv+i*4,&(ghostSizesSend(coords)),&(ghostSizesRecv(coords))); 
         }
         else {
            for (int j=0;j<4;j++) {
               requestSendRecv[i*4+j] = MPI_REQUEST_NULL;
            }
         }
      } 
   }
   MPI_Waitall(sendRecvArraySize, requestSendRecv, statusSendRecv);   

   gettimeofday(&t2,NULL);
   globals.wtTotalMPI.tv_sec += t2.tv_sec-t1.tv_sec;
   globals.wtTotalMPI.tv_usec += t2.tv_usec-t1.tv_usec;
      
   /*
    * all ghost data is now sent and received. Now loop thru the ghost
    * particles obtained from each neighbour and add them to procGhostParticles, At
    * the same time handle and special boundary conditions.
    */
   ghostIndicies = ghostSizesRecv/sizeof(CghostData);
   recvSizesGhosts = ghostIndicies;

   for (Array<int,NDIM>::iterator ap = globals.procNeighbrs.begin();ap != globals.procNeighbrs.end();ap++) {
      vectInt coords = ap.position();
      int neighbr = globals.procNeighbrs(coords);
      if (neighbr >= 0) {
         //work out ghost max h's
         double maxh = hmax;

         for (int z=0;z<ghostIndicies(coords);z++) {
            Cparticle p;
            p = ghostBuffersRecv(coords)[z];
            for (int i=0;i<NDIM;i++) {
               if ((PERIODIC[i])&&(coords[i]==0)&&(RMIN[i] >= globals.procDomain[i*2]-PSEP)&&(RMIN[i] <= globals.procDomain[i*2]+PSEP)) {
                  p.r[i] = RMIN[i]-RMAX[i]+p.r[i]; 
                  if (p.iam==sph)
                     p.dens += DENS_DROP[i];
                  #ifdef SLK
                  p.currR[i] = RMIN[i]-RMAX[i]+p.currR[i];
                  #endif
               }
               if ((PERIODIC[i])&&(coords[i]==2)&&(RMAX[i] >= globals.procDomain[i*2+1]-PSEP)&&(RMAX[i] <= globals.procDomain[i*2+1]+PSEP)) {
                  p.r[i] = RMAX[i]-RMIN[i]+p.r[i]; 
                  if (p.iam==sph)
                     p.dens -= DENS_DROP[i];
                  #ifdef SLK
                  p.currR[i] = RMAX[i]-RMIN[i]+p.currR[i];
                  #endif
               }
               if ((GHOST[2*i])&&(coords[i]==0)&&(RMIN[i] >= globals.procDomain[i*2]-PSEP)&&(RMIN[i] <= globals.procDomain[i*2]+PSEP)) {
                  p.r[i] = 2.0*RMIN[i]-p.r[i];
                  if (GHOST[2*i]==1) {
                     p.v[i] = -p.v[i];
                     p.vhat[i] = -p.vhat[i];
                  } else {
                     p.v = -p.v;
                     p.vhat = -p.vhat;
                  }
               }
               if ((GHOST[2*i+1])&&(coords[i]==2)&&(RMAX[i] >= globals.procDomain[i*2+1]-PSEP)&&(RMAX[i] <= globals.procDomain[i*2+1]+PSEP)) {
                  p.r[i] = 2.0*RMAX[i]-p.r[i];
                  if (GHOST[2*i+1]==1) {
                     p.v[i] = -p.v[i];
                     p.vhat[i] = -p.vhat[i];
                  } else {
                     p.v = -p.v;
                     p.vhat = -p.vhat;
                  }
               }
            }
            if (p.h>maxh) maxh = p.h;
            procGhostParticles.push_back(p);
         }

         //store max h for this neighbour
         procGhostMax2H(coords) = KERNAL_RADIUS*maxh;
      }
   }

   gettimeofday(&t22,NULL);
   globals.wtTotalUpdateDomain.tv_sec += t22.tv_sec-t11.tv_sec;
   globals.wtTotalUpdateDomain.tv_usec += t22.tv_usec-t11.tv_usec;
}

/*
 * non-blocking mpi send of the particle data in the buffer
 * particleBuffersSend. The function sends the number of particles to be sent before sending the data
 */
void CdataLL::sendRecvParticles(const vectInt coords,const int num,MPI_Request *request,int *particleSizesSend,int *particleSizesRecv) {
   int neighbr = globals.procNeighbrs(coords);
   vectInt middle = 1;
   vectInt split = 3;
   vectInt oppCoords = 2-coords;
   for (int i=0;i<NDIM;i++) {
      if (((coords[i]==2)&&(GHOST[2*i+1])&&(RMAX[i] >= globals.procDomain[i*2+1]-PSEP)&&(RMAX[i] <= globals.procDomain[i*2+1]+PSEP))||
          ((coords[i]==0)&&(GHOST[2*i])&&((RMIN[i] >= globals.procDomain[i*2]-PSEP)&&(RMIN[i] <= globals.procDomain[i*2]+PSEP)))) { 
         oppCoords[i] = coords[i];
      }
   }
   int oppNum = Nmisc::coordsToNum(oppCoords,split);
        
   MPI_Status status;
   MPI_Isend(particleSizesSend,1,MPI_INT,neighbr,300+num,MPI_COMM_WORLD,request);
   MPI_Irecv(particleSizesRecv,1,MPI_INT,neighbr,300+oppNum,MPI_COMM_WORLD,request+1);
   MPI_Isend(particleBuffersSend(coords),*particleSizesSend,MPI_BYTE,neighbr,400+num,MPI_COMM_WORLD,request+2);
   MPI_Irecv(particleBuffersRecv(coords),pBufferSizes(coords)*sizeof(Cparticle),MPI_BYTE,neighbr,400+oppNum,MPI_COMM_WORLD,request+3);
}

/*
 * non-blocking mpi send of the ghost data in the buffer
 * ghostBuffersSend. The function sends the number of particles to be sent before sending the data
 */
void CdataLL::sendRecvGhosts(const vectInt coords,const int num,MPI_Request *request,int *ghostSizesSend,int *ghostSizesRecv) {
   int neighbr = globals.procNeighbrs(coords);
   vectInt middle = 1;
   vectInt split = 3;
   vectInt oppCoords = 2-coords;
      for (int i=0;i<NDIM;i++) {
      if (((coords[i]==2)&&(GHOST[2*i+1])&&(RMAX[i] >= globals.procDomain[i*2+1]-PSEP)&&(RMAX[i] <= globals.procDomain[i*2+1]+PSEP))||
          ((coords[i]==0)&&(GHOST[2*i])&&((RMIN[i] >= globals.procDomain[i*2]-PSEP)&&(RMIN[i] <= globals.procDomain[i*2]+PSEP)))) { 
            oppCoords[i] = coords[i];
         }
      }

   int oppNum = Nmisc::coordsToNum(oppCoords,split);
        
   MPI_Status status;
   MPI_Isend(ghostSizesSend,1,MPI_INT,neighbr,100+num,MPI_COMM_WORLD,request);
   MPI_Irecv(ghostSizesRecv,1,MPI_INT,neighbr,100+oppNum,MPI_COMM_WORLD,request+1);
   MPI_Isend(ghostBuffersSend(coords),*ghostSizesSend,MPI_BYTE,neighbr,200+num,MPI_COMM_WORLD,request+2);
   MPI_Irecv(ghostBuffersRecv(coords),gBufferSizes(coords)*sizeof(CghostData),MPI_BYTE,neighbr,200+oppNum,MPI_COMM_WORLD,request+3);
}

/*
 * for the user-initiated syncing, the number of ghost particles is already
 * know. So no need to send sizes
 */
void CdataLL::sendRecvDataSync(const vectInt coords,const int num,MPI_Request *request,int sendSize, void *sendBuffer, int recvSize, void *recvBuffer) {
   int neighbr = globals.procNeighbrs(coords);
   vectInt middle = 1;
   vectInt split = 3;
   vectInt oppCoords = 2-coords;
      for (int i=0;i<NDIM;i++) {
      if (((coords[i]==2)&&(GHOST[2*i+1])&&(RMAX[i] >= globals.procDomain[i*2+1]-PSEP)&&(RMAX[i] <= globals.procDomain[i*2+1]+PSEP))||
          ((coords[i]==0)&&(GHOST[2*i])&&((RMIN[i] >= globals.procDomain[i*2]-PSEP)&&(RMIN[i] <= globals.procDomain[i*2]+PSEP)))) { 
            oppCoords[i] = coords[i];
         }
      }
   int oppNum = Nmisc::coordsToNum(oppCoords,split);
        
   MPI_Status status;
   MPI_Isend(sendBuffer,sendSize,MPI_BYTE,neighbr,500+num,MPI_COMM_WORLD,request);
   MPI_Irecv(recvBuffer,recvSize,MPI_BYTE,neighbr,500+oppNum,MPI_COMM_WORLD,request+1);
}


void CdataLL::setGlobalTimestep(double dt) {
   double sendbuf = dt;
   double recvbuf[globals.mpiSize];
   MPI_Allgather(&sendbuf, 1, MPI_DOUBLE, recvbuf, 1, MPI_DOUBLE, MPI_COMM_WORLD); 
   globals.dt = dt;
   for (int i=0;i<globals.mpiSize;i++) {
      if (recvbuf[i]<globals.dt) globals.dt = recvbuf[i];
   }
}


/*
 * TODO: should make a general one: "function" over procs
 * TODO: now have two functions for two types (vect and double). 
 * template this to one function
 */
void CdataLL::sumOverProcs(vect *data,int num) {
   
   int newNum = num*sizeof(vect);
   vect *sendbuf = data;
   vect recvbuf[globals.mpiSize*num];

   MPI_Gather(sendbuf,newNum,MPI_BYTE,recvbuf,newNum,MPI_BYTE,0,MPI_COMM_WORLD);
   if (globals.mpiRank==0) {
      for (int i=1;i<globals.mpiSize;i++) {
         for (int j=0;j<num;j++) {
            recvbuf[j] += recvbuf[i*num+j];
         }
      }
      for (int i=1;i<globals.mpiSize;i++) {
         for (int j=0;j<num;j++) {
            recvbuf[i*num+j] = recvbuf[j];
         }
      }
   }
   MPI_Scatter(recvbuf,newNum,MPI_BYTE,sendbuf,newNum,MPI_BYTE,0,MPI_COMM_WORLD);
}

void CdataLL::sumOverProcs(double *data,int num) {
   
   int newNum = num*sizeof(double);
   double *sendbuf = data;
   double recvbuf[globals.mpiSize*num];

   MPI_Gather(sendbuf,newNum,MPI_BYTE,recvbuf,newNum,MPI_BYTE,0,MPI_COMM_WORLD);
   if (globals.mpiRank==0) {
      for (int i=1;i<globals.mpiSize;i++) {
         for (int j=0;j<num;j++) {
            recvbuf[j] += recvbuf[i*num+j];
         }
      }
      for (int i=1;i<globals.mpiSize;i++) {
         for (int j=0;j<num;j++) {
            recvbuf[i*num+j] = recvbuf[j];
         }
      }
   }
   MPI_Scatter(recvbuf,newNum,MPI_BYTE,sendbuf,newNum,MPI_BYTE,0,MPI_COMM_WORLD);
}


unsigned long CdataLL::calcKey(vect r) {
   unsigned int rint[NDIM];
   unsigned long key = 0;
   for (int i=0;i<NDIM;i++) {
      rint[i] = (unsigned int)floor((r[i]-rmin[i])*KEY_DIM_MAX/(rmax[i]-rmin[i]));
   }
   unsigned int mask = 0x01<<(KEY_DIM_BITS-1);
   for (int i=0;i<KEY_DIM_BITS;i++) {
      for (int j=0;j<NDIM;j++) {
         key <<= 1;
         key |= (rint[j] & mask) >> (KEY_DIM_BITS-1-i);
      }
      mask >>= 1;
   } 
   return key; 
}

void CdataLL::updateParticles() {
   initPinfos();
   checkPinfos();
}

/*
 * this function clears the particle information contained in pInfos and
 * pInfoLinks and re-generates it from the particle list (ps)
 */
void CdataLL::initPinfos() {
   pInfos.clear();
   pInfoLinks.clear();
   pInfos.reserve(ps.capacity());
   pInfoLinks.reserve(ps.capacity());
   for (particleContainer::iterator thisP = ps.begin();thisP!=ps.end();thisP++) {
      CpInfo newpInfo;
      newpInfo.p = &(*thisP);
      newpInfo.key = calcKey(thisP->r);
      pInfos.push_back(newpInfo);

      CpInfoLink newpInfoLink;
      newpInfoLink.ppInfo = &(pInfos.back());
      pInfoLinks.push_back(newpInfoLink);
      pInfos.back().ppInfoLink = &(pInfoLinks.back());

      vector<CpInfo>::iterator ppInfo = pInfos.end()-1;
      vector<CpInfoLink>::iterator ppInfoLink = pInfoLinks.end()-1;
   }
}

/*
 * Sort pInfos using insertion sort. Make sure that the sorting
 * maintains consistency with pInfoLink
 */
void CdataLL::sortPInfos() {
   for (vector<CpInfo>::iterator i=pInfos.begin()+1;i!=pInfos.end();i++) {
      for (vector<CpInfo>::iterator j=pInfos.begin();j!=i;j++) {
         if (j->key > i->key) {
            CpInfo tmpInfo = *j;
            *j = *i;
            j->ppInfoLink->ppInfo = &(*j);
            vector<CpInfo>::iterator jplus1 = j+1;
            for (vector<CpInfo>::iterator k=i;k!=jplus1;k--) {
               *k = *(k-1);
               k->ppInfoLink->ppInfo = &(*k);
            }
            *jplus1 = tmpInfo;
            jplus1->ppInfoLink->ppInfo = &(*jplus1);
         }
      }
   }
}

/*
 * for each member of pInfoLinks, the particle information pointed to by ppInfo
 * should have a class member ppInfoLink that points back to the original entry in pInfoLinks
 */
void CdataLL::checkPinfos() {
   particleContainer::iterator p = ps.begin();
   for (vector<CpInfoLink>::iterator ppInfoLink = pInfoLinks.begin();ppInfoLink!=pInfoLinks.end();ppInfoLink++) {
      CpInfo *ppInfo = ppInfoLink->ppInfo;
      if (ppInfo->ppInfoLink != &(*ppInfoLink)) {
         cerr <<"Error: CdataLL::initTraverseOrder: ppInfo and ppInfoLink are not consistant!"<<endl;
         exit(-1);
      }
      if (ppInfo->p != &(*p)) {
         cerr <<"Error: CdataLL::initTraverseOrder: ppInfo, ppInfoLink and p are not consistant!"<<endl;
         exit(-1);
      }
      p++;
   }
}


void CdataLL::initTraverseOrder() {
   initPinfos();   
   sortPInfos();
   checkPinfos();
}

void CdataLL::calcTraverseOrder() {
   particleContainer::iterator p = ps.begin();
   for (vector<CpInfoLink>::iterator ppInfoLink = pInfoLinks.begin();ppInfoLink!=pInfoLinks.end();ppInfoLink++) {
      CpInfo *ppInfo = ppInfoLink->ppInfo;
      if (ppInfo->ppInfoLink != &(*ppInfoLink)) {
         cerr <<"Error: CdataLL::calcTraverseOrder: ppInfo and ppInfoLink are not consistant!"<<endl;
         exit(-1);
      }
      if (ppInfo->p != &(*p)) {
         cerr <<"Error: CdataLL::calcTraverseOrder: ppInfo, ppInfoLink and p are not consistant!"<<endl;
         exit(-1);
      }
      ppInfo->key = calcKey(ppInfo->p->r);
      ppInfo->neighbrs.clear();
      p++;
   }
   sortPInfos();
}

void CdataLL::calcDomainLimits() {
   for (int i=0;i<NDIM;i++) {
      rmin[i] = globals.procDomain[i*2];
      rmax[i] = globals.procDomain[i*2+1];
   }
   hmax = -100000;
   for (particleContainer::iterator thisP = ps.begin();thisP!=ps.end();thisP++) {
      for (int j=0;j<NDIM;j++) {
         if (thisP->r[j] < rmin[j]) { rmin[j] = thisP->r[j]; }
         if (thisP->r[j] > rmax[j]) { rmax[j] = thisP->r[j]; }
      }
     
      if (thisP->h > hmax) hmax = thisP->h;
   }
   for (int i=0;i<NDIM;i++) {
      vectInt coords = 1;
      coords[i] = 0;
      if (procGhostMax2H(coords)<hmax) procGhostMax2H(coords)=KERNAL_RADIUS*hmax;
      coords[i] = 2;
      if (procGhostMax2H(coords)<hmax) procGhostMax2H(coords)=KERNAL_RADIUS*hmax;
   }
}


void CdataLL::insert(Cparticle &p) {
   vectInt cellI = getCellI(p.r);
   cells(cellI).push_back(&p);
}

void CdataLL::addIfNeighbr(Cparticle &_p,vector<Cparticle *> &_neighbrs,vector<Cparticle *> &_toAdd) {
   for (int i=0;i<_toAdd.size();i++) {
      if (_p.is_neighbr(*(_toAdd[i]))) {
         _neighbrs.push_back(_toAdd[i]);
      }
   }
}

/*
 * finds bucket of given particle. Loops over neighbouring buckets in order to
 * find neighbouring particles.
 */
void CdataLL::calcNeighbours(vector<Cparticle *> &_neighbrs,Cparticle &_p) {
   vectInt cellI = getCellI(_p.r);
   vectInt tcellI = cellI;
   for (int i=cellI[0]-1;i<=cellI[0]+1;i++) {
      tcellI[0] = i;
      if (NDIM > 1) {
         for (int j=cellI[1]-1;j<=cellI[1]+1;j++) {
            tcellI[1] = j;
            if (NDIM > 2) {
               for (int k=cellI[2]-1;k<=cellI[2]+1;k++) {
                  tcellI[2] = k;
                  addIfNeighbr(_p,_neighbrs,cells(tcellI));
               }
            } else {
               addIfNeighbr(_p,_neighbrs,cells(tcellI));
            }
         }
      } else {
         addIfNeighbr(_p,_neighbrs,cells(tcellI));
      } 
   }
}


Cparticle *CdataLL::getParticleFromTag(int tag) {
   if (tagSortedParticlePointers.empty()) {
      cout <<"data struct empty, filling tagSortedParticlePointers"<<endl;
      tagSortedParticlePointers.resize(ps.size());
      for (vector<Cparticle>::reverse_iterator i = ps.rbegin();i!=ps.rend();i++) {
         if (tagSortedParticlePointers.size()<i->tag) {
            tagSortedParticlePointers.resize(i->tag*2);
         }
         tagSortedParticlePointers[i->tag-1].tag = i->tag;
         tagSortedParticlePointers[i->tag-1].p = &(*i);
      }
   }
   if (tagSortedParticlePointers[tag-1].tag != tag) {
      cerr << "Error in CdataLL::getParticleFromTag(): the tag aint here!!!"<<endl;
      exit(-1);
   }
   return tagSortedParticlePointers[tag-1].p;
}

CpInfo *CdataLL::getPinfoFromTag(int tag) {
   if (tagSortedPinfoPointers.empty()) {
      cout <<"data struct empty, filling tagSortedPinfoPointers"<<endl;
      tagSortedPinfoPointers.resize(pInfos.size());
      for (vector<CpInfo>::reverse_iterator i = pInfos.rbegin();i!=pInfos.rend();i++) {
         if (tagSortedPinfoPointers.size()<i->p->tag) {
            tagSortedPinfoPointers.resize(i->p->tag*2);
         }
         tagSortedPinfoPointers[i->p->tag-1].tag = i->p->tag;
         tagSortedPinfoPointers[i->p->tag-1].pInfo = &(*i);
      }
   }
   if (tagSortedPinfoPointers[tag-1].tag != tag) {
      cerr << "Error in CdataLL::getPinfoFromTag(): the tag aint here!!!"<<endl;
      exit(-1);
   }
   return tagSortedPinfoPointers[tag-1].pInfo;
}

int CdataLL::calcBufferSize(vectInt coords) {
   double bufferSize = 1.0;
   double bufferSizeDem = 1.0;
#ifdef HACK_TO_SAVE_MEM
   if (coords[2]!=-1) {
      const double nnxx = pow((1-POROSITY)/DEM_VOL,1.0/NDIM);
      const double demSep = 1.0/nnxx;
      bufferSizeDem *= pow(2.0*DEM_RADIUS,NDIM)/pow(demSep,NDIM);
   }
#endif
   for (int i=0;i<NDIM;i++) {
      if (coords[i]==0) {
#ifdef LIQ_DEM
         bufferSizeDem *= (RMAX[i]-RMIN[i])/(2.0*DEM_RADIUS);
#endif
         bufferSize *= (RMAX[i]-RMIN[i])/PSEP;
      } else {
#ifdef LIQ_DEM
         bufferSizeDem *= 2.0*H/(2.0*DEM_RADIUS);
#endif
         bufferSize *= 2.0*H/PSEP;
      }
   }
#ifdef LIQ_DEM
   return 1.5*(bufferSize+min(bufferSizeDem,double(MAX_NUM_DEM_PARTICLES)));
#else
   return 1.5*bufferSize;
#endif
}    
int CdataLL::calcPBufferSize(vectInt coords) {
   return int(ceil(calcBufferSize(coords)/4.0)); 
}