Updated Adaptive Mesh Refinement codes which compile against current
authorEric Bohm <ebohm@illinois.edu>
Sun, 22 Apr 2007 19:02:27 +0000 (19:02 +0000)
committerEric Bohm <ebohm@illinois.edu>
Sun, 22 Apr 2007 19:02:27 +0000 (19:02 +0000)
charm API.

13 files changed:
examples/amr/cfd/Makefile [new file with mode: 0644]
examples/amr/cfd/cfd.h [new file with mode: 0644]
examples/amr/cfd/cfdAMR.C [new file with mode: 0644]
examples/amr/cfd/cfdAMR.ci [new file with mode: 0644]
examples/amr/cfd/cfdAMR.h [new file with mode: 0644]
examples/amr/jacobi2D/Makefile [new file with mode: 0644]
examples/amr/jacobi2D/jacobi2DAMR.C [new file with mode: 0644]
examples/amr/jacobi2D/jacobi2DAMR.ci [new file with mode: 0644]
examples/amr/jacobi2D/jacobi2DAMR.h [new file with mode: 0644]
examples/amr/jacobi2Dthesis/Makefile [new file with mode: 0644]
examples/amr/jacobi2Dthesis/jacobi2DAMR.C [new file with mode: 0644]
examples/amr/jacobi2Dthesis/jacobi2DAMR.ci [new file with mode: 0644]
examples/amr/jacobi2Dthesis/jacobi2DAMR.h [new file with mode: 0644]

