GridHierarchy.h

Go to the documentation of this file.

// -*- C++ -*-

#ifndef _included_GridHierarchy_h
#define _included_GridHierarchy_h

#include "DAGHParams.h"
#include "DAGHDefaults.h"

#include "DAGHDistribution.h"

#include "GridUnitList.h"
#include "GridBoxList.h"
#include "BBoxList.h"

#include "DCoords.h"

#include "CommServer.h"
#include "CommIOServer.h"

#include "GridFunctionVoid.h"
#include "AllocError.h"

#include <iosfwd>

// John Shalf's TCP++ Library"
#ifdef DAGH_IOV
#include <RawTCP.H>
#endif

/*************************************************************************/
/* Space-Time Dimension Specs */
/*************************************************************************/
struct dimspec {
  double low; 
  double high; 
  double h; 
  double hfine; 
  unsigned n; 
  unsigned nfine; 
  inline dimspec():low(0),high(0),h(0),hfine(0),n(0),nfine(0) {}
};
/*************************************************************************/

/*************************************************************************/
/* Merged bbox + grid unit list combo */
/*************************************************************************/
struct MergedGridUnit { 
  BBox bb;
  GridUnitList gul;
};
/*************************************************************************/

/*************************************************************************/
/* Utility routines */
/*************************************************************************/
inline double toWorld(struct dimspec const &ds, const int lc)
        { return (ds.low + (lc*ds.hfine)); }
inline double toWorldLower(struct dimspec const &ds, const int lc)
        { return (ds.high*((1.0*lc)/(1.0*ds.nfine)) + 
                  ds.low*(1.0 - (1.0*lc)/(1.0*ds.nfine))); }
inline double toWorldUpper(struct dimspec const &ds, const int lc, 
                           const int step, const int olap)
        { return (ds.high*((1.0*(lc+(1-olap)*step))/(1.0*ds.nfine)) + 
                  ds.low*(1.0 - (1.0*(lc+(1-olap)*step))/(1.0*ds.nfine))); }
inline int toLocal(struct dimspec const &ds, double const wc)
        { return ((int) ((wc - ds.low)/ds.hfine)); }
/*************************************************************************/

/*************************************************************************/
/* IO ping functions */
/*************************************************************************/
/* A PingFunction for Ending IO */
void DAGHIO_EndIOPingFunction(GridHierarchy &GH);
/*************************************************************************/

#ifdef DAGH_IOV
/*************************************************************************/
/* For Wei's IOV Server */
/*************************************************************************/
void WriteStream(GridHierarchy& GH, RawTCPport* Server, const GridBoxList& cgbl);
/*************************************************************************/
#endif

#ifndef DAGHDefaultDistribution
#define DAGHDefaultDistribution (DAGHCompositeDistribution)
#endif

#ifndef DAGHDefaultPartMinGUWidth
#define DAGHDefaultPartMinGUWidth       (1)
#endif
#ifndef DAGHDefaultRefineMinGUWidth
#define DAGHDefaultRefineMinGUWidth     (1)
#endif

#ifndef DAGHDefaultBoundary
#define DAGHDefaultBoundary     (DAGHBoundaryRegular)
#endif

#ifndef DAGHDefaultAdaptBoundary
#define DAGHDefaultAdaptBoundary (DAGHAdaptBoundaryInterp)
#endif

#ifndef DAGHDefaultRefineFactor
#define DAGHDefaultRefineFactor (2)
#endif

#ifndef DAGHDefaultRefineLevel
#define DAGHDefaultRefineLevel  (1)
#endif

#ifndef DAGHMaxIOTypes
#define DAGHMaxIOTypes          (5)
#endif 

#ifndef DAGHChkPtTagNameSize
#define DAGHChkPtTagNameSize    (256)
#endif

