ParallelGravity.h

Go to the documentation of this file.

#ifndef PARALLELGRAVITY_H
#define PARALLELGRAVITY_H

#include "config.h"
#include <string>
#include <map>
#include <vector>
#include <algorithm>

#include "pup_stl.h"
#include "ckio.h"

#include "Vector3D.h"
#include "tree_xdr.h"
#include "TipsyFile.h"
#include "SFC.h"
#include "TreeNode.h"
#include "GenericTreeNode.h"
#include "Interval.h"
#include "parameters.h"
#include "param.h"
#include "dumpframe.h"
#include <liveViz.h>

#include "TaggedVector3D.h"

#include "codes.h"
#include "CacheInterface.h"

#ifdef SPCUDA
#include "EwaldCUDA.h"
#endif

#ifdef CUDA
#include "cuda_typedef.h"
#endif

#include "keytype.h"

PUPbytes(InDumpFrame);
PUPbytes(COOL);
PUPbytes(COOLPARAM);

#ifdef HPM_COUNTER
#include <libhpm.h>
#endif

#include <map>

#define MERGE_REMOTE_REQUESTS_VERBOSE(X) /*CkPrintf x*/

inline void CkMustAssert(bool cond, const char *err)
{
    if (!cond) CkAbort("%s", err);
}

using namespace std;

using namespace Tree;

enum LBStrategy{
  Null=0,
  Multistep,
  Orb3d,
  Multistep_notopo,
  MultistepNode_notopo,
  Orb3d_notopo,
  MultistepOrb,
  HierarchOrb
};
PUPbytes(LBStrategy);

#ifdef SELECTIVE_TRACING
enum TraceState {
  TraceNormal = 0,
  TraceSkip
};
#endif

enum DomainsDec {
    SFC_dec=0,  // Space Filling Curve with Morton ordering
    Oct_dec=1,  // Oct tree
    ORB_dec=2,  // Bisect the longest axis, balancing particles
    SFC_peano_dec=3,    // SFC with Peano-Hilbert ordering
    SFC_peano_dec_3D=4, // Joachim Stadel's implementation of P-H ordering
    SFC_peano_dec_2D=5, // 2D version of Peano-Hilbert ordering
    ORB_space_dec=6     // Bisect space
};

enum NborDir {
  LEFT = 0,
  RIGHT
};
PUPbytes(NborDir);

const double ddTolerance = 0.1;

inline void operator|(PUP::er &p,DomainsDec &d) {
  int di;
  if (p.isUnpacking()) {
    p | di;
    d = (DomainsDec)di;
  } else {
    di = (int)d;
    p | di;
  }
}

#include "GravityParticle.h"

class SmoothParams;

class NewMaxOrder
{
 public:
    int64_t nMaxOrderGas;
    int64_t nMaxOrderDark;
    int64_t nMaxOrder;
    void pup(PUP::er& p) {
        p| nMaxOrderGas;
        p| nMaxOrderDark;
        p| nMaxOrder;
        }
    };

#include "InOutput.h"

#include "ParallelGravity.decl.h"

extern CProxy_Main mainChare;
extern int verbosity;
extern bool _cache;
extern int _nocache;
extern int _cacheLineDepth;
extern unsigned int _yieldPeriod;
extern DomainsDec domainDecomposition;
extern double dExtraStore;
extern double dMaxBalance;
extern double dFracLoadBalance;
extern double dGlassDamper;
extern int bUseCkLoopPar;
extern GenericTrees useTree;
extern CProxy_TreePiece treeProxy;
#ifdef REDUCTION_HELPER
extern CProxy_ReductionHelper reductionHelperProxy;
#endif
extern CProxy_LvArray lvProxy;      // Proxy for the liveViz array
extern CProxy_LvArray smoothProxy;  // Proxy for smooth reduction
extern CProxy_LvArray gravityProxy; // Proxy for gravity reduction
extern CProxy_TreePiece streamingProxy;
extern CProxy_DataManager dMProxy;
extern CProxy_IntraNodeLBManager nodeLBMgrProxy;
extern unsigned int numTreePieces;
extern unsigned int particlesPerChare;
extern int nIOProcessor;

extern CProxy_DumpFrameData dfDataProxy;
extern CProxy_PETreeMerger peTreeMergerProxy;
extern CProxy_CkCacheManager<KeyType> cacheGravPart;
extern CProxy_CkCacheManager<KeyType> cacheSmoothPart;
extern CProxy_CkCacheManager<KeyType> cacheNode;

extern CkGroupID dataManagerID;

extern int boundaryEvaluationUE;
extern int weightBalanceUE;
extern int networkProgressUE;
extern int nodeForceUE;
extern int partForceUE;

extern int tbFlushRequestsUE;
extern int prefetchDoneUE;

extern int _prefetch;
extern int _randChunks;
extern int _numChunks;
extern unsigned int bucketSize;

//jetley
extern int localNodesPerReq;
extern int remoteNodesPerReq;
extern int remoteResumeNodesPerReq;
extern int localPartsPerReq;
extern int remotePartsPerReq;
extern int remoteResumePartsPerReq;

extern double largePhaseThreshold;

extern cosmoType theta;
extern cosmoType thetaMono;

extern int numInitDecompBins;
extern int octRefineLevel;

class dummyMsg : public CMessage_dummyMsg{
public:
};

#if COSMO_STATS > 0
class TreePieceStatistics {
  u_int64_t nodesOpenedLocal;
  u_int64_t nodesOpenedRemote;
  u_int64_t nodeInterLocal;
  u_int64_t nodeInterRemote;
  u_int64_t particleInterLocal;
  u_int64_t particleInterRemote;
  u_int64_t openCriterionCalls;
  int nActive;

  TreePieceStatistics() : nodesOpenedLocal(0), nodesOpenedRemote(0),
      nodeInterLocal(0), nodeInterRemote(0), particleInterLocal(0),
      particleInterRemote(0), openCriterionCalls(0), nActive(0) { }

 public:
  TreePieceStatistics(u_int64_t nol, u_int64_t nor, u_int64_t occ, u_int64_t nil, u_int64_t nir,
              u_int64_t pil, u_int64_t pir, int na) :
      nodesOpenedLocal(nol), nodesOpenedRemote(nor), nodeInterLocal(nil),
      nodeInterRemote(nir), particleInterLocal(pil), particleInterRemote(pir),
      openCriterionCalls(occ), nActive(na) { }

  void printTo(CkOStream &os) {
    os << "  TreePiece: " << nActive << " particles active." << endl;
    os << "  TreePiece: " << nodesOpenedLocal << " local nodes opened, ";
    os << nodesOpenedRemote << " remote" << endl;
    os << "  TreePiece: " << openCriterionCalls << " num of open criterion calls" << endl;
    os << "  TreePiece: " << nodeInterLocal << " local particle-node interactions, ";
    os << nodeInterRemote << " remote" << endl;
    os << "  TreePiece: " << particleInterLocal << " local particle-particle interactions, ";
    os << particleInterRemote << " remote" << endl;
    os << "  TreePiece: "
       << (particleInterLocal + particleInterRemote)/(double) nActive
       << " particles, "
       << (nodeInterLocal + nodeInterRemote)/(double) nActive
       << " nodes per particle" << endl;
  }

  static CkReduction::reducerType sum;

  static CkReductionMsg *sumFn(int nMsg, CkReductionMsg **msgs) {
    TreePieceStatistics ret;
    for (int i=0; i<nMsg; ++i) {
      CkAssert(msgs[i]->getSize() == sizeof(TreePieceStatistics));
      TreePieceStatistics *data = (TreePieceStatistics *)msgs[i]->getData();
      ret.nodesOpenedLocal += data->nodesOpenedLocal;
      ret.nodesOpenedRemote += data->nodesOpenedRemote;
      ret.openCriterionCalls += data->openCriterionCalls;
      ret.nodeInterLocal += data->nodeInterLocal;
      ret.nodeInterRemote += data->nodeInterRemote;
      ret.particleInterLocal += data->particleInterLocal;
      ret.particleInterRemote += data->particleInterRemote;
      ret.nActive += data->nActive;
    }
    return CkReductionMsg::buildNew(sizeof(TreePieceStatistics), &ret);
  }
};
#endif

