Merge branch 'charm' of charmgit:charm into charm
[charm.git] / examples / amr / cfd / cfdAMR.C
1 #include "cfdAMR.h"
2 #include "cfdAMR.decl.h"
3
4 /*
5 ********************************************
6 User Code for 2D
7 ********************************************
8 */
9 void AmrUserDataCfdProblemInit(void)
10 {
11   PUPable_reg(AmrUserData);
12   PUPable_reg(CfdProblem);
13 }
14
15 void CfdProblem::pup(PUP::er &p)
16 {
17   AmrUserData::pup(p);
18
19   p|gridH;
20   p|gridW;
21   p|stepCount;
22   if(p.isUnpacking()) {
23      grid = new CfdGrid();
24   }
25   grid->pup(p);
26   
27 }
28
29
30 AmrUserData* AmrUserData :: createData()
31 {
32   CfdProblem *instance = new CfdProblem;
33   return (AmrUserData *)instance;
34 }
35   
36 AmrUserData* AmrUserData :: createData(void *data, int dataSize)
37 {
38    CfdProblem *instance = new CfdProblem(data, dataSize);
39    return (AmrUserData *) instance;
40 }
41
42
43 void AmrUserData :: deleteNborData(void* data)
44 {
45   delete [](CfdLoc *) data;
46 }
47
48 void AmrUserData :: deleteChildData(void* data)
49 {
50   delete [](CfdLoc *) data;
51 }
52
53 void CfdProblem :: doComputation(void)
54 {
55   startStep();
56   stepCount++;
57   finishStep(false);
58 }
59
60 void ** CfdProblem ::getNborMsgArray(int* sizePtr)
61 {
62   CfdLoc ** dataArray = new CfdLoc*[4];
63   int x,y;
64   //allocate space for the +x and -x nbor data
65   for(int i=0; i<2; i++) {
66     dataArray[i] = new CfdLoc[gridH];
67     sizePtr[i] = gridH * sizeof(CfdLoc);
68   }
69   //allocate space for the +y and -y nbor data
70   for(int i=2; i<4; i++) {
71     dataArray[i] = new CfdLoc[gridW];
72     sizePtr[i] = gridW * sizeof(CfdLoc);
73   }
74
75   CfdLoc** mygrid = grid->getPrimary();
76   //For +ve X nbor & For -ve X nbor
77   //  x = gridW; x = 1;
78   for(int y=1; y< gridH+GHOSTWIDTH-1; y++) {
79     dataArray[0][y-1].copy(&mygrid[gridW][y]);
80     dataArray[1][y-1].copy(&mygrid[1][y]);
81   }
82   //For +ve Y Nbor & For -ve Y Nbor
83   //  y = gridH;  y = 1;
84   for(int x =1; x<gridW+GHOSTWIDTH-1;x++) {
85     dataArray[2][x-1].copy(&mygrid[x][gridH]);
86     dataArray[3][x-1].copy(&mygrid[x][1]);
87   }
88   
89   return (void **) dataArray;
90 }
91
92
93 void CfdProblem :: store(void* data, int dataSize, int neighborSide)
94 {
95   CfdLoc* dataArray = (CfdLoc*) data;
96   int x,y;
97   switch(neighborSide) {
98   case 0:
99     if(dataSize/sizeof(CfdLoc) == gridH) {
100       CfdLoc** mygrid = grid->getPrimary();
101       x=0;
102       for(int y=1;y<gridH+GHOSTWIDTH-1;y++)
103         mygrid[x][y].copy(&(dataArray[y-1]));
104     }
105     break;
106   case 1:
107     if(dataSize/sizeof(CfdLoc) == gridH) {
108       CfdLoc** mygrid = grid->getPrimary();
109       x=gridW+1;
110       for(int y=1;y<gridH+GHOSTWIDTH-1;y++)
111         mygrid[x][y].copy(&(dataArray[y-1]));
112     }
113     break;
114   case 2:
115     if(dataSize/sizeof(CfdLoc) == gridW) {
116       CfdLoc** mygrid = grid->getPrimary();
117       y=0;
118       for(int x=1;x<gridW+GHOSTWIDTH-1;x++)
119         mygrid[x][y].copy(&(dataArray[x-1]));
120     }
121     break;
122   case 3:
123     if(dataSize/sizeof(CfdLoc) == gridW) {
124       CfdLoc** mygrid = grid->getPrimary();
125       y=gridH+1;
126       for(int x=1;x<gridW+GHOSTWIDTH-1;x++)
127         mygrid[x][y].copy(&(dataArray[x-1]));
128     }
129     break;
130   }
131 }
132
133 bool CfdProblem :: refineCriterion() 
134 {
135   if(grid->averageP() > 10000.0)
136     return true;
137   else
138     return false;
139 }
140
141 void **CfdProblem :: fragmentNborData(void *data, int* sizePtr)
142 {
143   int elements = (*sizePtr)/sizeof(CfdLoc);
144   int newElements = elements/2;
145   CfdLoc **fragmentedArray = new CfdLoc* [2];
146   CfdLoc *indata = (CfdLoc *)data;
147   if(elements %2 == 0){
148     *sizePtr = newElements * sizeof(CfdLoc);
149     for(int i=0; i<2; i++) {
150       fragmentedArray[i] = new CfdLoc[newElements];
151       for(int j=0; j<newElements;j++)
152         fragmentedArray[i][j] = indata[i*newElements + j];
153     }
154   }
155   else {
156     *sizePtr =( ++newElements)*sizeof(CfdLoc);
157     for(int i=0; i<2; i++) {
158       fragmentedArray[i] = new CfdLoc[newElements];
159       for(int j=0; j<newElements-1;j++)
160         fragmentedArray[i][j] = indata[i*newElements + j];
161     }
162     fragmentedArray[1][newElements-1] = indata[elements -1];
163     fragmentedArray[0][newElements-1] = (fragmentedArray[0][newElements -2].average(fragmentedArray[1][0]));
164   }
165   return (void **)fragmentedArray;
166
167 }
168
169 void CfdProblem :: combineAndStore(void **dataArray, int dataSize,int neighborSide) 
170 {
171   int size = dataSize /sizeof(CfdLoc);
172   CfdLoc * buf = new CfdLoc[size];
173   CfdLoc ** data = (CfdLoc**)dataArray;
174   // memcpy((void *)buf, dataArray[0], dataSize);
175   //memcpy((void *)tmpbuf, dataArray[1], dataSize);
176   for(int i=0;i<size;i++)
177     buf[i] = data[0][i].average(data[1][i]);
178   DEBUGJ(("Calling store from combine and store msg size %d\n",dataSize));
179   store((void *)buf,(dataSize),neighborSide);
180   delete []buf;
181 }
182
183 void** CfdProblem :: fragmentForRefine(int *sizePtr)
184 {
185   int newXSize = gridW/2;
186   int newYSize = gridH/2;
187   //  sizePtr = newXSize*newYSize**sizeof(CfdLoc);
188
189   *sizePtr = ((gridW+2)*(gridH+2)+1)*sizeof(CfdLoc);
190   CfdLoc ** dataArray = new CfdLoc* [4];
191   CfdLoc** myGrid = grid->getPrimary();
192
193   for(int i=0;i<4;i++) {
194     /*dataArray[i] = new CfdLoc[newXSize*newYSize+1];
195       for(int x=1;x<=newXSize;x++){
196       for(int y=1;y<=newYSize;y++)
197       dataArray[i][(x-1)*newYSize+(y-1)] = dataGrid[((i/2)%2)*newXSize+x][(i%2)*newYSize+y];
198       }*/
199     dataArray[i] = new CfdLoc[(gridW+2)*(gridH+2)+1];
200
201     for(int y=0; y<gridH+2; y++)
202       for(int x=0; x<gridW+2; x++)
203         dataArray[i][y*(gridW+2)+x].copy(&myGrid[x][y]);
204     //Hack to get the step count to the children
205     dataArray[i][(gridW+2)*(gridH+2)] = CfdLoc((double)stepCount,0,0,0);
206   }
207   return  (void **)dataArray;
208 }
209
210 int CfdProblem :: readParentGrid(void* data, int dataSize)
211 {
212   CfdLoc *dataArray = (CfdLoc*) data;
213   CfdLoc **myGrid = grid->getPrimary();
214   int size = ((gridH+2)*(gridW+2)+1)*sizeof(CfdLoc);
215   if(size == dataSize) {
216     for(int y=0; y<gridH+2; y++)
217       for(int x=0; x<gridW+2; x++) 
218         myGrid[x][y].copy(&dataArray[y*(gridW+2)+x]);
219     //return the step count
220     return (int)(dataArray[(gridW+2)*(gridH+2)].P);
221   }
222   else {
223     CkError("readParentGrid: Error in the size of the parent grid\n");
224     return 0;
225   }
226 }
227
228 Properties* CfdProblem :: initParam() {
229   gridW = gridH =100;
230   double timeStep  = 0.00005;
231   //Physical and simulation constants
232   double meshSize=1.00;//Size of one side of a CFD cell, (m)
233   double diffuse=1.0e-2*meshSize/timeStep;//Diffusion constant, (m/s)
234   double density=1.0e3;//Mass density (Kg/m^3)
235   double modulus=2.0e9;//Bulk modulus (N/m^2)
236   double thik=0.01; //Plate thickness (m)
237   
238   return (new Properties(timeStep,meshSize,diffuse,density,modulus,thik));
239 }
240
241
242
243 void CfdProblem :: initGrid () {
244   //Initialize the interior
245   double varVel=1.0e-6;//Initial velocity variance (m/s)
246   double avgPres=100.0e3;//Initial pressure (N/m^2)
247   double varPres=  0.1e3;//Initial pressure variation (N/m^2)
248   double  backgroundT=300.0;//Normal temperature (K)
249
250   double Vx=varVel*0.05;//rand.nextGaussian();
251   double Vy=varVel*0.05;//rand.nextGaussian();
252   double P=avgPres+varPres*0.05;//rand.nextGaussian();
253   double T=backgroundT;
254
255   for (int y=0;y<gridH;y++)
256     for (int x=0;x<gridW;x++) {
257       CfdLoc* temp = new CfdLoc(P,Vx,Vy,T);
258       //  temp->init(P,Vx,Vy,T);
259       grid->set(x,y,temp);
260       //delete temp;
261     }
262  
263 }
264
265 void CfdProblem :: initBoundary()
266 {
267   double avgPres=100.0e3;//Initial pressure (N/m^2)
268   double  backgroundT=300.0;//Normal temperature (K)
269   //initialize the boundary
270   if(isOnNegXBoundary()) {
271     for (int y=0;y<gridH+GHOSTWIDTH;y++) 
272       grid->set(0,y,new CfdLoc(avgPres,0,0,backgroundT));
273   }
274   
275   if(isOnPosXBoundary()) {
276     for (int y=1;y<gridH+GHOSTWIDTH;y++)
277       grid->set(gridW+1,y,new CfdLoc(avgPres,0,0,backgroundT));
278   }
279   
280   if(isOnNegYBoundary()) {
281     for (int x=1;x<gridW+GHOSTWIDTH;x++) 
282       grid->set(x,0,new CfdLoc(avgPres,0,0,backgroundT));
283   }
284   
285   if(isOnPosYBoundary()) {
286     for (int x=1;x<gridW+GHOSTWIDTH;x++)
287       grid->set(x,gridH+1,new CfdLoc(avgPres,0,0,backgroundT));
288   }
289 }
290
291 void CfdProblem :: startStep()
292 {
293   //stepCount++;
294   //Impose the boundary conditions
295   setBoundaries();
296   //Do pressure calculations
297   grid->update(grid->getPrimary(),grid->getSecondary());
298   grid->setSecondary();
299 }
300
301 //Resample (and switch buffers)
302 void CfdProblem :: finishStep(bool resample)
303 {
304   int presSteps=4;//Number of pressure-only steps before resampling
305   if (resample && stepCount%presSteps==presSteps-1)
306     //Resample new values by new velocities
307     grid->resample(grid->getSecondary(),grid->getPrimary(),
308                    presSteps);
309   else
310     grid->swap();
311   grid->setPrimary();
312 }
313
314 void CfdProblem ::setBoundaries()
315 {
316   int x,y;
317   if(isOnNegYBoundary()) {
318     y=0;for (x=0;x<gridW+GHOSTWIDTH;x++) {//Top
319       grid->at(x,y)->Vy=0;
320       grid->at(x,y)->P=grid->at(x,y+1)->P;
321     }
322   }
323   
324   if(isOnPosYBoundary()) {
325     y=gridH-1;for (x=0;x<gridW+GHOSTWIDTH;x++) {//Bottom (out of bounds)
326       grid->at(x,y)->Vy=0;
327       grid->at(x,y)->P=grid->at(x,y-1)->P;
328     }
329   }
330
331   int xL=0,xR=gridW-1;
332
333   if(isOnNegXBoundary()) {
334     for (y=0;y<gridH+GHOSTWIDTH;y++) { 
335       grid->at(xL,y)->Vx=0;
336       grid->at(xL,y)->P=grid->at(xL+1,y)->P;
337     }
338   }
339
340   if(isOnPosYBoundary()) { 
341     for (y=0;y<gridH+GHOSTWIDTH;y++) { 
342       grid->at(xR,y)->Vx=0;
343       grid->at(xR,y)->P=grid->at(xR-1,y)->P;
344     }
345   }
346 }
347
348 void CfdGrid :: update(CfdLoc **src,CfdLoc **dest)
349 {
350   for(int x=1; x<gridW+GHOSTWIDTH-1; x++)
351     for(int y=1; y<gridH+GHOSTWIDTH-1; y++) {
352       dest[x][y].update(properties,&src[x][y-1],
353                         &src[x-1][y],&src[x][y],
354                         &src[x+1][y],&src[x][y+1]);
355     }
356
357 }
358
359 void CfdGrid :: resample(CfdLoc **src,CfdLoc **dest,double velScale)
360 {
361   double scale=0.5*velScale*properties.AGlob;
362
363   for(int destX=1; destX<gridW+GHOSTWIDTH-1; destX++)
364     for(int destY=1; destY<gridH+GHOSTWIDTH-1; destY++) {
365       CfdLoc *srcCur=&(src[destX][destY]);
366       double srcX=destX-scale*(srcCur->Vx+src[destX+1][destY].Vx);
367       double srcY=destY-scale*(srcCur->Vy+src[destX][destY+1].Vy);
368       int ix=(int)srcX;
369       int iy=(int)srcY;
370       if (ix>=0 && iy>=0 && ix<gridW+GHOSTWIDTH-1 && iy<gridH+GHOSTWIDTH-1)
371         dest[destX][destY].interpolate(srcX-ix,srcY-iy,&(src[destX][destY]),
372                                        &(src[ix][iy+1]),&(src[ix+1][iy+1]),
373                                        &(src[ix][iy]),&(src[ix+1][iy]));
374       else
375         dest[destX][destY].copy(srcCur);
376     }
377 }
378
379 void CfdLoc::update(Properties prop,
380                     CfdLoc* t,
381                     CfdLoc* l,CfdLoc* c,CfdLoc* r,
382                     CfdLoc* b) 
383 {
384   //Diffuse out the pressure
385   double neighborP=t->P+b->P+l->P+r->P;
386   double Pd=prop.diffuseGlobL*c->P+prop.diffuseGlobR*neighborP;
387   
388   //Pressure -> velocity
389   Vx=c->Vx+prop.VGlob*(l->P-r->P);
390   Vy=c->Vy+prop.VGlob*(t->P-b->P);
391   
392   //Velocity -> pressure
393   P=Pd+prop.PGlob*(l->Vx-r->Vx  +  t->Vy-b->Vy);                
394 }
395
396 double CfdLoc::interpolate(double dx,double dy,
397                            double tl,double tr,
398                            double bl,double br)
399 {
400   double t=tl+dx*(tr-tl);
401   double b=bl+dx*(br-bl);
402   return b+dy*(t-b);
403 }
404
405 void CfdLoc::interpolate(double dx,double dy,
406                          CfdLoc* src,
407                          CfdLoc* tl,CfdLoc* tr,
408                          CfdLoc* bl,CfdLoc* br) 
409 {
410   Vx=interpolate(dx,dy,tl->Vx,tr->Vx,bl->Vx,br->Vx);
411   Vy=interpolate(dx,dy,tl->Vy,tr->Vy,bl->Vy,br->Vy);
412   P =interpolate(dx,dy,tl->P ,tr->P ,bl->P ,br->P );
413   T =interpolate(dx,dy,tl->T ,tr->T ,bl->T ,br->T );
414 }
415
416 PUPable_def(AmrUserData);
417 PUPable_def(CfdProblem);
418
419
420 #include "cfdAMR.def.h"