class GridHierarchy 
  {
   friend std::ostream& operator<<(std::ostream& os, const GridHierarchy& );
   friend std::ofstream& operator<<(std::ofstream& ofs, const GridHierarchy&);
   friend std::ifstream& operator>>(std::ifstream& ifs, GridHierarchy&);
   friend std::stringstream& operator<<(std::stringstream& ofs, 
       const GridHierarchy&);
   friend std::stringstream& operator>>(std::stringstream& ifs, 
       GridHierarchy&);

public:
   short const rank;

private:
   struct dimspec dspecs[DAGHMaxRank];  
   struct dimspec tspecs;
   BBox basebbox;

   short daghtype;
   short overlap[DAGHMaxRank];

   short bndrytype;
   short adaptbndrytype;
   short bndrywidth;
   short bndryperiodic[DAGHMaxRank];  
   short periodicidx[DAGHMaxRank];  

   short extghostwidth;

   /* Refinement Information */
   short levels;        
   short maxlev;        
   short crslev;
   short *refby;        

   short finelev;

   BBox intbbox;                
   BBoxList* intbaselist;    
   double* bndrycoords;          
   int nbndrycoords;

   /* Boxes of cut-off regions */
   BBox* cutbbox;              
   int ncutbbox;

   int *curtime;        

   int updatedstep;     

   /* number of buffer cells for proper nesting of higher levels */
   int nestbuffer;     

   /* coarsening factor  when constructing GridUnits */
   int cfactor;         
   int guwidth;

   /* Maximal Gridbox size */
   int maxgridsize;         

public:
   short distribution_type;
   DAGHDistribution* distribution;
   
public:
   /* Global Information */
   GridUnitList *complist;
   GridBoxList **globalgbl, **oldglobalgbl;

   /* Local Information */
   GridUnitList *locallist;     
   GridBoxList **localgbl, **oldlocalgbl;
   BBoxList** localbboxarray;

public:
   short gfnum;
   GridFunctionVoid **gflist;

   short io_type;
   DAGHIOServerRcv *io_write_server;
   DAGHIOServerSnd *io_read_server;
   DAGHIOServerPing *io_flush_server;
   DAGHIOServerPing *io_close_server;
   DAGHIOServerPing *io_end_server;

private:
   short chkptflag;
   int chkptpnum;
   char chkptname[DAGHChkPtTagNameSize];

#ifdef DAGH_IOV
private:
   RawTCPserver *ServerPort;
   RawTCPport *Server;
#endif

private:
   /*************************************************************************/
   /* Disable the following */
   /*************************************************************************/
   GridHierarchy(GridHierarchy const &other);
   GridHierarchy const &operator= (GridHierarchy const &other);

public:
   /*************************************************************************/
   /* Constructors */
   /*************************************************************************/
   GridHierarchy(const int ndim, const int type, const int max_level);

   /*************************************************************************/
   /* Destructor */
   /*************************************************************************/
   ~GridHierarchy(void);
    
   /*************************************************************************/
   /* Setup DAGH */
   /*************************************************************************/
public:
   void DAGH_SetBaseGrid(double const *bbox, const int *shape,
                         const int& ncuts, const int *cuts);
   void DAGH_SetTimeSpecs(double const starttime, 
                          double const stoptime, 
                          const int ntsteps);
   inline void DAGH_GetTimeSpecs(double& starttime, double& stoptime) {
     starttime = tspecs.low;
     stoptime  = tspecs.high;
   }

   inline void DAGH_SetBoundaryType(const int type) 
     { bndrytype = type; }
   inline void DAGH_SetAdaptBoundaryType(const int type) 
     { adaptbndrytype = type; }
   inline void DAGH_SetBoundaryWidth(const int width) 
     { bndrywidth = width; }
   inline void DAGH_SetPeriodicBoundaries(const int dir, const int value) 
     { if (dir>=0 && dir<rank) bndryperiodic[dir] = value; }

   inline void DAGH_SetExternalGhostWidth(const int width) 
     { extghostwidth = width; }

   void DAGH_SetRefineFactor(const int rfactor);
   void DAGH_SetRefineFactor(const int lev, const int rfactor);

   inline void DAGH_SetDistributionType(const int type)
        { distribution_type = type; if (distribution) distribution->set_dist(type); }
   inline void DAGH_SetDistributionType(const int type, BBoxList& bbl)
        { distribution_type = type; if (distribution) distribution->set_dist(type,bbl); }

   inline void DAGH_SetWorkFunction(void *wf)
        { if (distribution) distribution->set_workfunc(wf); }

   inline void DAGH_SetNestingBuffer(const int& nb)
        { nestbuffer = nb; }

   /*************************************************************************/
   /* Extend GridBoxLists in case of periodic boundaries */
   /*************************************************************************/
private:
   void DAGH_PeriodicAppendTogbl(GridBoxList *gbl, GridBox* gb, const BBox& bbox, 
                                 const BBox& bbtest, const int dir, const int base,
                                 const int maxlength);
public:
   void DAGH_PeriodicBoundaries(GridBoxList *gbl, GridBoxList *agbl);

   /*************************************************************************/
   /* Creation of GridUnitLists from BBoxLists */
   /*************************************************************************/
private:
   BBox shrinktoroot(BBox& bb);
   BBox growtoroot(BBox& bb);
   BBoxList* shrinktoroot(BBoxList& bbl);
   BBoxList* growtoroot(BBoxList& bbl);
   //what to count: cells (what==0) or workload (what!=0)
   int countcells(BBox& bb, BBoxList** bbla, int count);
public:
   void DAGH_CreateGridUnitList(GridUnitList*& gul, BBoxList** bbla, int count=DAGHFalse);
   void DAGH_CreateGridUnitList(GridUnitList*& gul, const BBox& _bb);
   void DAGH_CreateBBoxList(BBoxList*& bbl, GridUnitList& gul);
   
public:
   void DAGH_ComposeHierarchy(void);

   /*************************************************************************/
   /* Setup GridFunctions */
   /*************************************************************************/
public:
   int DAGH_AddGridFunction(GridFunctionVoid *gfv);
   inline void DAGH_DelGridFunction(const int gfid)
        { 
         gflist[gfid] = (GridFunctionVoid *) 0; 
         if (comm_service::dce()) comm_service::delete_comm(gfid);
         gfnum--; 
        }

   inline void DAGH_ComposeGridFunctions(void)
        { 
         for (register int i=0; i<gfnum; i++) 
          { if (gflist[i]) gflist[i]->GF_Compose(); }
         if (chkpt_restart()) {
          for (register int i=0; i<gfnum; i++) 
           { if (gflist[i] && gflist[i]->checkpoint()) 
             gflist[i]->GF_CheckpointRestart(); }
         }
        }

   inline int DAGH_GridFunctionTemplate(GridFunctionVoid& gfv)
        { return (DAGHTemplateComm + gfv.gfid); }

   /*************************************************************************/
   /* Recomposition */
   /*************************************************************************/
public:
   void DAGH_RecomposeHierarchy(int Partition); 
   inline void DAGH_RecomposeHierarchy() { DAGH_RecomposeHierarchy(DAGHTrue); }

public:
   void DAGH_Refine(const int lev);
   void DAGH_Refine(BBoxList &bblist, const int lev);

   /*************************************************************************/
   /* Redistribution methods */
   /*************************************************************************/
public:
   inline void DAGH_RedistributeHierarchy(const int type)
        { if (distribution) distribution->set_dist(type); DAGH_RecomposeHierarchy(); }
   inline void DAGH_RedistributeHierarchy(const int type, BBoxList& bbl)
        { if (distribution) distribution->set_dist(type,bbl); DAGH_RecomposeHierarchy(); }

   /*************************************************************************/
   /* For Checkpointing.... */
   /*************************************************************************/
public:
   void DAGH_Checkpoint(const char* name);
   void DAGH_Checkpoint(std::stringstream& ofs);
   int  DAGH_Checkpoint_StrStream_Memory();
   
   void DAGH_ComposeHierarchy(const char* name);
   bool DAGH_RecomposeHierarchy(const char* name);
   void DAGH_RecomposeHierarchy(std::stringstream& ifs);

   void DAGH_InitChkpt(const char* name) { 
     chkptflag = DAGHTrue; 
     std::strncpy(chkptname,name,DAGHChkPtTagNameSize); 
     chkptname[DAGHChkPtTagNameSize-1] = '\0'; 
   }
   void DAGH_InitChkpt()
     { chkptflag = DAGHTrue; }
   void DAGH_StopChkpt()
     { chkptflag = DAGHFalse; chkptname[0] = '\0'; chkptpnum = 0; }

   int DAGH_OpenChkptIStream(const int p,
                             std::ifstream& ifs);
   void DAGH_CloseChkptIStream(std::ifstream& ifs);

   int DAGH_OpenChkptIStream(const int p,
                             std::streampos*& chkptdir, 
                             char*& chkptnamedir,
                             std::ifstream& ifs);
   int DAGH_OpenChkptStrStream(std::streampos*& chkptdir, 
                               char*& chkptnamedir,
                               std::stringstream& ifs);

   void DAGH_CloseChkptIStream(std::streampos*& chkptdir, 
                               char*& chkptnamedir,
                               std::ifstream& ifs);

   void DAGH_CloseChkptStrStream(std::streampos*& chkptdir, 
                                 char*& chkptnamedir,
                                 std::stringstream& ifs)
     { if (chkptdir) delete [] chkptdir; 
       if (chkptnamedir) delete [] chkptnamedir; }

   int DAGH_GetGFChkptIStream(const int p, const char* gfname, 
                              const int gfid,
                              std::ifstream& ifs);
   int DAGH_GetGFChkptStrStream(const char* gfname, 
                                const int gfid,
                                std::stringstream& ifs);

   int DAGH_OpenChkptOStream(const int p,
                             std::ofstream& ofs);
   void DAGH_CloseChkptOStream(std::ofstream& ofs);

   /*************************************************************************/
   /* Communications */
   /*************************************************************************/
private:
   void DAGH_GlbConcat(void *snddata, int sndsize, void *&rcvdata, 
                       int &rcvsize, MPI_Comm Cin=0);
   //void DAGH_GlbConcat(void *snddata, const int size, 
   //                    void *&rcvdata, MPI_Comm Cin=0);
   //void DAGH_PingIONode(const int tag, int flg = DAGHTrue);
   void DAGH_PingIONode(const int tag, int flg);
public:
   int DAGH_CommInit(MPI_Comm c = 0);
   void DAGH_CommKill(void);

    /*************************************************************************/
    /* Initialize the overlap flag array */
    /*************************************************************************/ 
public:
   static void DAGH_InitOverlaps(const int Type, short* olap)
   { if (Type == DAGHCellCentered) 
       {olap[0] = 0; olap[1] = 0; olap[2] = 0; }
     else if (Type == DAGHVertexCentered) 
       {olap[0] = 1; olap[1] = 1; olap[2] = 1; } 
     else if (Type == DAGHFaceCentered_X) 
       {olap[0] = 1; olap[1] = 0; olap[2] = 0; }  
     else if (Type == DAGHFaceCentered_Y) 
       {olap[0] = 0; olap[1] = 1; olap[2] = 0; }  
     else if (Type == DAGHFaceCentered_Z) 
       {olap[0] = 0; olap[1] = 0; olap[2] = 1; } }
   /*************************************************************************/

   /*************************************************************************/
   /* Private query operators */
   /*************************************************************************/
private:   
   inline int daghoverlapCC_or_NCC() const 
     { return ((daghtype==DAGHCellCentered) ? 0 : 1); }

   inline const short* daghoverlap() const 
     { return (overlap); }

   inline int daghoverlap(const int dir) const 
     { return (overlap[dir]); }
   /*************************************************************************/
   /* Static Generic public query operators */
   /*************************************************************************/
public:
   static int daghoverlapCC_or_NCC(const int type) 
      { return((type==DAGHCellCentered) ? 0 : 1); }
   
   /*************************************************************************/
   /* Generic public query operators */
   /*************************************************************************/
public:
   inline int default_alignment(const int r) const 
        { return ((r == rank) ? DAGH_All : 
                  (r == 3) ? DAGH_All : 
                  (r == 2) ? DAGH_XY : DAGH_X); }
   /*************************************************************************/
   /* GridHierarchy public query operators */
   /*************************************************************************/
public:
   inline int dagh_type(void) const { return daghtype; }

   inline int boundarywidth(void) const { return bndrywidth; }
   inline int boundarytype(void) const { return bndrytype; }
   inline int adaptboundarytype(void) const { return adaptbndrytype; }
   inline int periodicboundary(const int d) const { return bndryperiodic[d]; }
   inline int periodicindex(const int d) const { return periodicidx[d]; }

   inline int externalghostwidth(void) const { return extghostwidth; }

   inline int totallevels(void) const { return levels; }

   inline int coarselevel(void) const { return crslev; }
   inline int maxlevel(void) const { return maxlev; }
   inline int finelevel(void) const { return finelev; }

   inline int stepsize(const int lev) const { 
     int ret=1; for (register int l=Max(0,lev); l<levels-1; l++) ret *= refby[l]; return (ret); 
   }
   inline int timestep(const int lev) const { 
     int ret=1; for (register int l=Max(0,lev); l<levels-1; l++) ret *= refby[l]; return (ret); 
   }
   inline int stepstaken(const int t, const int l) const
     { return (t/timestep(l)); }
   inline double delta_t(const int l) const
     { return (tspecs.h/refinedby(l)); }
   inline double delta_x(const int dim, const int l) const
     { return (dspecs[dim].h/refinedby(l)); }

   inline int refinefactor(const int lev) const 
     { return (lev>=0 && lev<levels ? refby[lev] : 1); }
   inline int refinedby(const int lev) const {
     int ret=1; for (register int l=0; l<Min(lev,levels); l++) ret *= refby[l]; return (ret);
   }
   int refinedby(const int l1, const int l2) const { 
     int ret=1;
     if (l1>=0 && l2>=0 && l1<levels && l2<levels)
       if (l2<l1) 
         for (register int l=l2; l<l1; l++) 
           ret *= refby[l]; 
       else if (l2>l1) {
         ret = -1;
         for (register int l=l1; l<l2; l++) 
           ret *= refby[l]; 
       }
     return (ret);
   }

   inline int buffer() const { return nestbuffer; }

   inline const int& gucfactor() const { return cfactor; }
   inline const int& guminwidth() const { return guwidth; }   
   inline void setgucfactor(const int& cfac) { 
     cfactor = cfac; 
     guwidth = (int)(std::log((double)cfactor)/std::log((double)2.0));
     if (guwidth<std::log((double)cfactor)/std::log((double)2.0)) guwidth++;
   }

   inline const int& maxGridBoxSize() const { return maxgridsize; }
   inline void setmaxGridBoxSize(const int& mgsize) { maxgridsize = mgsize; }

   /*************************************************************************/
   /* Set/Query Updated Timestep Info */
   /*************************************************************************/
   inline int updatedvaluestep() const { return updatedstep; } 
   inline void setupdatedvaluestep(const int ustep) { updatedstep = ustep; } 

   /*************************************************************************/
   /* Checkpoint info... */
   /*************************************************************************/
   void chkpt_filename(std::ostringstream& obuf, 
                       const char* name, const int p);
   inline int chkpt_restart(void) const 
     { return (chkptflag ==  DAGHTrue); }
   inline const char* chkpt_name(void) const 
     { return chkptname; }
   inline int chkpt_pnum(void) const
     { return chkptpnum; }
   inline void chkpt_reset(void)
     { chkptflag = DAGHFalse; chkptname[0] = '\0'; chkptpnum = 0; }

   /*************************************************************************/
   /* Set/Get current time */
   /*************************************************************************/
public:
   inline int getCurrentTime(const int lev) const 
     { return curtime[lev] ; }
   inline void setCurrentTime(const int ctime, const int lev) 
     { curtime[lev] = ctime; }
   inline void incrCurrentTime(const int lev) 
     { curtime[lev] += timestep(lev); }

   inline int getPreviousTime(const int lev) const 
     { return curtime[lev] - timestep(lev); }
   inline int getNextTime(const int lev) const 
     { return curtime[lev] + timestep(lev); }

   /*************************************************************************/
   /* DAGH structure & coordinate information */
   /*************************************************************************/
public:
   inline int GlobalLength(int lev) { return (globalgbl[lev] ? globalgbl[lev]->number() : 0); }
   inline int GlobalMaxIndex(int lev) { return (globalgbl[lev] ? globalgbl[lev]->maxindex() : 0); }
   inline int LocalLength(int lev)  { return (localgbl[lev] ? localgbl[lev]->number() : 0); }

   inline GridUnitList* clist() { return (complist); }
   inline GridUnitList* llist() { return (locallist); }

   inline GridBoxList** lgbl() { return (localgbl); }
   inline GridBoxList** ggbl() { return (globalgbl); }
   inline GridBoxList** oldlgbl() { return (oldlocalgbl); }
   inline GridBoxList** oldggbl() { return (oldglobalgbl); }

   inline GridBoxList* lgbl(int lev) { return (localgbl[lev]); }
   inline GridBoxList* ggbl(int lev) { return (globalgbl[lev]); }
   inline GridBoxList* oldlgbl(int lev) { return (oldlocalgbl[lev]); }
   inline GridBoxList* oldggbl(int lev) { return (oldglobalgbl[lev]); }

   /* Integer coords */
   void llb(int *lc) const;
   void lub(int *lc) const;

   inline Coords llb(void) const { return (basebbox.lower()); }
   inline Coords lub(void) const { return (basebbox.upper()); }

   /* World coords */
   void wlb(double *wc) const;
   void wub(double *wc) const;
  
   DCoords worldCoords(const int *lc, const int *ls) const;
   DCoords worldCoordsUpper(const int *lc, const int *ls) const;
   DCoords worldCoordsLower(const int *lc, const int *ls) const;
   DCoords worldStep(const int *ls) const;

   void worldCoords(const int *lc, const int *ls, double *wc) const;
   void worldStep(const int *ls, double *ws) const;

   /* To integer coords */
   Coords localCoords(double const *wc) const;
   Coords localStep(double const *ws) const;

   void localCoords(double const *wc, int *lc) const;
   void localStep(double const *ws, int *ls) const;

   /*************************************************************************/
   /* Bounding box information */
   /*************************************************************************/
public:
   inline BBox bbbox(void) const { return basebbox; }
   inline BBox glbbbox(void) const { return intbbox; }
   inline BBox wholebbox(void) const 
     { return shrinkupperbydim(intbbox,daghoverlap()); }
   inline BBoxList* wholebaselist(void) const { return intbaselist; };
   inline const double* wholebndry(void) const { return bndrycoords; }
   inline const int& nbndry(void) const { return nbndrycoords; }

   /* Global BBox Lists */
   void glb_bboxlist(BBoxList& bbl,
                     const int l, 
                     const int type=DAGHCellCentered) const;
   /*************************************************************************/
   /* IO Support */
   /*************************************************************************/
public:
   void DAGH_IOServe(void);
   inline void DAGH_SetIOType(const short type) { io_type = type; }
   inline const short& DAGH_IOType(void) { return io_type; }
   inline DAGHIO_WriteFunc DAGH_IOWrite(void) { return DAGHIO_Write[io_type]; }
   inline DAGHIO_ReadFunc DAGH_IORead(void) { return DAGHIO_Read[io_type]; }
   void DAGH_IOFlush(void);
   void DAGH_IOClose(void);
   void DAGH_IOEnd(void);
   friend void DAGHIO_EndIOPingFunction(GridHierarchy &GH);

   /*************************************************************************/
   /* Ghost Communications */
   /*************************************************************************/
public:
   inline void DAGH_SyncGhostRegions(const int t, const int l)
     { 
      register int i;
      for (i=0; i<gfnum; i++)
        if (gflist[i] && gflist[i]->comm()) 
          gflist[i]->GF_WriteGhosts(t,l);
      for (i=0; i<gfnum; i++)
        if (gflist[i] && gflist[i]->comm()) 
          gflist[i]->GF_ReadGhosts(t,l);
     }

   inline void DAGH_SyncGhostRegions(const int t, const int l, 
                                     const int axis, const int dir)
     { 
      register int i;
      for (i=0; i<gfnum; i++) 
        if (gflist[i] && gflist[i]->comm()) 
          gflist[i]->GF_WriteGhosts(t,l,axis,dir);
      for (i=0; i<gfnum; i++) 
        if (gflist[i] && gflist[i]->comm()) 
          gflist[i]->GF_ReadGhosts(t,l,axis,dir);
     }
  };

std::ostream& operator<<(std::ostream&, const GridHierarchy&);
std::ostream& operator<<(std::ostream&, const struct dimspecs&);
std::ostream& operator<<(std::ostream&, const struct MergedGridUnit&);

std::ofstream& operator<<(std::ofstream&, const GridHierarchy&);
std::ifstream& operator>>(std::ifstream&, GridHierarchy&);
std::stringstream& operator<<(std::stringstream&, const GridHierarchy&);
std::stringstream& operator>>(std::stringstream&, GridHierarchy&);


#endif

---