class ComputeChunkMsg : public CMessage_ComputeChunkMsg {
  ComputeChunkMsg() {} // not available
 public:
  int chunkNum;

  ComputeChunkMsg(int i) : chunkNum(i) {
  }
};

class ORBSplittersMsg : public CMessage_ORBSplittersMsg{
public:
  int length;
  double *pos;
  char *dim;
  CkCallback cb;

  ORBSplittersMsg(int len, CkCallback callback): length (len), cb(callback) {}

};

class ParticleShuffleMsg : public CMessage_ParticleShuffleMsg{
public:
    int nloads;
    int n;
    int nSPH;
    int nStar;
    double *loads;
    unsigned int *parts_per_phase;
    GravityParticle *particles;
    extraSPHData *pGas;
    extraStarData *pStar;
    ParticleShuffleMsg(int nload, int npart, int nsph, int nstar): 
      nloads(nload), n(npart), nSPH(nsph), nStar(nstar) {}
};

#ifdef PUSH_GRAVITY
#include "ckmulticast.h"

struct BucketMsg : public CkMcastBaseMsg, public CMessage_BucketMsg {
  GenericTreeNode *buckets;
  int numBuckets;
  ExternalGravityParticle *particles;
  int numParticles;
  int whichTreePiece;
};
#endif
    
class CountSetPart 
{
 public:
    int index;          /* chare index */
    int nAddGas;
    int nDelGas;
    int nAddDark;
    int nDelDark;
    int nAddStar;
    int nDelStar;

    void pup(PUP::er& p) {
    p | index;
    p | nAddGas;
    p | nDelGas;
    p | nAddDark;
    p | nDelDark;
    p | nAddStar;
    p | nDelStar;
    }
    };

/*
 * Multistepping routines
 *
 * Each major timestep can have MAXSUBSTEPS substeps where MAXSUBSTEPS
 * is a large power of 2 that can fit in an integer: 1 << MAXRUNG,
 * with MAXRUNG something like 30.
 * A given particle is on a "Rung" such that it takes 1 << Rung
 * substeps per major timestep.  That is, it's force is updated every
 * 1 << (MAXRUNG - Rung) smallest substeps.
 *
 * Needed routines:
 * DtToRung(): take an ideal timestep, the major timestep and
 * determine a rung.
 */

const int MAXRUNG = 30;
const int MAXSUBSTEPS = 1 << MAXRUNG;
const double MAXSUBSTEPS_INV = 1 / (double)MAXSUBSTEPS;

inline int RungToSubsteps(int iRung) {
  CkAssert(iRung <= MAXRUNG);
  return 1 << (MAXRUNG - iRung);
}

inline int DtToRung(double dDelta, double dTideal) {
  int iSteps = (int) ceil(dDelta/dTideal);

  int iRung = 0;
  iSteps--;
  while(iSteps > 0) {
    iRung++;
    iSteps >>= 1;
  }
  return iRung;
}

inline double RungToDt(double dDelta, int iRung) {
  return dDelta*RungToSubsteps(iRung)*MAXSUBSTEPS_INV;
}

const int PHASE_FEEDBACK = MAXRUNG + 1;


inline double PoverRhoFloorJeans(double dResolveJeans, GravityParticle *p)
{
    /*
     * Add pressure floor to keep Jeans Mass
     * resolved.  In comparison with Agertz et
     * al. 2009, dResolveJeans should be 3.0:
     * P_min = 3*G*max(h,eps)^2*rho^2
     * Note that G = 1 in our code
     */
#ifdef JEANSSOFTONLY
    double l2 = p->soft*p->soft;
#else
    double l2 = 0.25*p->fBall*p->fBall;
#ifdef JEANSSOFT
    double e2 = p->soft*p->soft; 
    if (l2 < e2) l2 = e2; /* Jeans scale can't be smaller than softening */
#endif
#endif
    return l2*dResolveJeans*p->fDensity;
}

const double GAMMA_JEANS = 2.0;
const double GAMMA_NONCOOL = 5.0/3.0;


class Main : public CBase_Main {
    CkArgMsg *args;
    std::string basefilename;
        OutputParams *pOutput;
    // NChilada file names used to generate the XML description
    CkVec<std::string> *NCgasNames;
    CkVec<std::string> *NCdarkNames;
    CkVec<std::string> *NCstarNames;
    CkCallback cbIO;
        Ck::IO::File fIOFile;
    CProxy_Sorter sorter;
    int64_t nTotalParticles;
    int64_t nTotalSPH;
    int64_t nTotalDark;
    int64_t nTotalStar;
    int64_t nSink;
    int64_t nMaxOrderGas;  /* Maximum iOrders */
    int64_t nMaxOrderDark;
    int64_t nMaxOrder;
    
    double dTime;       /* Simulation time */
    double dTime0;      
    double dEcosmo;     /* variables for integrating
                   Lazer-Irvine eq. */
    double dUOld;
    double dTimeOld;
    PRM prm;        /* parameter parsing info */
    Parameters param; /* actual parameters */
    CkVec<double> vdOutTime; // Desired output times
    int iOut;
    /*
     ** Tracking for frame dumping function
     */
    int bDumpFrame;
    struct DumpFrameContext **df;
    int bIsRestarting;
        int bHaveAlpha;
    int bChkFirst;      /* alternate between 0 and 1 for checkpoint */
    double dSimStartTime;   // Start time for entire simulation
    int iStop;      /* indicate we're stopping the
                   simulation early */
    int64_t nActiveGrav;
    int64_t nActiveSPH;

#ifdef CUDA
          double localNodesPerReqDouble;
          double remoteNodesPerReqDouble;
          double remoteResumeNodesPerReqDouble;
          double localPartsPerReqDouble;
          double remotePartsPerReqDouble;
          double remoteResumePartsPerReqDouble;
#endif

#ifdef CHECK_TIME_WITHIN_BIGSTEP
       double wallTimeStart;
#endif

       class timing_fields {
       public:
           int count;           
           double tGrav;        
           double tuDot;        
           double tDD;          
           double tLoadB;       
           double tTBuild;      
           double tAdjust;      
           double tEmergAdjust; 
           double tKick;        
           double tDrift;       
           double tCache;       
       public:
           void clear() {
               count = 0;
               tGrav = tuDot = tDD = tLoadB = tTBuild = tAdjust
                   = tEmergAdjust = tKick = tDrift = tCache = 0.0;
               }
           };
       
       CkVec<timing_fields> timings;  
       void writeTimings(int iStep);

#ifdef SELECTIVE_TRACING
       int monitorRung;
       int monitorStart;
       int numTraceIterations;
       int numSkipIterations;
       int numMaxTrace;
       int traceIteration;
       int traceState;
       bool projectionsOn;

       void turnProjectionsOn(int activeRung);
       void turnProjectionsOff();
#endif


public:

    Main(CkArgMsg* m);
    Main(CkMigrateMessage *m);