diff --git a/examples/amr/cfd/Makefile b/examples/amr/cfd/Makefile
new file mode 100644 (file)
index 0000000..c0ffe4b
--- /dev/null
@@ -0,0 +1,132 @@
+
+CHARMDIR=../../..
+CHARMC=$(CHARMDIR)/bin/charmc $(OPTS)
+AMRLIB=$(CHARMDIR)/lib/libamr.a
+
+DEST=cfdAMR
+OBJS =         cfdAMR.o 
+MODS = -module RefineLB
+
+all: $(DEST)
+
+$(DEST): $(OBJS)
+       $(CHARMC) -lm -language charm++ -o $(DEST) $(OBJS) $(AMRLIB) $(MODS)
+       chmod 755 $(DEST)
+
+cfdAMR.o: cfdAMR.C cfdAMR.decl.h cfdAMR.def.h
+       $(CHARMC) -c cfdAMR.C
+
+cfdAMR.decl.h cfdAMR.def.h: cfdAMR.ci
+       $(CHARMC) cfdAMR.ci
+
+clean:
+       rm -f *.o *~ *.de*.h $(DEST) conv-host charmrun
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/examples/amr/cfd/cfd.h b/examples/amr/cfd/cfd.h
new file mode 100644 (file)
index 0000000..97c5d73
--- /dev/null
@@ -0,0 +1,241 @@
+#include <iostream.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <charm++.h>
+#define PI 3.14
+#define GHOSTWIDTH 2
+
+class Properties {
+ public:
+  double timeStep; //Delta-t, (s)
+  double meshSize;//Size of one side of a CFD cell, (m)
+  double diffuse;//Diffusion constant, (m/s)
+  double density;//Mass density (Kg/m^3)
+  double modulus;//Bulk modulus (N/m^2)
+  double thickness; //Plate thickness (m)
+  
+  //These are the constant terms in the various state equations
+  double diffuseGlobL,diffuseGlobR; //diffusion equation
+  double VGlob; //velocity update
+  double PGlob; //pressure update
+  double AGlob; //advection equation (see CfdGrid)
+  Properties(){}
+  Properties(double dt,double s,double df,
+            double d,double m,double thik) 
+    {
+      timeStep=dt;
+      meshSize=s;
+      diffuse=df;
+      density=d;
+      modulus=m;
+      
+      double alpha=timeStep/meshSize*diffuse;//Diffusion rate, cells/step
+      diffuseGlobL=1-alpha;
+      diffuseGlobR=0.25*alpha;
+      
+      thickness=thik;
+      double A=meshSize*thickness;//Area of cell-cell interface (m)
+      double V=A*meshSize;//Volume of cell (m^3)
+      VGlob=0.5*timeStep*(1.0/density)*A/V;
+      PGlob=0.5*timeStep*(modulus)*A/V;
+      AGlob=timeStep/meshSize;
+    }
+  void pup(PUP::er& p) {
+    p|timeStep; 
+    p|meshSize;
+    p|diffuse;
+    p|density;
+    p|modulus;
+    p|thickness; 
+    p|diffuseGlobL;
+    p|diffuseGlobR; 
+    p|VGlob; 
+    p|PGlob; 
+    p|AGlob;
+  }
+};
+
+class CfdLoc {
+ public:
+  //Everything's measured at the center of the cell
+  double Vx,Vy;//Velocities (m/s)
+  double P;//Pressure (N/m^2)
+  double T;//Temperature (K)
+  CfdLoc(double P0,double Vx0, double Vy0,double T0)
+    {
+      P=P0;
+      Vx=Vx0;
+      Vy=Vy0;
+      T=T0;
+    }
+  
+  CfdLoc() {
+    P=0.0;
+    Vx=0.0;
+    Vy=0.0;
+    T=0.0;
+  }
+
+  /*Perform one physics step, updating our values based on
+    our old values and those of our neighbors.
+  */
+  void update(Properties prop,
+             CfdLoc* t,
+             CfdLoc* l,CfdLoc* c,CfdLoc* r,
+             CfdLoc* b);
+  
+  //Utility: interpolate given values to find (dx,dy)
+  double interpolate(double dx,double dy,double tl,double tr,
+                    double bl,double br);
+  /*Interpolate all quantities from the given corners*/
+  void interpolate(double dx,double dy,
+                  CfdLoc* src,
+                  CfdLoc* tl,CfdLoc* tr,
+                  CfdLoc* bl,CfdLoc* br);
+  
+  /*Copy all quantities from the given*/
+  void copy(CfdLoc* src) {
+    Vx=src->Vx;
+    Vy=src->Vy;
+    P =src->P;
+    T =src->T;
+    // q =src.q;
+  }
+
+  CfdLoc average(CfdLoc src) {
+    double vx,vy,p,t;
+    vx = (Vx+src.Vx)/2;
+    vy = (Vy+src.Vy)/2;
+    p = (P+src.P)/2;
+    t = (T+src.T)/2;
+    return CfdLoc(p,vx,vy,t);
+  }
+
+  void pup(PUP::er& p) {
+    p|Vx;
+    p|Vy;
+    p|P;
+    p|T;
+  }
+  //Print out debugging info.
+  void printLoc() {
+    printf("P=%d  kpa  Vx=%lf  m/s  Vy=%lf m/s    T=%d K per n \n",
+          (int)(P/1000),Vx,Vy,(int)T);
+  }
+};
+
+class CfdGrid
+{
+ private:
+  //The grid size and storage
+  int gridW,gridH;
+  CfdLoc **cfdNow,**cfdNext, **cur;
+  Properties properties;
+ public:
+  int getWidth() {return gridW;}
+  int getHeight() {return gridH;}
+
+  CfdGrid(){}
+
+  CfdGrid(int w,int h,Properties* prop) {
+    gridW=w;gridH=h;
+    properties=*prop;
+    cfdNow= new CfdLoc* [gridW+GHOSTWIDTH];
+    cfdNext= new CfdLoc* [gridW+GHOSTWIDTH];
+    for(int x=0; x<gridW+GHOSTWIDTH;x++) {
+      cfdNow[x] = new CfdLoc [gridH +GHOSTWIDTH];
+      cfdNext[x] = new CfdLoc [gridH +GHOSTWIDTH];
+    }
+
+    cur=cfdNow;
+    //    delete prop;
+  }
+
+  Properties getProperties() {return properties;}
+  
+  //Get the primary grid
+  CfdLoc** getPrimary() {return cfdNow;}
+  CfdLoc** getSecondary() {return cfdNext;}
+
+  void setPrimary() {cur=cfdNow;}
+  void setSecondary() {cur=cfdNext;}
+  
+  //Swap the primary and secondary grids
+  void swap(){
+    CfdLoc** tmp=cfdNow;
+    cfdNow=cfdNext;
+    cfdNext=tmp; 
+  }
+
+  //Get location x,y from the primary grid
+  CfdLoc* at(int x,int y){
+    return &(cur[x][y]);
+  }
+  
+  //Set for initial conditions
+  void set(int x,int y,CfdLoc* n) {
+    cfdNow[x][y].copy(n);
+    cfdNext[x][y].copy(n);
+    delete n;
+  }
+
+  double averageP(){
+    /*double sumP =0.0;
+    int count;
+    for(int x=1;x<gridW+1;x++)
+      for(int y=1;y<gridH+1;y++)
+      sumP += cur[x][y].P;*/
+    int sumP =0;
+    int count;
+   
+    for(int y=1;y<gridH+1;y++)
+      for(int x=1;x<gridW+1;x++)
+       sumP += (int)cur[x][y].P;
+    count = gridH*gridW;
+    return (double)(sumP/count);
+  }
+
+  //Determine grid values for the next timestep
+  void update(CfdLoc **src,CfdLoc **dest);
+
+  //Resample src into dest based on src velocities
+  void resample(CfdLoc **src,CfdLoc **dest,double velScale);
+
+  void printGrid() {
+    for(int y=0; y<gridH; y++)
+      for(int x=0;x<gridW;x++)
+       cur[x][y].printLoc();
+  }
+  
+  void pup(PUP::er& p) {
+    p|gridW;
+    p|gridH;
+    properties.pup(p);
+    if(p.isUnpacking()) { 
+      cfdNow= new CfdLoc* [gridW+GHOSTWIDTH];
+      cfdNext= new CfdLoc* [gridW+GHOSTWIDTH];
+      for(int x=0; x<gridW+GHOSTWIDTH;x++) {
+       cfdNow[x] = new CfdLoc [gridH +GHOSTWIDTH];
+       cfdNext[x] = new CfdLoc [gridH +GHOSTWIDTH];
+      }
+      cur = cfdNow;
+    }
+   
+    for(int y=0; y<gridH+GHOSTWIDTH; y++)
+      for(int x=0; x<gridW+GHOSTWIDTH;x++){
+       cfdNow[x][y].pup(p);
+       cfdNext[x][y].pup(p);
+      }
+  }
+
+  ~CfdGrid() {
+     for(int x=0;x<gridW+GHOSTWIDTH;x++) {
+      delete []cfdNow[x];
+      delete []cfdNext[x];
+      }
+    delete []cfdNow;
+    delete []cfdNext;
+  }
+       
+};
diff --git a/examples/amr/cfd/cfdAMR.C b/examples/amr/cfd/cfdAMR.C
new file mode 100644 (file)
index 0000000..b4b5db0
--- /dev/null
@@ -0,0 +1,420 @@
+#include "cfdAMR.h"
+#include "cfdAMR.decl.h"
+
+/*
+********************************************
+User Code for 2D
+********************************************
+*/
+void AmrUserDataCfdProblemInit(void)
+{
+  PUPable_reg(AmrUserData);
+  PUPable_reg(CfdProblem);
+}
+
+void CfdProblem::pup(PUP::er &p)
+{
+  AmrUserData::pup(p);
+
+  p|gridH;
+  p|gridW;
+  p|stepCount;
+  if(p.isUnpacking()) {
+     grid = new CfdGrid();
+  }
+  grid->pup(p);
+  
+}
+
+
+AmrUserData* AmrUserData :: createData()
+{
+  CfdProblem *instance = new CfdProblem;
+  return (AmrUserData *)instance;
+}
+  
+AmrUserData* AmrUserData :: createData(void *data, int dataSize)
+{
+   CfdProblem *instance = new CfdProblem(data, dataSize);
+   return (AmrUserData *) instance;
+}
+
+
+void AmrUserData :: deleteNborData(void* data)
+{
+  delete [](CfdLoc *) data;
+}
+
+void AmrUserData :: deleteChildData(void* data)
+{
+  delete [](CfdLoc *) data;
+}
+
+void CfdProblem :: doComputation(void)
+{
+  startStep();
+  stepCount++;
+  finishStep(false);
+}
+
+void ** CfdProblem ::getNborMsgArray(int* sizePtr)
+{
+  CfdLoc ** dataArray = new CfdLoc*[4];
+  int x,y;
+  //allocate space for the +x and -x nbor data
+  for(int i=0; i<2; i++) {
+    dataArray[i] = new CfdLoc[gridH];
+    sizePtr[i] = gridH * sizeof(CfdLoc);
+  }
+  //allocate space for the +y and -y nbor data
+  for(int i=2; i<4; i++) {
+    dataArray[i] = new CfdLoc[gridW];
+    sizePtr[i] = gridW * sizeof(CfdLoc);
+  }
+
+  CfdLoc** mygrid = grid->getPrimary();
+  //For +ve X nbor & For -ve X nbor
+  //  x = gridW; x = 1;
+  for(int y=1; y< gridH+GHOSTWIDTH-1; y++) {
+    dataArray[0][y-1].copy(&mygrid[gridW][y]);
+    dataArray[1][y-1].copy(&mygrid[1][y]);
+  }
+  //For +ve Y Nbor & For -ve Y Nbor
+  //  y = gridH;  y = 1;
+  for(int x =1; x<gridW+GHOSTWIDTH-1;x++) {
+    dataArray[2][x-1].copy(&mygrid[x][gridH]);
+    dataArray[3][x-1].copy(&mygrid[x][1]);
+  }
+  
+  return (void **) dataArray;
+}
+
+
+void CfdProblem :: store(void* data, int dataSize, int neighborSide)
+{
+  CfdLoc* dataArray = (CfdLoc*) data;
+  int x,y;
+  switch(neighborSide) {
+  case 0:
+    if(dataSize/sizeof(CfdLoc) == gridH) {
+      CfdLoc** mygrid = grid->getPrimary();
+      x=0;
+      for(int y=1;y<gridH+GHOSTWIDTH-1;y++)
+       mygrid[x][y].copy(&(dataArray[y-1]));
+    }
+    break;
+  case 1:
+    if(dataSize/sizeof(CfdLoc) == gridH) {
+      CfdLoc** mygrid = grid->getPrimary();
+      x=gridW+1;
+      for(int y=1;y<gridH+GHOSTWIDTH-1;y++)
+       mygrid[x][y].copy(&(dataArray[y-1]));
+    }
+    break;
+  case 2:
+    if(dataSize/sizeof(CfdLoc) == gridW) {
+      CfdLoc** mygrid = grid->getPrimary();
+      y=0;
+      for(int x=1;x<gridW+GHOSTWIDTH-1;x++)
+       mygrid[x][y].copy(&(dataArray[x-1]));
+    }
+    break;
+  case 3:
+    if(dataSize/sizeof(CfdLoc) == gridW) {
+      CfdLoc** mygrid = grid->getPrimary();
+      y=gridH+1;
+      for(int x=1;x<gridW+GHOSTWIDTH-1;x++)
+       mygrid[x][y].copy(&(dataArray[x-1]));
+    }
+    break;
+  }
+}
+
+bool CfdProblem :: refineCriterion() 
+{
+  if(grid->averageP() > 10000.0)
+    return true;
+  else
+    return false;
+}
+
+void **CfdProblem :: fragmentNborData(void *data, int* sizePtr)
+{
+  int elements = (*sizePtr)/sizeof(CfdLoc);
+  int newElements = elements/2;
+  CfdLoc **fragmentedArray = new CfdLoc* [2];
+  CfdLoc *indata = (CfdLoc *)data;
+  if(elements %2 == 0){
+    *sizePtr = newElements * sizeof(CfdLoc);
+    for(int i=0; i<2; i++) {
+      fragmentedArray[i] = new CfdLoc[newElements];
+      for(int j=0; j<newElements;j++)
+       fragmentedArray[i][j] = indata[i*newElements + j];
+    }
+  }
+  else {
+    *sizePtr =( ++newElements)*sizeof(CfdLoc);
+    for(int i=0; i<2; i++) {
+      fragmentedArray[i] = new CfdLoc[newElements];
+      for(int j=0; j<newElements-1;j++)
+       fragmentedArray[i][j] = indata[i*newElements + j];
+    }
+    fragmentedArray[1][newElements-1] = indata[elements -1];
+    fragmentedArray[0][newElements-1] = (fragmentedArray[0][newElements -2].average(fragmentedArray[1][0]));
+  }
+  return (void **)fragmentedArray;
+
+}
+
+void CfdProblem :: combineAndStore(void **dataArray, int dataSize,int neighborSide) 
+{
+  int size = dataSize /sizeof(CfdLoc);
+  CfdLoc * buf = new CfdLoc[size];
+  CfdLoc ** data = (CfdLoc**)dataArray;
+  // memcpy((void *)buf, dataArray[0], dataSize);
+  //memcpy((void *)tmpbuf, dataArray[1], dataSize);
+  for(int i=0;i<size;i++)
+    buf[i] = data[0][i].average(data[1][i]);
+  DEBUGJ(("Calling store from combine and store msg size %d\n",dataSize));
+  store((void *)buf,(dataSize),neighborSide);
+  delete []buf;
+}
+
+void** CfdProblem :: fragmentForRefine(int *sizePtr)
+{
+  int newXSize = gridW/2;
+  int newYSize = gridH/2;
+  //  sizePtr = newXSize*newYSize**sizeof(CfdLoc);
+
+  *sizePtr = ((gridW+2)*(gridH+2)+1)*sizeof(CfdLoc);
+  CfdLoc ** dataArray = new CfdLoc* [4];
+  CfdLoc** myGrid = grid->getPrimary();
+
+  for(int i=0;i<4;i++) {
+    /*dataArray[i] = new CfdLoc[newXSize*newYSize+1];
+      for(int x=1;x<=newXSize;x++){
+      for(int y=1;y<=newYSize;y++)
+      dataArray[i][(x-1)*newYSize+(y-1)] = dataGrid[((i/2)%2)*newXSize+x][(i%2)*newYSize+y];
+      }*/
+    dataArray[i] = new CfdLoc[(gridW+2)*(gridH+2)+1];
+
+    for(int y=0; y<gridH+2; y++)
+      for(int x=0; x<gridW+2; x++)
+       dataArray[i][y*(gridW+2)+x].copy(&myGrid[x][y]);
+    //Hack to get the step count to the children
+    dataArray[i][(gridW+2)*(gridH+2)] = CfdLoc((double)stepCount,0,0,0);
+  }
+  return  (void **)dataArray;
+}
+
+int CfdProblem :: readParentGrid(void* data, int dataSize)
+{
+  CfdLoc *dataArray = (CfdLoc*) data;
+  CfdLoc **myGrid = grid->getPrimary();
+  int size = ((gridH+2)*(gridW+2)+1)*sizeof(CfdLoc);
+  if(size == dataSize) {
+    for(int y=0; y<gridH+2; y++)
+      for(int x=0; x<gridW+2; x++) 
+       myGrid[x][y].copy(&dataArray[y*(gridW+2)+x]);
+    //return the step count
+    return (int)(dataArray[(gridW+2)*(gridH+2)].P);
+  }
+  else {
+    CkError("readParentGrid: Error in the size of the parent grid\n");
+    return 0;
+  }
+}
+
+Properties* CfdProblem :: initParam() {
+  gridW = gridH =100;
+  double timeStep  = 0.00005;
+  //Physical and simulation constants
+  double meshSize=1.00;//Size of one side of a CFD cell, (m)
+  double diffuse=1.0e-2*meshSize/timeStep;//Diffusion constant, (m/s)
+  double density=1.0e3;//Mass density (Kg/m^3)
+  double modulus=2.0e9;//Bulk modulus (N/m^2)
+  double thik=0.01; //Plate thickness (m)
+  
+  return (new Properties(timeStep,meshSize,diffuse,density,modulus,thik));
+}
+
+
+
+void CfdProblem :: initGrid () {
+  //Initialize the interior
+  double varVel=1.0e-6;//Initial velocity variance (m/s)
+  double avgPres=100.0e3;//Initial pressure (N/m^2)
+  double varPres=  0.1e3;//Initial pressure variation (N/m^2)
+  double  backgroundT=300.0;//Normal temperature (K)
+
+  double Vx=varVel*0.05;//rand.nextGaussian();
+  double Vy=varVel*0.05;//rand.nextGaussian();
+  double P=avgPres+varPres*0.05;//rand.nextGaussian();
+  double T=backgroundT;
+
+  for (int y=0;y<gridH;y++)
+    for (int x=0;x<gridW;x++) {
+      CfdLoc* temp = new CfdLoc(P,Vx,Vy,T);
+      //  temp->init(P,Vx,Vy,T);
+      grid->set(x,y,temp);
+      //delete temp;
+    }
+}
+
+void CfdProblem :: initBoundary()
+{
+  double avgPres=100.0e3;//Initial pressure (N/m^2)
+  double  backgroundT=300.0;//Normal temperature (K)
+  //initialize the boundary
+  if(isOnNegXBoundary()) {
+    for (int y=0;y<gridH+GHOSTWIDTH;y++) 
+      grid->set(0,y,new CfdLoc(avgPres,0,0,backgroundT));
+  }
+  
+  if(isOnPosXBoundary()) {
+    for (int y=1;y<gridH+GHOSTWIDTH;y++)
+      grid->set(gridW+1,y,new CfdLoc(avgPres,0,0,backgroundT));
+  }
+  
+  if(isOnNegYBoundary()) {
+    for (int x=1;x<gridW+GHOSTWIDTH;x++) 
+      grid->set(x,0,new CfdLoc(avgPres,0,0,backgroundT));
+  }
+  
+  if(isOnPosYBoundary()) {
+    for (int x=1;x<gridW+GHOSTWIDTH;x++)
+      grid->set(x,gridH+1,new CfdLoc(avgPres,0,0,backgroundT));
+  }
+}
+
+void CfdProblem :: startStep()
+{
+  //stepCount++;
+  //Impose the boundary conditions
+  setBoundaries();
+  //Do pressure calculations
+  grid->update(grid->getPrimary(),grid->getSecondary());
+  grid->setSecondary();
+}
+
+//Resample (and switch buffers)
+void CfdProblem :: finishStep(bool resample)
+{
+  int presSteps=4;//Number of pressure-only steps before resampling
+  if (resample && stepCount%presSteps==presSteps-1)
+    //Resample new values by new velocities
+    grid->resample(grid->getSecondary(),grid->getPrimary(),
+                  presSteps);
+  else
+    grid->swap();
+  grid->setPrimary();
+}
+
+void CfdProblem ::setBoundaries()
+{
+  int x,y;
+  if(isOnNegYBoundary()) {
+    y=0;for (x=0;x<gridW+GHOSTWIDTH;x++) {//Top
+      grid->at(x,y)->Vy=0;
+      grid->at(x,y)->P=grid->at(x,y+1)->P;
+    }
+  }
+  
+  if(isOnPosYBoundary()) {
+    y=gridH-1;for (x=0;x<gridW+GHOSTWIDTH;x++) {//Bottom (out of bounds)
+      grid->at(x,y)->Vy=0;
+      grid->at(x,y)->P=grid->at(x,y-1)->P;
+    }
+  }
+
+  int xL=0,xR=gridW-1;
+
+  if(isOnNegXBoundary()) {
+    for (y=0;y<gridH+GHOSTWIDTH;y++) { 
+      grid->at(xL,y)->Vx=0;
+      grid->at(xL,y)->P=grid->at(xL+1,y)->P;
+    }
+  }
+
+  if(isOnPosYBoundary()) { 
+    for (y=0;y<gridH+GHOSTWIDTH;y++) { 
+      grid->at(xR,y)->Vx=0;
+      grid->at(xR,y)->P=grid->at(xR-1,y)->P;
+    }
+  }
+}
+
+void CfdGrid :: update(CfdLoc **src,CfdLoc **dest)
+{
+  for(int x=1; x<gridW+GHOSTWIDTH-1; x++)
+    for(int y=1; y<gridH+GHOSTWIDTH-1; y++) {
+      dest[x][y].update(properties,&src[x][y-1],
+                       &src[x-1][y],&src[x][y],
+                       &src[x+1][y],&src[x][y+1]);
+    }
+
+}
+
+void CfdGrid :: resample(CfdLoc **src,CfdLoc **dest,double velScale)
+{
+  double scale=0.5*velScale*properties.AGlob;
+
+  for(int destX=1; destX<gridW+GHOSTWIDTH-1; destX++)
+    for(int destY=1; destY<gridH+GHOSTWIDTH-1; destY++) {
+      CfdLoc *srcCur=&(src[destX][destY]);
+      double srcX=destX-scale*(srcCur->Vx+src[destX+1][destY].Vx);
+      double srcY=destY-scale*(srcCur->Vy+src[destX][destY+1].Vy);
+      int ix=(int)srcX;
+      int iy=(int)srcY;
+      if (ix>=0 && iy>=0 && ix<gridW+GHOSTWIDTH-1 && iy<gridH+GHOSTWIDTH-1)
+       dest[destX][destY].interpolate(srcX-ix,srcY-iy,&(src[destX][destY]),
+                                      &(src[ix][iy+1]),&(src[ix+1][iy+1]),
+                                      &(src[ix][iy]),&(src[ix+1][iy]));
+      else
+       dest[destX][destY].copy(srcCur);
+    }
+}
+
+void CfdLoc::update(Properties prop,
+                   CfdLoc* t,
+                   CfdLoc* l,CfdLoc* c,CfdLoc* r,
+                   CfdLoc* b) 
+{
+  //Diffuse out the pressure
+  double neighborP=t->P+b->P+l->P+r->P;
+  double Pd=prop.diffuseGlobL*c->P+prop.diffuseGlobR*neighborP;
+  
+  //Pressure -> velocity
+  Vx=c->Vx+prop.VGlob*(l->P-r->P);
+  Vy=c->Vy+prop.VGlob*(t->P-b->P);
+  
+  //Velocity -> pressure
+  P=Pd+prop.PGlob*(l->Vx-r->Vx  +  t->Vy-b->Vy);               
+}
+
+double CfdLoc::interpolate(double dx,double dy,
+                          double tl,double tr,
+                          double bl,double br)
+{
+  double t=tl+dx*(tr-tl);
+  double b=bl+dx*(br-bl);
+  return b+dy*(t-b);
+}
+
+void CfdLoc::interpolate(double dx,double dy,
+                        CfdLoc* src,
+                        CfdLoc* tl,CfdLoc* tr,
+                        CfdLoc* bl,CfdLoc* br) 
+{
+  Vx=interpolate(dx,dy,tl->Vx,tr->Vx,bl->Vx,br->Vx);
+  Vy=interpolate(dx,dy,tl->Vy,tr->Vy,bl->Vy,br->Vy);
+  P =interpolate(dx,dy,tl->P ,tr->P ,bl->P ,br->P );
+  T =interpolate(dx,dy,tl->T ,tr->T ,bl->T ,br->T );
+}
+
+PUPable_def(AmrUserData);
+PUPable_def(CfdProblem);
+
+
+#include "cfdAMR.def.h"
diff --git a/examples/amr/cfd/cfdAMR.ci b/examples/amr/cfd/cfdAMR.ci
new file mode 100644 (file)
index 0000000..12bd238
--- /dev/null
@@ -0,0 +1,11 @@
+mainmodule cfdAMR {
+ extern module amr;
+ //extern module CommLB;       
+ //extern module HeapCentLB;
+ extern module RefineLB;
+ initcall void AmrUserDataCfdProblemInit(void);
+ mainchare main
+ {
+   entry main();
+ };
+};
diff --git a/examples/amr/cfd/cfdAMR.h b/examples/amr/cfd/cfdAMR.h
new file mode 100644 (file)
index 0000000..0374a96
--- /dev/null
@@ -0,0 +1,89 @@
+
+#include "amr.h"
+#include "GreedyCommLB.h"
+#include "GreedyLB.h"
+#include "RefineLB.h"
+#include "cfd.h"
+
+#define LB_FREQUENCY 15
+
+class CfdProblem:public AmrUserData {
+ private:
+  int gridW, gridH;
+  CfdGrid *grid;
+   
+  //Split physics
+  int stepCount;
+ public:
+  CfdProblem() {
+    stepCount =0;
+    Properties* properties = initParam();
+    
+    grid=new CfdGrid(gridW,gridH,properties);
+    initGrid();
+  }
+
+  CfdProblem(void* data, int dataSize) {
+    
+    Properties* properties = initParam();
+    grid=new CfdGrid(gridW,gridH,properties);
+    stepCount = readParentGrid(data, dataSize);
+  }
+  
+  void init() {
+    initBoundary();
+    if(stepCount == 0)
+      createInitialPressure();
+  }
+  
+  void createInitialPressure() {
+    if (myIndex.vec[0] == 0 && myIndex.vec[1] == 0 ) {
+      for (int y=0;y<gridH;y++)
+       for (int x=0;x<gridW;x++) 
+         grid->at(x,y)->P+=20.0e3;
+    }
+  }
+  Properties* initParam();
+  void initGrid();
+  void initBoundary();
+  int readParentGrid(void* data, int dataSize);
+  void startStep();
+  void finishStep(bool resample);
+  void setBoundaries();
+  
+  CfdProblem(CkMigrateMessage *m): AmrUserData(m){}
+  PUPable_decl(CfdProblem);
+  virtual void doComputation(void);
+  virtual void **fragmentNborData(void* data, int* sizePtr);
+  virtual void **getNborMsgArray(int *sizeptr);
+  virtual void store(void* data, int dataSize, int neighborSide);
+  virtual void combineAndStore(void **dataArray, int dataSize,int neighborSide);
+  virtual bool refineCriterion(void);
+  virtual void **fragmentForRefine(int *sizePtr);
+  virtual void pup(PUP::er &p);
+  ~CfdProblem() {
+    delete grid;
+  }
+};
+
+
+
+
+class main : public Chare {
+ public:
+  main(CkArgMsg* args) {
+
+    StartUpMsg *msg;
+    msg = new StartUpMsg;
+    msg->synchInterval = 100;
+    msg->depth = 4;
+    msg->dimension = 2;
+    msg->totalIterations = 399;
+    msg->statCollection = 0;
+    //CreateCommLB();
+    //CreateHeapCentLB();
+    CreateRefineLB();
+    CProxy_AmrCoordinator::ckNew(msg,0);
+
+  }
+};
diff --git a/examples/amr/jacobi2D/Makefile b/examples/amr/jacobi2D/Makefile
new file mode 100644 (file)
index 0000000..5485a38
--- /dev/null
@@ -0,0 +1,131 @@
+
+CHARMDIR=../../..
+CHARMC=$(CHARMDIR)/bin/charmc $(OPTS)
+AMRLIB=$(CHARMDIR)/lib/libamr.a
+
+DEST=jacobi
+OBJS =         jacobi2DAMR.o 
+
+all: $(DEST)
+
+$(DEST): $(OBJS)
+       $(CHARMC) -lm -language charm++ -o $(DEST) $(OBJS) $(AMRLIB)
+       chmod 755 $(DEST)
+
+jacobi2DAMR.o: jacobi2DAMR.C jacobi2DAMR.decl.h jacobi2DAMR.def.h
+       $(CHARMC) -g -c jacobi2DAMR.C
+
+jacobi2DAMR.decl.h jacobi2DAMR.def.h: jacobi2DAMR.ci
+       $(CHARMC) jacobi2DAMR.ci
+
+clean:
+       rm -f *.o *~ *.de*.h $(DEST) conv-host charmrun
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/examples/amr/jacobi2D/jacobi2DAMR.C b/examples/amr/jacobi2D/jacobi2DAMR.C
new file mode 100644 (file)
index 0000000..b835fdc
--- /dev/null
@@ -0,0 +1,167 @@
+#include "jacobi2DAMR.h"
+#include "jacobi2DAMR.decl.h"
+/*
+********************************************
+User Code for 2D
+********************************************
+*/
+
+AmrUserData* AmrUserData :: createData()
+{
+  Jacobi2DAMR *instance = new Jacobi2DAMR;
+  return (AmrUserData *)instance;
+}
+
+AmrUserData* AmrUserData :: createData(void *data, int dataSize)
+{
+  Jacobi2DAMR *instance = new Jacobi2DAMR(data , dataSize);
+  return (AmrUserData *) instance;
+}
+
+
+void AmrUserData :: deleteNborData(void* data)
+{
+  delete (double *) data;
+}
+
+void AmrUserData :: deleteChildData(void* data)
+{
+  delete (double *) data;
+}
+
+
+
+void Jacobi2DAMR :: doComputation(void)
+{
+  for(int i=1; i <= cellSize ;i++) 
+    for(int j=1; j<=cellSize;j++) 
+      newDataGrid[i][j] = 0.2 * (dataGrid[i][j-1] + dataGrid[i][j+1] 
+                                +dataGrid[i][j] +dataGrid[i-1][j] +dataGrid[i+1][j]);
+  copyGrid();
+  
+}
+
+void Jacobi2DAMR :: copyGrid(void) 
+{
+  for(int i=0; i< cellSize+2; i++)
+    for(int j=0; j<cellSize+2; j++){
+      dataGrid[i][j] = newDataGrid[i][j];
+      newDataGrid[i][j] = 10.0;
+    }
+}
+
+void ** Jacobi2DAMR :: getNborMsgArray(int* sizePtr)
+{
+  //gives the size of each individual column in bytes
+  *sizePtr = cellSize* sizeof(double);
+  //since we are using 2D mesh so have an array of size 4
+  double ** dataArray = new double* [4];
+  for(int i =0;i<4;i++) {
+    dataArray[i] = new double[cellSize];
+  }
+  //To my Right neighbor
+  copyColumn(dataArray[0], cellSize);
+  //To my left neighbor
+  copyColumn(dataArray[1], 1);
+  //To my Down neighbor
+  copyRow(dataArray[2], cellSize);
+  //To my Up neighbor
+  copyRow(dataArray[3], 1);
+  return (void **) dataArray;
+}
+
+void **Jacobi2DAMR :: fragmentNborData(void *data, int* sizePtr)
+{
+  int elements = (*sizePtr)/sizeof(double);
+  int newElements = elements/2;
+  double **fragmentedArray = new double* [2];
+  double *indata = (double *)data;
+  if(elements %2 == 0){
+    *sizePtr = newElements * sizeof(double);
+    for(int i=0; i<2; i++) {
+      fragmentedArray[i] = new double[newElements];
+      for(int j=0; j<newElements;j++)
+       fragmentedArray[i][j] = indata[i*newElements + j];
+    }
+  }
+  else {
+    *sizePtr =( ++newElements)*sizeof(double);
+    for(int i=0; i<2; i++) {
+      fragmentedArray[i] = new double[newElements];
+      for(int j=0; j<newElements-1;j++)
+       fragmentedArray[i][j] = indata[i*newElements + j];
+    }
+    fragmentedArray[1][newElements-1] = indata[elements -1];
+    fragmentedArray[0][newElements-1] = (fragmentedArray[0][newElements -2] + fragmentedArray[1][0])/2;
+  }
+  return (void **)fragmentedArray;
+    
+}
+
+void Jacobi2DAMR :: store(void* data , int dataSize , int neighborSide)
+{
+  //  if(dataSize>0)
+  //    CkPrintf("neighborSide %d, size %d\n",neighborSide, dataSize);
+  if(dataSize/sizeof(double) == cellSize) {
+    switch(neighborSide) {
+    case 0:
+      copyToColumn((double*) data, 0);
+      break;
+    case 1:
+      copyToColumn((double *) data, cellSize+1);
+      break;
+    case 2:
+      copyToRow((double *) data, 0);
+      break;
+    case 3:
+      copyToRow((double *) data, cellSize+1);
+      break;
+    }
+  }
+  //  else
+    //    CkError("Error: Jacobi::store...wrong sized message size %d cellsize %d\n",
+    //     dataSize/sizeof(double), cellSize);
+}
+
+void Jacobi2DAMR :: combineAndStore(void *data1, void* data2, int dataSize,int neighborSide) {
+  int size = dataSize /sizeof(double);
+  double * buf = new double[2*size];
+  double *tmpbuf = buf + size;
+  memcpy((void *)buf, data1, dataSize);
+  memcpy((void *)tmpbuf, data2, dataSize);
+  DEBUGJ(("Calling store from combine and store msg size %d\n",dataSize));
+  store((void *)buf,(2*dataSize),neighborSide);
+  delete buf;
+}
+
+bool Jacobi2DAMR :: refineCriterion(void) 
+{
+  double average = sumofGrid()/(cellSize*cellSize);
+  //  CkPrintf("refineCriterion: average %lf \n",average);
+  if(average < 15.0 && cellSize >= 4)
+    return true;
+  else
+    return false;
+}
+
+void** Jacobi2DAMR :: fragmentForRefine(int *sizePtr)
+{
+  // CkError("Error: Jacobi ::fragment for refine\n");
+  //return (void **) NULL;
+   int newSize = cellSize/2;
+  *sizePtr = newSize*newSize*sizeof(double);
+
+  double ** dataArray = new double* [4];
+  for(int i=0;i<4;i++) {
+    dataArray[i] = new double[newSize*newSize];
+    for(int j=1;j<=newSize;j++){
+      for(int k=1;k<=newSize;k++)
+       dataArray[i][(j-1)*newSize+(k-1)] = dataGrid[((i/2)%2)*newSize+j][(i%2)*newSize+k];
+    }
+  }
+  return (void **)dataArray;
+  
+}
+PUPable_def(AmrUserData);
+
+#include "jacobi2DAMR.def.h"
diff --git a/examples/amr/jacobi2D/jacobi2DAMR.ci b/examples/amr/jacobi2D/jacobi2DAMR.ci
new file mode 100644 (file)
index 0000000..40231ae
--- /dev/null
@@ -0,0 +1,7 @@
+mainmodule jacobi2DAMR {
+ extern module amr;
+ mainchare main
+ {
+   entry main();
+ };
+};
\ No newline at end of file
diff --git a/examples/amr/jacobi2D/jacobi2DAMR.h b/examples/amr/jacobi2D/jacobi2DAMR.h
new file mode 100644 (file)
index 0000000..0296536
--- /dev/null
@@ -0,0 +1,102 @@
+#include "amr.h"
+
+class Jacobi2DAMR:public AmrUserData {
+ private:
+  int cellSize;
+  double **dataGrid;
+  double **newDataGrid;
+  void copyGrid(void);
+  void copyColumn(double *buf , int colNum) {
+    for(int i=1; i<=cellSize;i++)
+      buf[i-1] = dataGrid[i][colNum];   
+  }
+  
+  void copyRow(double *buf , int rowNum) {
+    for(int i=1; i<=cellSize;i++)
+      buf[i-1] = dataGrid[rowNum][i];   
+  }
+
+  void copyToColumn(double *buf , int colNum) {
+    for(int i=1; i<=cellSize;i++)
+      dataGrid[i][colNum]=  buf[i-1];   
+  }
+  
+  void copyToRow(double *buf , int rowNum) {
+    for(int i=1; i<=cellSize;i++)
+      dataGrid[rowNum][i]= buf[i-1];   
+  }
+  double sumofGrid(void) {
+    double sum = 0.0;
+    for(int i=1;i<cellSize+1;i++)
+      for(int j=1;j<cellSize+1;j++)
+       sum += dataGrid[i][j];
+    return sum;
+  }
+
+ public:
+  Jacobi2DAMR() {
+    cellSize = 64;
+    dataGrid = new double* [cellSize+2];
+    newDataGrid = new double* [cellSize+2];
+    for(int i=0;i< cellSize+2;i++) {
+      dataGrid[i] = new double [cellSize+2];
+      newDataGrid[i] = new double [cellSize +2];
+      for(int k = 0;  k < cellSize+2; k++) {
+       newDataGrid[i][k]=10.0;
+       dataGrid[i][k] = (i+k) *1.0;
+      }
+    }
+  }
+
+  Jacobi2DAMR(int size) {
+    cellSize = size;
+    dataGrid = new double* [cellSize+2];
+    newDataGrid = new double* [cellSize+2];
+    for(int i=0;i< cellSize+2;i++) {
+      dataGrid[i] = new double [cellSize+2];
+      newDataGrid[i] = new double [cellSize +2];
+      for(int k = 0; k < cellSize+2; k++) {
+       newDataGrid[i][k]= 10.0;
+       dataGrid[i][k] = (i+k) *1.0;
+      }
+    }
+  }
+  Jacobi2DAMR(void *data,int dataSize)
+  {
+    double *indata = (double*) data;
+    cellSize = (int) sqrt((double) (dataSize/sizeof(double)));
+    //    cellSize = cellSize/sizeof(double);
+    dataGrid = new double* [cellSize+2];
+    newDataGrid = new double* [cellSize+2];
+    for(int i=0;i< cellSize+2;i++) {
+      dataGrid[i] = new double [cellSize+2];
+      newDataGrid[i] = new double [cellSize +2];
+      for(int k = 0;  k < cellSize+2; k++) {
+       newDataGrid[i][k]= 10.0;
+       if(i== 0 || i == cellSize+1 || k==0 || k== cellSize + 1)
+         dataGrid[i][k] = (i+k) *1.0;
+       else
+         dataGrid[i][k] = indata[(i-1) * cellSize + (k-1)];
+      }
+    }
+  
+  }
+  virtual void doComputation(void);
+  virtual void **fragmentNborData(void* data, int* sizePtr);
+  virtual void **getNborMsgArray(int *sizeptr);
+  virtual void store(void* data, int dataSize, int neighborSide);
+  virtual void combineAndStore(void *data1, void *data2, int dataSize,int neighborSide);
+  virtual bool refineCriterion(void);
+  virtual void **fragmentForRefine(int *sizePtr);
+  ~Jacobi2DAMR() {
+    delete newDataGrid;
+    delete dataGrid;
+  }
+};
+
+class main : public Chare {
+ public:
+  main(CkArgMsg* args) {
+    CProxy_AmrCoordinator::ckNew( new _DMsg);
+  }
+};
diff --git a/examples/amr/jacobi2Dthesis/Makefile b/examples/amr/jacobi2Dthesis/Makefile
new file mode 100644 (file)
index 0000000..95c3d7f
--- /dev/null
@@ -0,0 +1,131 @@
+
+CHARMDIR=../../..
+CHARMC=$(CHARMDIR)/bin/charmc $(OPTS)
+AMRLIB=$(CHARMDIR)/lib/libamr.a
+
+DEST=jacobi
+OBJS =         jacobi2DAMR.o 
+
+all: $(DEST)
+
+$(DEST): $(OBJS)
+       $(CHARMC) -module EveryLB -lm -language charm++ -o $(DEST) $(OBJS) $(AMRLIB) 
+       chmod 755 $(DEST)
+
+jacobi2DAMR.o: jacobi2DAMR.C jacobi2DAMR.decl.h jacobi2DAMR.def.h
+       $(CHARMC) -c jacobi2DAMR.C
+
+jacobi2DAMR.decl.h jacobi2DAMR.def.h: jacobi2DAMR.ci
+       $(CHARMC) jacobi2DAMR.ci
+
+clean:
+       rm -f *.o *~ *.de*.h $(DEST) conv-host charmrun
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/examples/amr/jacobi2Dthesis/jacobi2DAMR.C b/examples/amr/jacobi2Dthesis/jacobi2DAMR.C
new file mode 100644 (file)
index 0000000..a4883f7
--- /dev/null
@@ -0,0 +1,190 @@
+#include "jacobi2DAMR.h"
+#include "jacobi2DAMR.decl.h"
+/*
+********************************************
+User Code for 2D
+********************************************
+*/
+void AmrUserDataJacobiInit(void)
+{
+  PUPable_reg(AmrUserData);
+  PUPable_reg(Jacobi2DAMR);
+}
+
+AmrUserData* AmrUserData :: createData()
+{
+  Jacobi2DAMR *instance = new Jacobi2DAMR;
+  return (AmrUserData *)instance;
+}
+  
+AmrUserData* AmrUserData :: createData(void *data, int dataSize)
+{
+   Jacobi2DAMR *instance = new Jacobi2DAMR(data , dataSize);
+   return (AmrUserData *) instance;
+}
+
+void AmrUserData :: deleteNborData(void* data)
+{
+  delete [](double *) data;
+}
+
+void AmrUserData :: deleteChildData(void* data)
+{
+  delete [](double *) data;
+}
+
+void Jacobi2DAMR :: doComputation(void)
+{
+  for(int i=1; i <= cellSize ;i++) 
+    for(int j=1; j<=cellSize;j++) 
+      newDataGrid[i][j] = 0.2 * (dataGrid[i][j-1] + dataGrid[i][j+1] 
+                                +dataGrid[i][j] +dataGrid[i-1][j] +
+                                 dataGrid[i+1][j]);
+  copyGrid();
+  
+}
+
+void Jacobi2DAMR :: copyGrid(void) 
+{
+  for(int i=0; i< cellSize+2; i++)
+    for(int j=0; j<cellSize+2; j++){
+      dataGrid[i][j] = newDataGrid[i][j];
+      newDataGrid[i][j] = 10.0;
+    }
+}
+
+void ** Jacobi2DAMR :: getNborMsgArray(int* sizePtr)
+{
+  //gives the size of each individual column in bytes
+  *sizePtr = cellSize* sizeof(double);
+  //since we are using 2D mesh so have an array of size 4
+  double ** dataArray = new double* [4];
+  for(int i =0;i<4;i++) {
+    dataArray[i] = new double[cellSize];
+  }
+  //To my Right neighbor
+  copyColumn(dataArray[0], cellSize);
+  //To my left neighbor
+  copyColumn(dataArray[1], 1);
+  //To my Down neighbor
+  copyRow(dataArray[2], cellSize);
+  //To my Up neighbor
+  copyRow(dataArray[3], 1);
+  return (void **) dataArray;
+}
+
+void **Jacobi2DAMR :: fragmentNborData(void *data, int* sizePtr)
+{
+  int elements = (*sizePtr)/sizeof(double);
+  int newElements = elements/2;
+  double **fragmentedArray = new double* [2];
+  double *indata = (double *)data;
+  if(elements %2 == 0){
+    *sizePtr = newElements * sizeof(double);
+    for(int i=0; i<2; i++) {
+      fragmentedArray[i] = new double[newElements];
+      for(int j=0; j<newElements;j++)
+       fragmentedArray[i][j] = indata[i*newElements + j];
+    }
+  }
+  else {
+    *sizePtr =( ++newElements)*sizeof(double);
+    for(int i=0; i<2; i++) {
+      fragmentedArray[i] = new double[newElements];
+      for(int j=0; j<newElements-1;j++)
+       fragmentedArray[i][j] = indata[i*newElements + j];
+    }
+    fragmentedArray[1][newElements-1] = indata[elements -1];
+    fragmentedArray[0][newElements-1] = (fragmentedArray[0][newElements -2] 
+                                       + fragmentedArray[1][0])/2;
+  }
+  return (void **)fragmentedArray;
+    
+}
+
+void Jacobi2DAMR :: store(void* data , int dataSize , int neighborSide)
+{
+  if(dataSize/sizeof(double) == cellSize) {
+    switch(neighborSide) {
+    case 0:
+      copyToColumn((double*) data, 0);
+      break;
+    case 1:
+      copyToColumn((double *) data, cellSize+1);
+      break;
+    case 2:
+      copyToRow((double *) data, 0);
+      break;
+    case 3:
+      copyToRow((double *) data, cellSize+1);
+      break;
+    }
+  }
+  //  else
+  //    CkError("Error: Jacobi::store...wrong sized message size %d size/sizeof(double) %d cellsize %d\n",
+  //       dataSize,dataSize/sizeof(double), cellSize);
+}
+
+void Jacobi2DAMR :: combineAndStore(void **dataArray, 
+                                       int dataSize,
+                                       int neighborSide) {
+  int size = dataSize /sizeof(double);
+  double * buf = new double[2*size];
+  double *tmpbuf = buf + size;
+  memcpy((void *)buf, dataArray[0], dataSize);
+  memcpy((void *)tmpbuf, dataArray[1], dataSize);
+  store((void *)buf,(2*dataSize),neighborSide);
+  delete []buf;
+}
+
+bool Jacobi2DAMR :: refineCriterion(void) 
+{
+  double average = sumofGrid()/(cellSize*cellSize);
+  if(average < 15.0 && cellSize >= 4)
+    return true;
+  else
+    return false;
+}
+
+void** Jacobi2DAMR :: fragmentForRefine(int *sizePtr)
+{
+  int newSize = cellSize/2;
+  *sizePtr = newSize*newSize*sizeof(double);
+
+  double ** dataArray = new double* [4];
+  for(int i=0;i<4;i++) {
+    dataArray[i] = new double[newSize*newSize];
+    for(int j=1;j<=newSize;j++){
+      for(int k=1;k<=newSize;k++)
+       dataArray[i][(j-1)*newSize+(k-1)] = 
+               dataGrid[((i/2)%2)*newSize+j][(i%2)*newSize+k];
+    }
+  }
+  return (void **)dataArray;
+  
+}
+
+void Jacobi2DAMR::pup(PUP::er &p)
+{
+  AmrUserData::pup(p);
+  p(cellSize);
+
+  //have to pup dataGrid and newDataGrid here
+  if(p.isUnpacking()) {
+    dataGrid = new double*[cellSize+2];
+    newDataGrid = new double*[cellSize+2];
+    for(int i=0;i<cellSize+2;i++) {
+      dataGrid[i] = new double[cellSize+2];
+      newDataGrid[i] = new double[cellSize+2];
+    }
+  }
+  for(int i=0; i<cellSize+2;i++) {
+    p(dataGrid[i],cellSize+2);
+    p(newDataGrid[i],cellSize+2);
+  }
+}
+
+PUPable_def(AmrUserData);
+PUPable_def(Jacobi2DAMR);
+
+#include "jacobi2DAMR.def.h"
diff --git a/examples/amr/jacobi2Dthesis/jacobi2DAMR.ci b/examples/amr/jacobi2Dthesis/jacobi2DAMR.ci
new file mode 100644 (file)
index 0000000..c96b2ac
--- /dev/null
@@ -0,0 +1,9 @@
+mainmodule jacobi2DAMR {
+ extern module amr;
+// extern module GreedyLB;
+ initcall void AmrUserDataJacobiInit(void);
+ mainchare main
+ {
+   entry main();
+ };
+};
diff --git a/examples/amr/jacobi2Dthesis/jacobi2DAMR.h b/examples/amr/jacobi2Dthesis/jacobi2DAMR.h
new file mode 100644 (file)
index 0000000..9befa00
--- /dev/null
@@ -0,0 +1,136 @@
+
+#include "amr.h"
+//#include "GreedyLB.h"
+#define LB_FREQUENCY 15
+
+class Jacobi2DAMR:public AmrUserData {
+ private:
+  int cellSize;
+  double **dataGrid;
+  double **newDataGrid;
+  /*Utility Functions*/
+  void copyGrid(void);
+  void copyColumn(double *buf , int colNum) {
+    for(int i=1; i<=cellSize;i++)
+      buf[i-1] = dataGrid[i][colNum];   
+  }
+  
+  void copyRow(double *buf , int rowNum) {
+    for(int i=1; i<=cellSize;i++)
+      buf[i-1] = dataGrid[rowNum][i];   
+  }
+
+  void copyToColumn(double *buf , int colNum) {
+    for(int i=1; i<=cellSize;i++)
+      dataGrid[i][colNum]=  buf[i-1];   
+  }
+  
+  void copyToRow(double *buf , int rowNum) {
+    for(int i=1; i<=cellSize;i++)
+      dataGrid[rowNum][i]= buf[i-1];   
+  }
+  double sumofGrid(void) {
+    double sum = 0.0;
+    for(int i=1;i<cellSize+1;i++)
+      for(int j=1;j<cellSize+1;j++)
+       sum += dataGrid[i][j];
+    return sum;
+  }
+
+ public:
+  /*Default Constructor: Called in the initial setup of the tree*/
+  Jacobi2DAMR() {
+    cellSize = 32;
+    dataGrid = new double* [cellSize+2];
+    newDataGrid = new double* [cellSize+2];
+    for(int i=0;i< cellSize+2;i++) {
+      dataGrid[i] = new double [cellSize+2];
+      newDataGrid[i] = new double [cellSize +2];
+      for(int k = 0;  k < cellSize+2; k++) {
+       newDataGrid[i][k]=10.0;
+       dataGrid[i][k] = (i+k) *1.0;
+      }
+    }
+  }
+
+  Jacobi2DAMR(int size) {
+    cellSize = size;
+    dataGrid = new double* [cellSize+2];
+    newDataGrid = new double* [cellSize+2];
+    for(int i=0;i< cellSize+2;i++) {
+      dataGrid[i] = new double [cellSize+2];
+      newDataGrid[i] = new double [cellSize +2];
+      for(int k = 0; k < cellSize+2; k++) {
+       newDataGrid[i][k]= 10.0;
+       dataGrid[i][k] = (i+k) *1.0;
+
+
+
+      }
+    }
+  }
+
+  /*This constructor is called after refinement with data from te parent*/
+  Jacobi2DAMR(void *data,int dataSize)
+  {
+    double *indata = (double*) data;
+    cellSize = (int) sqrt((double) (dataSize/sizeof(double)));
+    //    cellSize = cellSize/sizeof(double);
+    dataGrid = new double* [cellSize+2];
+    newDataGrid = new double* [cellSize+2];
+    for(int i=0;i< cellSize+2;i++) {
+      dataGrid[i] = new double [cellSize+2];
+      newDataGrid[i] = new double [cellSize +2];
+      for(int k = 0;  k < cellSize+2; k++) {
+       newDataGrid[i][k]= 10.0;
+       if(i== 0 || i == cellSize+1 || k==0 || k== cellSize + 1)
+         dataGrid[i][k] = (i+k) *1.0;
+       else
+         dataGrid[i][k] = indata[(i-1) * cellSize + (k-1)];
+      }
+    }
+  
+  }
+
+  Jacobi2DAMR(CkMigrateMessage *m): AmrUserData(m){}
+
+  PUPable_decl(Jacobi2DAMR);
+  /*Mandatory Library Interface functions*/
+  virtual void doComputation(void);
+  virtual void **fragmentNborData(void* data, int* sizePtr);
+  virtual void **getNborMsgArray(int *sizeptr);
+  virtual void store(void* data, int dataSize, int neighborSide);
+  virtual void combineAndStore(void **dataArray, int dataSize,int neighborSide);
+  virtual bool refineCriterion(void);
+  virtual void **fragmentForRefine(int *sizePtr);
+
+  /*If load balancing is required*/
+  virtual void pup(PUP::er &p);
+  /*Destructor*/
+  ~Jacobi2DAMR() {
+    for (int i=0; i< cellSize+2;i++)
+      delete [] newDataGrid[i];
+    delete [] newDataGrid;
+    for (int i=0; i< cellSize+2;i++)
+      delete [] dataGrid[i];
+    delete[] dataGrid;
+  }
+};
+
+/*Main Chare*/
+class main : public Chare {
+ public:
+  /*Constructor: Library is created from here*/
+  main(CkArgMsg* args) {
+
+    StartUpMsg *msg;
+    msg = new StartUpMsg;
+    msg->synchInterval = 200;
+    msg->depth = 2;
+    msg->dimension = 2;
+    msg-> totalIterations = 500;
+    //    CreateGreedyLB();
+    CProxy_AmrCoordinator::ckNew(msg,0);
+
+  }
+};