    void niceExit();
    void setupICs();
    void initialForces();
    void doSimulation();
    void restart(CkCheckpointStatusMsg *msg);
    void waitForGravity(const CkCallback &cb, double startTime,
            int activeRung);
        void advanceBigStep(int);
        void domainDecomp(int iPhase);
        void loadBalance(int iPhase);
        void buildTree(int iPhase);
        void startGravity(const CkCallback& cbGravity, int iActiveRung,
            double *startTime) ;
        void externalGravity(int iActiveRung);
        void updateuDot(int iActiveRung, const double duKick[],
            const double dStartTime[], int bUpdateState, int bAll);
        void kick(bool bClosing, int iActiveRung, int nextMaxRung,
            const CkCallback &cbGravity, double gravStartTime);
    int adjust(int iKickRung);
    void rungStats();
    void countActive(int activeRung);
        void emergencyAdjust(int iRung);
    void starCenterOfMass();
    void calcEnergy(double, double, const char *);
    void getStartTime();
    void getOutTimes();
    int bOutTime();
    void writeOutput(int iStep) ;
        void outputBinary(OutputParams& params, int bParaWrite,
                          const CkCallback& cb);
        void cbOpen(Ck::IO::FileReadyMsg *msg);
        void cbIOReady(Ck::IO::SessionReadyMsg *msg);
        void cbIOComplete(CkMessage *msg);
        void cbIOClosed(CkMessage *msg);
        std::string getNCNextOutput(OutputParams& params);
    void writeNCXML(std::string filename);
    void NCXMLattrib(ofstream *desc, CkVec<std::string> *names, std::string family);
    void updateSoft();
    void growMass(double dTime, double dDelta);
    void initSph();
    void initCooling();
    void initStarLog();
    int ReadASCII(char *extension, int nDataPerLine, double *dDataOut);
        void restartGas();
    void doSph(int activeRung, int bNeedDensity = 1);
    void AGORAfeedbackPreCheck(double dTime, double dDelta, double dTimeToSF);
    void FormStars(double dTime, double dDelta);
    void StellarFeedback(double dTime, double dDelta);
    void outputBlackHoles(double dTime);
    void SetSink();
    void FormSinks(double dTime, double dDelta, int iKickRung);
    void doSinks(double dTime, double dDelta, int iKickRung);
    int DumpFrameInit(double dTime, double dStep, int bRestart);
    void DumpFrame(double dTime, double dStep);
    int nextMaxRungIncDF(int nextMaxRung);
    void addDelParticles();
    void memoryStats();
    void memoryStatsCache();
    void pup(PUP::er& p);
    void liveVizImagePrep(liveVizRequestMsg *msg);
        void doSIDM(double dTime,double dDelta, int activeRung); /* SIDM */
};

/* IBM brain damage */
#undef hz
typedef struct ewaldTable {
  double hx,hy,hz;
  double hCfac,hSfac;
} EWT;

// jetley
class MissRecord;
class State;

typedef struct OffsetNodeStruct
{
      GenericTreeNode *node;
      int offsetID;
}OffsetNode;

#if INTERLIST_VER > 0
typedef CkQ<OffsetNode> CheckList;
typedef CkVec<OffsetNode> UndecidedList;
typedef CkVec<UndecidedList> UndecidedLists;
#endif

typedef struct particlesInfoR{
    ExternalGravityParticle* particles;
    int numParticles;
    Vector3D<cosmoType> offset;
#if defined CHANGA_REFACTOR_PRINT_INTERACTIONS || defined CHANGA_REFACTOR_WALKCHECK_INTERLIST || defined CUDA
    NodeKey key;
#endif
#if COSMO_DEBUG > 1
    GenericTreeNode *nd;
#endif
} RemotePartInfo;

typedef struct particlesInfoL{
    GravityParticle* particles;
    int numParticles;
    Vector3D<cosmoType> offset;
#if defined CHANGA_REFACTOR_PRINT_INTERACTIONS || defined CHANGA_REFACTOR_WALKCHECK_INTERLIST || defined CUDA
    NodeKey key;
#endif
#if COSMO_DEBUG > 1 || defined CUDA
    GenericTreeNode *nd;
#endif
} LocalPartInfo;

typedef struct LoopParDataStruct {
  CkVec<GenericTreeNode*> lowNodes;  
  CkVec<int> bucketids;              
  CkVec<int> chunkids;               
  CkVec<CkVec<OffsetNode> > clists;  
  CkVec<CkVec<RemotePartInfo> > rpilists;  
  CkVec<CkVec<LocalPartInfo> > lpilists;   
  TreePiece* tp;                           
} LoopParData;



#ifdef CUDA
struct BucketActiveInfo{
  int start;
  int size;
};
#endif

class SmoothCompute;
#include "Compute.h"

#if INTERLIST_VER > 0 && defined CUDA
template<typename T> class GenericList;
#endif

struct NonLocalMomentsClient {
  TreePiece *clientTreePiece;
  GenericTreeNode *clientNode;

  NonLocalMomentsClient() :
    clientTreePiece(NULL),
    clientNode(NULL)
  {}

  NonLocalMomentsClient(TreePiece *tp, GenericTreeNode *node) : 
    clientTreePiece(tp),
    clientNode(node)
  {}
};

struct NonLocalMomentsClientList {
  GenericTreeNode *targetNode;
  CkVec<NonLocalMomentsClient> clients;

  NonLocalMomentsClientList() : 
    targetNode(NULL)
  {}

  NonLocalMomentsClientList(GenericTreeNode *node) :
    targetNode(node)
  {}

  void addClient(const NonLocalMomentsClient &cli){
    clients.push_back(cli);
  }
};

class TreePiece : public CBase_TreePiece {
   // jetley
   friend class PrefetchCompute;
   friend class GravityCompute;
   friend class SmoothCompute;
   friend class KNearestSmoothCompute;
   friend class ReSmoothCompute;
   friend class MarkSmoothCompute;
   friend class ListCompute;
   friend class NearNeighborState;
   friend class ReNearNeighborState;
   friend class MarkNeighborState;
   friend class BottomUpTreeWalk;
#if INTERLIST_VER > 0 && defined CUDA
   friend class DataManager;
   template<typename T> friend class GenericList;
#endif

   friend class RemoteTreeBuilder; 
   friend class LocalTreeBuilder; 

   TreeWalk *sTopDown;
   TreeWalk *twSmooth;
#if INTERLIST_VER > 0
   TreeWalk *sInterListWalk;
   // clearable, used for resumed walks
   State *sInterListStateRemoteResume;
#endif
   Compute *sGravity, *sPrefetch;
   SmoothCompute *sSmooth;
   
   Opt *sLocal, *sRemote, *sPref;
   Opt *optSmooth;

   State *sPrefetchState;
   State *sLocalGravityState, *sRemoteGravityState, *sSmoothState;
   typedef std::map<KeyType, CkVec<int>* > SmPartRequestType;
   // buffer of requests for smoothParticles.
   SmPartRequestType smPartRequests;

   CkVec<ActiveWalk> activeWalks;
   int completedActiveWalks; // XXX this should be part of the gravity
                 // walk state.
   int nCacheAccesses; // keep track of outstanding cache accesses to
               // know when writebacks complete.  XXX this
               // should be part of the smooth state
   
   unsigned int nPrevActiveParts;
   std::vector<double> savedPhaseLoad;
   std::vector<unsigned int> savedPhaseParticle;
   std::vector<double> savedPhaseLoadTmp;
   std::vector<unsigned int> savedPhaseParticleTmp;

   int memWithCache, memPostCache;  // store memory usage.
   int nNodeCacheEntries, nPartCacheEntries;  // store memory usage.

#ifdef PUSH_GRAVITY
   bool doMerge;
   bool createdSpanningTree;
   CProxySection_TreePiece allTreePieceSection;
   CkVec<GravityParticle> foreignParticles;
   CkVec<double> foreignParticleAccelerations;

   map<int,CkSectionInfo> cookieJar;
  
   BucketMsg *createBucketMsg();
   void unpackBuckets(BucketMsg *, GenericTreeNode *&foreignBuckets, int &numForeignBuckets);
   void calculateForces(GenericTreeNode *foreignBuckets, int numForeignBuckets);

#endif

 public:

#ifdef PUSH_GRAVITY
  void startPushGravity(int am, double myTheta);
  void recvPushBuckets(BucketMsg *);
  void recvPushAccelerations(CkReductionMsg *);
#endif

#if COSMO_PRINT_BK > 1
  State *getSRemoteGravityState(){ return sRemoteGravityState; }
  State *getSLocalGravityState(){ return sLocalGravityState; }
#endif
  void memCacheStats(const CkCallback &cb);
  void addActiveWalk(int iAwi, TreeWalk *tw, Compute *c, Opt *o, State *s);

  void markWalkDone();
  void finishWalk();
  void markSmoothWalkDone();
  void finishSmoothWalk();

  int getIndex() {
    return thisIndex;
  }

  void addToNodeInterRemote(int chunk, int howmany){
    nodeInterRemote[chunk] += howmany;
  }

  void addToParticleInterRemote(int chunk, int howmany){
    particleInterRemote[chunk] += howmany;
  }

  void addToNodeInterLocal(int howmany){
    nodeInterLocal += howmany;
  }

  void addToParticleInterLocal(int howmany){
    particleInterLocal += howmany;
  }

  void initiatePrefetch(int chunk);
  void startRemoteChunk();

  int getNumParticles(){
    return myNumParticles;
  }

  GravityParticle *getParticles(){return myParticles;}


#ifdef CUDA
        // this variable holds the number of buckets active at
        // the start of an iteration
        // it is used to ascertain how many buckets still need to 
        // be processed via the stateReady function with regard
        // to their local and remote-no-resume walks. 
        // if all numActiveBuckets
        // have been processed but we still have leftover nodes/particles
        // in the list of interations to the sent to the gpu, we flush
        // the list
        int numActiveBuckets; 
        int myNumActiveParticles;
        // First and Last indices of GPU particle
        int FirstGPUParticleIndex;
        int LastGPUParticleIndex;
        int NumberOfGPUParticles;
        BucketActiveInfo *bucketActiveInfo;

        int getNumBuckets(){
            return numBuckets;
        }

        void callFreeRemoteChunkMemory(int chunk);

        int getActiveRung(){ return activeRung; }
#ifdef HAPI_INSTRUMENT_WRS
        int getInstrumentId(){ return instrumentId; }
#endif
        // returns either all particles or only active particles,
        // depending on fraction of active particles to their
        // total count.
        int getDMNumParticles(){
          if(largePhase()){
            return myNumParticles;
          }
          else{
            return myNumActiveParticles; 
          }
        }

        int getNumActiveParticles(){
          return myNumActiveParticles;
        }

        void calculateNumActiveParticles(){ 
          myNumActiveParticles = 0;
          for(int i = 1; i <= myNumParticles; i++){
            if(myParticles[i].rung >= activeRung){
              myNumActiveParticles++;
            }
          }
        }

        bool largePhase(){
          return (1.0*myNumActiveParticles/myNumParticles) >= largePhaseThreshold;
        }

        void getDMParticles(CompactPartData *fillArray, int &fillIndex){
          NumberOfGPUParticles = 0;
          FirstGPUParticleIndex = fillIndex;//This is for the GPU Ewald
          if(largePhase()){
            for(int b = 0; b < numBuckets; b++){
              GenericTreeNode *bucket = bucketList[b];
              int buckstart = bucket->firstParticle;
              int buckend = bucket->lastParticle;
              GravityParticle *buckparts = bucket->particlePointer;
              bucket->bucketArrayIndex = fillIndex;
              for(int i = buckstart; i <= buckend; i++){
                fillArray[fillIndex] = buckparts[i-buckstart];
                fillIndex++;
              }
            }
          }
          else{
            for(int b = 0; b < numBuckets; b++){
              GenericTreeNode *bucket = bucketList[b];
              if(bucket->rungs < activeRung){
                continue;
              }
              BucketActiveInfo *binfo = &(bucketActiveInfo[b]);
              
              int buckstart = bucket->firstParticle;
              int buckend = bucket->lastParticle;
              GravityParticle *buckparts = bucket->particlePointer;

              binfo->start = fillIndex;
              for(int i = buckstart; i <= buckend; i++){
                if(buckparts[i-buckstart].rung >= activeRung){
                  fillArray[fillIndex] = buckparts[i-buckstart];
                  fillIndex++;
                }
              }
              binfo->size = fillIndex-binfo->start;
            }
          }
          //This is for the GPU Ewald
          if(FirstGPUParticleIndex == fillIndex){
            //This means no particle is on GPU
            FirstGPUParticleIndex = -1;
            LastGPUParticleIndex = -1;
            NumberOfGPUParticles = 0;
          }
          else{
            LastGPUParticleIndex = fillIndex - 1;
            NumberOfGPUParticles = LastGPUParticleIndex - FirstGPUParticleIndex + 1;
          }
        }

        bool isActive(int partNum){
          return myParticles[partNum].rung >= activeRung;
        }

        void clearMarkedBuckets(CkVec<GenericTreeNode *> &markedBuckets);
        void clearMarkedBucketsAll();

#ifdef HAPI_TRACE
        long long localNodeInteractions;
        long long localPartInteractions;
        long long remoteNodeInteractions;
        long long remotePartInteractions;
        long long remoteResumeNodeInteractions;
        long long remoteResumePartInteractions;
#endif

#ifdef HAPI_INSTRUMENT_WRS
        int instrumentId;

        double localNodeListConstructionTime;
        double remoteNodeListConstructionTime;
        double remoteResumeNodeListConstructionTime;
        double localPartListConstructionTime;
        double remotePartListConstructionTime;
        double remoteResumePartListConstructionTime;
        
        int nLocalNodeReqs;
        int nRemoteNodeReqs;
        int nRemoteResumeNodeReqs;
        int nLocalPartReqs;
        int nRemotePartReqs;
        int nRemoteResumePartReqs;

#endif

#endif

        void continueStartRemoteChunk(int chunk);
#ifdef CUDA
        void updateParticles(intptr_t data, int partIndex);
#endif
        void continueWrapUp();

#if INTERLIST_VER > 0
        void getBucketsBeneathBounds(GenericTreeNode *&node, int &start, int &end);
        void updateBucketState(int start, int end, int n, int chunk, State *state);
        void updateUnfinishedBucketState(int start, int end, int n, int chunk, State *state);
#endif

#if defined CHANGA_REFACTOR_WALKCHECK || defined CHANGA_REFACTOR_WALKCHECK_INTERLIST
        void addToBucketChecklist(int bucketIndex, NodeKey k){
          bucketcheckList[bucketIndex].insert(k);
          if(bucketIndex == TEST_BUCKET && thisIndex == TEST_TP)
            CkPrintf("[%d] add %ld\n", thisIndex, k);
        }
#endif

#if INTERLIST_VER > 0
        GenericTreeNode *getStartAncestor(int current, int previous, GenericTreeNode *dflt);
#endif
    inline GenericTreeNode *keyToNode(const Tree::NodeKey k){
          NodeLookupType::iterator iter = nodeLookupTable.find(k);
          if (iter != nodeLookupTable.end()) return iter->second;
          else return NULL;
        }

        GenericTreeNode *getBucket(int i){
          return bucketList[i];
        }

    double dStartTime;

private:        
    // liveViz 
    liveVizRequestMsg * savedLiveVizMsg;

        LBStrategy foundLB;
        // jetley - saved first internal node
        Vector3D<float> savedCentroid;
        int iPrevRungLB;
        int iActiveRungLB;

    CkCallback callback;
    CkCallback cbGravity;
    CkCallback cbSmooth;
  CkCallback after_dd_callback;
    unsigned int myNumParticles;
    unsigned int numActiveParticles;
    GravityParticle* myParticles;
  int nbor_msgs_count_;
    int nStore;
        bool bBucketsInited;

  // Temporary location to hold the particles that have come from outside this
  // TreePiece. This is used in the case where we migrate the particles and
  // detect completion of migration using QD.
    std::vector<GravityParticle> myTmpShuffleParticle;
    std::vector<extraSPHData> myTmpShuffleSphParticle;
    std::vector<extraStarData> myTmpShuffleStarParticle;
  ParticleShuffleMsg* myShuffleMsg;
    int myTreeParticles;
 public:
    int64_t nTotalParticles;
    int64_t nTotalSPH;
    int64_t nTotalDark;
    int64_t nTotalStar;
 private:
    unsigned int myNumSPH;
    extraSPHData *mySPHParticles;
    int nStoreSPH;
    unsigned int myNumStar;
    extraStarData *myStarParticles;
    int nStoreStar;
    int64_t nMaxOrder;
        int64_t nStartRead;
    int myIOParticles;

    std::vector<GenericTreeNode *> bucketList;

    std::vector<GravityParticle> mySortedParticles;
    std::vector<extraSPHData> mySortedParticlesSPH;
    std::vector<extraStarData> mySortedParticlesStar;
    std::vector<ParticleShuffleMsg*> incomingParticlesMsg;
        int incomingParticlesArrived;
        bool incomingParticlesSelf;

    double piecemass;

    int cnt;
    int packed;

          // the index assigned by the CacheManager upon registration
          int localIndex;

    //unsigned int numSplitters;
        //SFC::Key* splitters;
    CProxy_TreePiece pieces;
    //CProxy_DataManager dataManager; // unused...
    DataManager* dm;

#ifndef REDUCTION_HELPER
    CkVec<int64_t> myBinCounts;
#endif
    std::vector<int> myBinCountsORB;
    int myPlace;
    SFC::Key leftSplitter, rightSplitter;

    std::string basefilename;
    // Bounding box of the entire simulation
    OrientedBox<float> boundingBox;
    unsigned iterationNo;
    GenericTreeNode* root;
    NodePool *pTreeNodes;

    typedef std::map<NodeKey, CkVec<int>* >   MomentRequestType;
    MomentRequestType momentRequests;

    //double theta; -- moved as readonly
    //double thetaMono; -- moved as readonly
        int activeRung;

    int bPeriodic;
    int bComove;
    double dRhoFac;
    Vector3D<cosmoType> fPeriod;
    int nReplicas;
    int bEwald;     /* Perform Ewald */
    double fEwCut;
    double dEwhCut;
    EWT *ewt;
    int nMaxEwhLoop;
    int nEwhLoop;
#ifdef HEXADECAPOLE
    MOMC momcRoot;      /* complete moments of root */
#endif
        bool bEwaldInited;

    int bGasCooling;
#ifndef COOLING_NONE
    clDerivsData *CoolData;
#endif
        std::vector<int> iSeTab;
    int nSetupWriteStage;
    int64_t nStartWrite;    // Particle number at which this piece starts
                // to write file.

    NodeLookupType nodeLookupTable;

    //u_int64_t prefetchWaiting;
    Tree::NodeKey *prefetchRoots;
        OrientedBox<cosmoType> prefetchReq;

    int numChunks;

#if COSMO_STATS > 0
    u_int64_t myNumMACChecks;
    u_int64_t nodesOpenedLocal;
    u_int64_t nodesOpenedRemote;
    u_int64_t numOpenCriterionCalls;
#endif
    u_int64_t nodeInterLocal;
    u_int64_t *nodeInterRemote;
    u_int64_t particleInterLocal;
    u_int64_t *particleInterRemote;

    int nActive;        // number of particles that are active

    unsigned int numBuckets;
#if INTERLIST_VER > 0
    int prevRemoteBucket;
    int prevBucket;
#endif
    unsigned int ewaldCurrentBucket;
    BucketGravityRequest *bucketReqs;

        //CacheManager *localCache;

    // Entries used for the CacheManager
    EntryTypeGravityParticle gravityParticleEntry;
    EntryTypeSmoothParticle smoothParticleEntry;
    EntryTypeGravityNode gravityNodeEntry;

#if COSMO_DEBUG > 1 || defined CHANGA_REFACTOR_WALKCHECK || defined CHANGA_REFACTOR_WALKCHECK_INTERLIST
  typedef std::vector< std::multiset<Tree::NodeKey> > DebugList;
  DebugList bucketcheckList;
#endif


  bool firstTime;

  std::vector<int> tempBinCounts;

  std::list<GravityParticle *> orbBoundaries;

  int myExpectedCount;
  int myExpectedCountSPH;
  int myExpectedCountStar;

  unsigned int chunkRootLevel;

  OrientedBox<float>* boxes;
  char* splitDims;

  int phase;

  double myTotalMass;

 #if INTERLIST_VER > 0


 public:
#endif
  // called when a chunk has been used completely (chunkRemaining[chunk] == 0)
  void finishedChunk(int chunk);

  bool (*compFuncPtr[3])(GravityParticle,GravityParticle);
  double tmpTime;
  double totalTime;
 public:

#ifdef SPCUDA
  EwaldData *h_idata;
  CkCallback *cbEwaldGPU;
#endif
  void EwaldGPU(); 
  void EwaldGPUComplete();

#if COSMO_DEBUG > 1 || defined CHANGA_REFACTOR_WALKCHECK || defined CHANGA_REFACTOR_WALKCHECK_INTERLIST
  void checkWalkCorrectness();
  void combineKeys(Tree::NodeKey key,int bucket);
#endif

  /* DEBUGGING */
  void quiescence();
  /*END DEBUGGING */

  NodeLookupType &getNodeLookupTable() {
      return nodeLookupTable;
  }

  int getNumNodes() {
      return nodeLookupTable.size();
  }

  void deleteTree() {
    if(pTreeNodes != NULL) {
        delete pTreeNodes;
        pTreeNodes = NULL;
        root = NULL;
        nodeLookupTable.clear();
        }
    else {
        if (root != NULL) {
            root->fullyDelete();
            delete root;
            root = NULL;
            nodeLookupTable.clear();
            }
        }
    }

    void calculateRemoteMoments(GenericTreeNode* node);

    void checkTree(GenericTreeNode* node);

    bool nodeOwnership(const Tree::NodeKey nkey, int &firstOwner, int &lastOwner);


    void initBuckets();
    template <class Tsmooth>
    void initBucketsSmooth(Tsmooth tSmooth);
    void smoothNextBucket();
    void reSmoothNextBucket();
    void markSmoothNextBucket();
    void smoothBucketComputation();
    void startNextBucket();
    void doAllBuckets();
    void reconstructNodeLookup(GenericTreeNode *node);
    //void rebuildSFCTree(GenericTreeNode *node,GenericTreeNode *parent,int *);

public:
 TreePiece() : pieces(thisArrayID), root(0),
            iPrevRungLB (-1), sTopDown(0), sGravity(0),
            sPrefetch(0), sLocal(0), sRemote(0), sPref(0), sSmooth(0), 
            nPrevActiveParts(0) {
      dm = NULL;
      foundLB = Null; 
      iterationNo=0;
      usesAtSync = true;
      pTreeNodes = NULL;
      bucketReqs=NULL;
      nCacheAccesses = 0;
      memWithCache = 0;
      memPostCache = 0;
      nNodeCacheEntries = 0;
      nPartCacheEntries = 0;
      completedActiveWalks = 0;
      myPlace = -1;
      nSetupWriteStage = -1;
    //openingDiffCount=0;
    chunkRootLevel=0;
    //splitters = NULL;
#if COSMO_STATS > 0
      nodesOpenedLocal = 0;
      nodesOpenedRemote = 0;
      nodeInterLocal = 0;
      particleInterLocal = 0;

      numOpenCriterionCalls=0;

      piecemass = 0.0;
#endif
      packed=0;         /* output accelerations in x,y,z
                       packed format */
      cnt=0;            /* Counter over x,y,z when not packed */
      particleInterRemote = NULL;
      nodeInterRemote = NULL;

#if INTERLIST_VER > 0
      sInterListWalk = NULL;
#endif
#ifdef CUDA
          numActiveBuckets = -1;
#ifdef HAPI_TRACE
          localNodeInteractions = 0;
          localPartInteractions = 0;
          remoteNodeInteractions = 0;
          remotePartInteractions = 0;
          remoteResumeNodeInteractions = 0;
          remoteResumePartInteractions = 0;
#endif
#endif

      tmpTime=0.0;
      totalTime=0.0;
      // temporarely set to -1, it will updated after the tree is built
      numChunks=-1;
      prefetchRoots = NULL;
      ewt = NULL;
      nMaxEwhLoop = 100;
          bEwaldInited = false;

          incomingParticlesMsg.clear();
          incomingParticlesArrived = 0;
          incomingParticlesSelf = false;

          myParticles = NULL;
          mySPHParticles = NULL;
          myStarParticles = NULL;
      myNumParticles = myNumSPH = myNumStar = 0;
      nStore = nStoreSPH = nStoreStar = 0;
          bBucketsInited = false;
      myTreeParticles = -1;
      orbBoundaries.clear();
      boxes = NULL;
      splitDims = NULL;
      bGasCooling = 0;

#ifdef PUSH_GRAVITY
          createdSpanningTree = false;
#endif

          localTreeBuildComplete = false;
    }

    TreePiece(CkMigrateMessage* m): pieces(thisArrayID) {

      usesAtSync = true;
      //localCache = NULL;
      dm = NULL;
      bucketReqs = NULL;
      numChunks=-1;
      nCacheAccesses = 0;
      memWithCache = 0;
      memPostCache = 0;
      nNodeCacheEntries = 0;
      nPartCacheEntries = 0;
      completedActiveWalks = 0;
      prefetchRoots = NULL;
      //remaining Chunk = NULL;
          ewt = NULL;
          bEwaldInited = false;
      root = NULL;
      pTreeNodes = NULL;

      sTopDown = 0;
      sGravity = NULL;
      sPrefetch = NULL;
      sSmooth = NULL;
#if INTERLIST_VER > 0
      sInterListWalk = NULL;
#endif

          incomingParticlesMsg.clear();
          incomingParticlesArrived = 0;
          incomingParticlesSelf = false;

      cnt=0;
      nodeInterRemote = NULL;
          particleInterRemote = NULL;

      orbBoundaries.clear();
      boxes = NULL;
      splitDims = NULL;
          bBucketsInited = false;
      myTreeParticles = -1;


          localTreeBuildComplete = false;
    }

        private:
        void freeWalkObjects();
    void allocateStars() {
        nStoreStar = (int) (myNumStar*(1.0 + dExtraStore));
        // Stars tend to form out of gas, so make sure there is
        // enough room.
        nStoreStar += 12 + (int) (myNumSPH*dExtraStore);
        myStarParticles = new extraStarData[nStoreStar];
        }

        public:
    ~TreePiece() {
      if (verbosity>1) ckout <<"Deallocating treepiece "<<thisIndex<<endl;
      if(nStore > 0) delete[] myParticles;
      if(nStoreSPH > 0) delete[] mySPHParticles;
      if(nStoreStar > 0) delete[] myStarParticles;
      delete[] nodeInterRemote;
      delete[] particleInterRemote;
      delete[] bucketReqs;
          delete[] ewt;

      deleteTree();

      if(boxes!= NULL ) delete[] boxes;
      if(splitDims != NULL) delete[] splitDims;

#ifndef COOLING_NONE
      if(bGasCooling)
          CoolDerivsFinalize(CoolData);
#endif
          if (verbosity>1) ckout <<"Finished deallocation of treepiece "<<thisIndex<<endl;
    }

    void setPeriodic(int nReplicas, Vector3D<cosmoType> fPeriod, int bEwald,
             double fEwCut, double fEwhCut, int bPeriod,
                         int bComove, double dRhoFac);
    void BucketEwald(GenericTreeNode *req, int nReps,double fEwCut);
    void EwaldInit();
    void calculateEwald(dummyMsg *m);
  void calculateEwaldUsingCkLoop(dummyMsg *msg, int yield_num);
  void callBucketEwald(int id);
  void doParallelNextBucketWork(int id, LoopParData* lpdata);
    void initCoolingData(const CkCallback& cb);
    // Scale velocities (needed to convert to canonical momenta for
    // comoving coordinates.)
    void velScale(double dScale, const CkCallback& cb);

        void loadNChilada(const std::string& filename, const double dTuFac,
                          const CkCallback& cb);
        void readFloatBinary(OutputParams& params, int bParaRead,
            const CkCallback& cb);
    void loadTipsy(const std::string& filename, const double dTuFac,
                       const bool bDoublePos,
                       const bool bDoubleVel,
               const CkCallback& cb);
        void readTipsyArray(OutputParams& params, const CkCallback& cb);
        void resetMetals(const CkCallback& cb);
        void getMaxIOrds(const CkCallback& cb);
        void RestartEnergy(double dTuFac, const CkCallback& cb);
        void findTotalMass(const CkCallback &cb);
        void recvTotalMass(CkReductionMsg *msg);

    // Write a Tipsy file
    void writeTipsy(Tipsy::TipsyWriter& w,
            const double dvFac, // scale velocities
            const double duTFac, // convert temperature
                        const bool bDoublePos,
                        const bool bDoubleVel,
            const int bCool);
    // Find position in the file to start writing
    void setupWrite(int iStage, u_int64_t iPrevOffset,
            const std::string& filename, const double dTime,
            const double dvFac, const double duTFac,
                        const bool bDoublePos,
                        const bool bDoubleVel,
            const int bCool, const CkCallback& cb);
    // control parallelism in the write
    void parallelWrite(int iPass, const CkCallback& cb,
               const std::string& filename, const double dTime,
               const double dvFac, // scale velocities
               const double duTFac, // convert temperature
                           const bool bDoublePos,
                           const bool bDoubleVel,
               const int bCool);
    // serial output
    void serialWrite(u_int64_t iPrevOffset, const std::string& filename,
             const double dTime,
             const double dvFac, const double duTFac,
                         const bool bDoublePos,
                         const bool bDoubleVel,
             const int bCool, const CkCallback& cb);
    // setup for serial output
    void oneNodeWrite(int iIndex,
                   int iOutParticles,
                   int iOutSPH,
                   int iOutStar,
                   GravityParticle *particles, // particles to
                             // write
                   extraSPHData *pGas, // SPH data
                   extraStarData *pStar, // Star data
                   int *piSPH, // SPH data offsets
                   int *piStar, // Star data offsets
                   const u_int64_t iPrevOffset,
                   const std::string& filename,  // output file
                   const double dTime,      // time or expansion
                   const double dvFac,  // velocity conversion
                 const double duTFac, // temperature
                          // conversion
                           const bool bDoublePos,
                           const bool bDoubleVel,
              const int bCool,
              const CkCallback &cb);
    // Reorder for output
        void reOrder(int64_t nMaxOrder, const CkCallback& cb);
    // move particles around for output
    void ioShuffle(CkReductionMsg *msg);
    void ioAcceptSortedParticles(ParticleShuffleMsg *);
        void resetObjectLoad(const CkCallback& cb);
    // Assign keys after loading tipsy file and finding Bounding box
    void assignKeys(CkReductionMsg* m);
    void evaluateBoundaries(SFC::Key* keys, const int n, int isRefine, const CkCallback& cb);
    void unshuffleParticles(CkReductionMsg* m);
    void acceptSortedParticles(ParticleShuffleMsg *);
  void shuffleAfterQD();
  void unshuffleParticlesWoDD(const CkCallback& cb);
  void acceptSortedParticlesFromOther(ParticleShuffleMsg *);
  void setNumExpectedNeighborMsgs();

  /*****ORB Decomposition*******/
  void initORBPieces(const CkCallback& cb);
  void initBeforeORBSend(unsigned int myCount, unsigned int myCountGas,
             unsigned int myCountStar,
             const CkCallback& cb, const CkCallback& cback);
  void sendORBParticles();
  void acceptORBParticles(const GravityParticle* particles, const int n);
  void acceptORBParticles(const GravityParticle* particles, const int n,
              const extraSPHData *pGas, const int nGasIn,
              const extraStarData *pStar, const int nStarIn);
  void finalizeBoundaries(ORBSplittersMsg *splittersMsg);
  void evaluateParticleCounts(ORBSplittersMsg *splittersMsg);
  /*****************************/

  void kick(int iKickRung, double dDelta[MAXRUNG+1], int bClosing,
        int bNeedVPred, int bGasIsothermal, double dMaxEnergy, double duDelta[MAXRUNG+1],
        double gammam1, double dThermalCondSatCoeff,
        double dMultiPhaseMaxTime, double dMultiPhaseMinTemp, double dEvapCoeff, const CkCallback& cb);
  void drift(double dDelta, int bNeedVPred, int bGasIsothermal, double dvDelta,
             double duDelta, int nGrowMass, bool buildTree, double dMaxEnergy,
         const CkCallback& cb);
  void initAccel(int iKickRung, const CkCallback& cb);
#ifdef COOLING_MOLECULARH
  void distribLymanWerner(const CkCallback& cb);
#endif /*COOLING_MOLECULARH*/

  void applyFrameAcc(int iKickRung, Vector3D<double> frameAcc, const CkCallback& cb);
  void externalGravity(int activeRung, const ExternalGravity exGrav,
                       const CkCallback& cb);
  void adjust(int iKickRung, int bEpsAccStep, int bGravStep,
          int bSphStep, int bViscosityLimitdt,
          double dEta, double dEtaCourant, double dEtauDot,
              double dDiffCoeff, double dEtaDiffusion,
          double dDelta, double dAccFac,
          double dCosmoFac, double dhMinOverSoft,
              double dResolveJeans,
          int bDoGas,
          const CkCallback& cb);
  void truncateRung(int iCurrMaxRung, const CkCallback& cb);
  void rungStats(const CkCallback& cb);
  void countActive(int activeRung, const CkCallback& cb);
  void countType(int iType, const CkCallback& cb);
  void outputBlackHoles(const std::string& pszFileName, double dvFac,
                        long lFPos, const CkCallback &cb);
  void SetSink(double dSinkMassMin, const CkCallback &cb);
  void SinkStep(int iCurrSinkRung, int iKickRung, const CkCallback &cb);
  void formSinks(int bJeans, double dJConst2, int bDensity,
         double dDensityCut, double dTime, int iKickRung, int bSimple,
         const CkCallback &cb);
  void emergencyAdjust(int iRung, double dDelta, double dDeltaThresh,
               const CkCallback &cb);
  void assignDomain(const CkCallback& cb);
  void starCenterOfMass(const CkCallback& cb);
  void calcEnergy(const CkCallback& cb);
  void newParticle(GravityParticle *p);
  void adjustTreePointers(GenericTreeNode *node, GravityParticle *newParts);
  void colNParts(const CkCallback &cb);
  void newOrder(const NewMaxOrder *nStarts, const int n, const CkCallback &cb) ;
  
  void setNParts(int64_t _nTotalSPH, int64_t _nTotalDark,
         int64_t _nTotalStar, const CkCallback &cb);
  void setSoft(const double dSoft, const CkCallback &cb);
    void physicalSoft(const double dSoftMax, const double dFac,
              const int bSoftMaxMul, const CkCallback& cb);
    void growMass(int nGrowMass, double dDeltaM, const CkCallback& cb);
    void InitEnergy(double dTuFac, double z, double dTime, double gammam1,
            const CkCallback& cb);
    void updateuDot(int activeRung, double duDelta[MAXRUNG+1],
            double dStartTime[MAXRUNG+1], int bCool, int bAll,
            int bUpdateState, double gammam1, double dResolveJeans, const CkCallback& cb);
    void ballMax(int activeRung, double dFac, const CkCallback& cb);
    void sphViscosityLimiter(int bOn, int activeRung, const CkCallback& cb);
    void getAdiabaticGasPressure(double gamma, double gammam1, double dTuFac, double dThermalCondCoeff,
        double dThermalCond2Coeff, double dThermalCondSatCoeff, double dThermalCond2SatCoeff,
        double dEvapMinTemp, double dDtCourantFac, double dResolveJeans, const CkCallback &cb);
    void getCoolingGasPressure(double gamma, double gammam1, double dThermalCondCoeff,
        double dThermalCond2Coeff, double dThermalCondSatCoeff, double dThermalCond2SatCoeff,
        double dEvapMinTemp, double dDtCourantFac, double dResolveJeans, const CkCallback &cb);
#ifdef SPLITGAS
    void SplitGas(double dInitGasMass, const CkCallback& cb);
#endif
    inline COOL* Cool() {return dm->Cool;}
    void initRand(int iRand, const CkCallback &cb);
    void FormStars(Stfm param, double dTime, double dDelta, double dCosmoFac,
               const CkCallback& cb);
    void flushStarLog(const CkCallback& cb);
        void Feedback(const Fdbk &fb, double dTime, double dDelta,
                      const CkCallback& cb);
    void massMetalsEnergyCheck(int bPreDist, const CkCallback& cb);
    void SetTypeFromFileSweep(int iSetMask, char *file,
       struct SortStruct *ss, int nss, int *pniOrder, int *pnSet);
    void setTypeFromFile(int iSetMask, char *file, const CkCallback& cb);
    void getCOM(const CkCallback& cb, int bLiveViz);
    void getCOMByType(int iType, const CkCallback& cb, int bLiveViz);
    void DumpFrame(InDumpFrame in, const CkCallback& cb, int liveVizDump) ;
    void liveVizDumpFrameInit(liveVizRequestMsg *msg);
    void setProjections(int bOn);

#ifdef PUSH_GRAVITY
    void buildTree(int bucketSize, const CkCallback& cb, bool merge);
#else
    void buildTree(int bucketSize, const CkCallback& cb);
#endif

    void startOctTreeBuild(CkReductionMsg* m);
  void recvBoundary(SFC::Key key, NborDir dir);
    void recvdBoundaries(CkReductionMsg* m);

  /********ORB Tree**********/
  //void receiveBoundingBoxes(BoundingBoxes *msg);
  void startORBTreeBuild(CkReductionMsg* m);
  OrientedBox<float> constructBoundingBox(GenericTreeNode *node,int level, int numChild);
  void buildORBTree(GenericTreeNode * node, int level);
  /**************************/

    bool sendFillReqNodeWhenNull(CkCacheRequestMsg<KeyType> *msg);
    void requestRemoteMoments(const Tree::NodeKey key, int sender);
    void receiveRemoteMoments(const Tree::NodeKey key, Tree::NodeType type,
    int firstParticle, int numParticles, int remIdx,
    const MultipoleMoments& moments, const OrientedBox<double>& box,
            const OrientedBox<double>& boxBall,
            const unsigned int iParticleTypes, const int64_t nSPH);

    void calculateGravityLocal();
    void commenceCalculateGravityLocal();

    void calculateGravityRemote(ComputeChunkMsg *msg);
  void executeCkLoopParallelization(LoopParData *lpdata, int startbucket,
      int yield_num, int chunkNum, State* gravityState);
  int doBookKeepingForTargetActive(int curbucket, int end, int chunkNum, bool
    updatestate, State* gravityState);
  int doBookKeepingForTargetInactive(int chunkNum, bool updatestate,
    State* gravityState);

    void calculateSmoothLocal();
    void calculateReSmoothLocal();
    void calculateMarkSmoothLocal();
    void nextBucketSmooth(dummyMsg *msg);
    void nextBucketReSmooth(dummyMsg *msg);
    void nextBucketMarkSmooth(dummyMsg *msg);
#if INTERLIST_VER > 0
#if 0
  void calculateForceRemoteBucket(int bucketIndex, int chunk);
  void calculateForceLocalBucket(int bucketIndex);
#endif
#endif

  void startGravity(int am, double myTheta, const CkCallback& cb);
  void setupSmooth();
  void startSmooth(SmoothParams *p, int iLowhFix, int nSmooth,
           double dfBall2OverSoft2, const CkCallback &cb);
  void startReSmooth(SmoothParams *p, const CkCallback& cb);
  void startMarkSmooth(SmoothParams *p, const CkCallback& cb);

  void finishNodeCache(const CkCallback& cb);

    GenericTreeNode* requestNode(int remoteIndex, Tree::NodeKey lookupKey,
                                 int chunk, int reqID, int awi, void *source);
    void fillRequestNode(CkCacheRequestMsg<KeyType> *msg);
#if 0
    void receiveNode(GenericTreeNode &node, int chunk, unsigned int reqID);
#endif
    const GenericTreeNode* lookupNode(Tree::NodeKey key);
    const GravityParticle* lookupParticles(int begin);

    void finishBucket(int iBucket);

    void cachedWalkBucketTree(GenericTreeNode* node, int chunk, int reqID);

#if INTERLIST_VER > 0
    void calculateForcesNode(OffsetNode node, GenericTreeNode *myNode,
                 int level,int chunk);
    void calculateForces(OffsetNode node, GenericTreeNode *myNode,
                 int level,int chunk);
#endif

    ExternalGravityParticle *requestParticles(Tree::NodeKey key,int chunk,
                                              int remoteIndex,int begin,
                                              int end,int reqID, int awi,
                                              void *source);
    GravityParticle *requestSmoothParticles(Tree::NodeKey key, int chunk,
                                            int remoteIndex, int begin,int end,
                                            int reqID, int awi, void *source);
    void fillRequestParticles(CkCacheRequestMsg<KeyType> *msg);
    void fillRequestSmoothParticles(CkCacheRequestMsg<KeyType> *msg);
    void flushSmoothParticles(CkCacheFillMsg<KeyType> *msg);
    void processReqSmoothParticles();

    void getParticleInfoForLB(int64_t active_part, int64_t total_part);
        void startlb(const CkCallback &cb, int activeRung);
  void setTreePieceLoad(int activeRung);
  void populateSavedPhaseData(int phase, double tpload, unsigned int activeparts);
  bool havePhaseData(int phase);
  void savePhaseData(std::vector<double> &loads, std::vector<unsigned int>
    &parts_per_phase, double* shuffleloads, unsigned int *shuffleparts,
    int shufflelen);
    void ResumeFromSync();

    void outputASCII(OutputParams& params, int bParaWrite,
             const CkCallback& cb);
    void oneNodeOutVec(OutputParams& params, Vector3D<double>* avOut,
               int nPart, int iIndex, int bDone,
               const CkCallback& cb) ;
    void oneNodeOutArr(OutputParams& params, double* adOut,
               int nPart, int iIndex, int bDone,
               const CkCallback& cb) ;
    void oneNodeOutIntArr(OutputParams& params, int *aiOut,
                              int nPart, int iIndex, const CkCallback& cb);
        void outputBinaryStart(OutputParams& params, int64_t nStart,
                               const CkCallback& cb);
    void outputBinary(Ck::IO::Session, OutputParams& params);
        void minmaxNCOut(OutputParams& params, const CkCallback& cb);

    void outputStatistics(const CkCallback& cb);
        void collectStatistics(const CkCallback &cb);

        void nextBucket(dummyMsg *m);
  void nextBucketUsingCkLoop(dummyMsg *m);

    void report();
    void printTreeViz(GenericTreeNode* node, std::ostream& os);
    void printTree(GenericTreeNode* node, std::ostream& os);
    void pup(PUP::er& p);

        // jetley
        // need this in TreeWalk
        GenericTreeNode *getRoot() {return root;}
        // need this in Compute
    inline Vector3D<cosmoType> decodeOffset(int reqID) {
        int offsetcode = reqID >> 22;
        int x = (offsetcode & 0x7) - 3;
        int y = ((offsetcode >> 3) & 0x7) - 3;
        int z = ((offsetcode >> 6) & 0x7) - 3;

        Vector3D<cosmoType> offset(x*fPeriod.x, y*fPeriod.y, z*fPeriod.z);

        return offset;
        }

        void receiveNodeCallback(GenericTreeNode *node, int chunk, int reqID, int awi, void *source);
        void receiveParticlesCallback(ExternalGravityParticle *egp, int num, int chunk, int reqID, Tree::NodeKey &remoteBucket, int awi, void *source);
        void receiveParticlesFullCallback(GravityParticle *egp, int num, int chunk, int reqID, Tree::NodeKey &remoteBucket, int awi, void *source);
        void doAtSync();

        void balanceBeforeInitialForces(const CkCallback &cb);

        // For merging of remote moment requests
        // before sending messages during tree building
        public:
        void sendRequestForNonLocalMoments(GenericTreeNode *pickedNode);
        void mergeNonLocalRequestsDone();
        //void addTreeBuildMomentsClient(GenericTreeNode *targetNode, TreePiece *client, GenericTreeNode *clientNode);
        std::map<NodeKey,NonLocalMomentsClientList>::iterator createTreeBuildMomentsEntry(GenericTreeNode *pickedNode);


        private:
        // XXX - hashtable instead of map
        std::map<NodeKey,NonLocalMomentsClientList> nonLocalMomentsClients;
        bool localTreeBuildComplete;
        int getResponsibleIndex(int first, int last);
        

        GenericTreeNode *boundaryParentReady(GenericTreeNode *parent);
        void accumulateMomentsFromChild(GenericTreeNode *parent, GenericTreeNode *child);

        void deliverMomentsToClients(GenericTreeNode *);
        void deliverMomentsToClients(const std::map<NodeKey,NonLocalMomentsClientList>::iterator &it);
        void treeBuildComplete();
        void processRemoteRequestsForMoments();
        void sendParticlesDuringDD(bool withqd);
        void mergeAllParticlesAndSaveCentroid();
        bool otherIdlePesAvail();

};

class LvArray : public CBase_LvArray {
 public:
    LvArray() {}
    LvArray(CkMigrateMessage* m) {}
    } ;

int decodeReqID(int);
int encodeOffset(int reqID, int x, int y, int z);
bool bIsReplica(int reqID);
void printGenericTree(GenericTreeNode* node, std::ostream& os) ;
//bool compBucket(GenericTreeNode *ln,GenericTreeNode *rn);

#ifdef REDUCTION_HELPER

class TreePieceCounter : public CkLocIterator { 
  public:
    TreePieceCounter() { reset(); }
    void addLocation(CkLocation &loc);
    void reset();

  public:
    CkVec<TreePiece *> presentTreePieces;
};



class ReductionHelper : public CBase_ReductionHelper {
  public:
  ReductionHelper();
  ReductionHelper(CkMigrateMessage *);
  void pup(PUP::er &p);

  void countTreePieces(const CkCallback &cb);
  void reduceBinCounts(int nBins, int64_t *binCounts, const CkCallback &cb);
  void evaluateBoundaries(SFC::Key *keys, const int n, int isRefine, const CkCallback& cb);
  void evaluateBoundaries(const CkBitVector &binsToSplit, const CkCallback& cb);

  private:
  void senseLocalTreePieces();

  private:

  CkVec<int64_t> myBinCounts;
  int numTreePiecesCheckedIn;

  TreePieceCounter localTreePieces;
  std::vector<SFC::Key> splitters;
};
#endif

class NonEmptyTreePieceCounter : public CkLocIterator {            
  public:
    NonEmptyTreePieceCounter() { reset(); }
    void addLocation(CkLocation &loc);
    void reset();

  public:
    int count;
};

#endif //PARALLELGRAVITY_H