Merge branch 'charm' of charmgit:charm into charm
[charm.git] / examples / amr / cfd / cfd.h
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <charm++.h>
5 #define GHOSTWIDTH 2
6
7 class Properties {
8  public:
9   double timeStep; //Delta-t, (s)
10   double meshSize;//Size of one side of a CFD cell, (m)
11   double diffuse;//Diffusion constant, (m/s)
12   double density;//Mass density (Kg/m^3)
13   double modulus;//Bulk modulus (N/m^2)
14   double thickness; //Plate thickness (m)
15   
16   //These are the constant terms in the various state equations
17   double diffuseGlobL,diffuseGlobR; //diffusion equation
18   double VGlob; //velocity update
19   double PGlob; //pressure update
20   double AGlob; //advection equation (see CfdGrid)
21   Properties(){}
22   Properties(double dt,double s,double df,
23              double d,double m,double thik) 
24     {
25       timeStep=dt;
26       meshSize=s;
27       diffuse=df;
28       density=d;
29       modulus=m;
30       
31       double alpha=timeStep/meshSize*diffuse;//Diffusion rate, cells/step
32       diffuseGlobL=1-alpha;
33       diffuseGlobR=0.25*alpha;
34       
35       thickness=thik;
36       double A=meshSize*thickness;//Area of cell-cell interface (m)
37       double V=A*meshSize;//Volume of cell (m^3)
38       VGlob=0.5*timeStep*(1.0/density)*A/V;
39       PGlob=0.5*timeStep*(modulus)*A/V;
40       AGlob=timeStep/meshSize;
41     }
42   void pup(PUP::er& p) {
43     p|timeStep; 
44     p|meshSize;
45     p|diffuse;
46     p|density;
47     p|modulus;
48     p|thickness; 
49     p|diffuseGlobL;
50     p|diffuseGlobR; 
51     p|VGlob; 
52     p|PGlob; 
53     p|AGlob;
54   }
55 };
56
57 class CfdLoc {
58  public:
59   //Everything's measured at the center of the cell
60   double Vx,Vy;//Velocities (m/s)
61   double P;//Pressure (N/m^2)
62   double T;//Temperature (K)
63   CfdLoc(double P0,double Vx0, double Vy0,double T0)
64     {
65       P=P0;
66       Vx=Vx0;
67       Vy=Vy0;
68       T=T0;
69     }
70   
71   CfdLoc() {
72     P=0.0;
73     Vx=0.0;
74     Vy=0.0;
75     T=0.0;
76   }
77
78   /*Perform one physics step, updating our values based on
79     our old values and those of our neighbors.
80   */
81   void update(Properties prop,
82               CfdLoc* t,
83               CfdLoc* l,CfdLoc* c,CfdLoc* r,
84               CfdLoc* b);
85   
86   //Utility: interpolate given values to find (dx,dy)
87   double interpolate(double dx,double dy,double tl,double tr,
88                      double bl,double br);
89   /*Interpolate all quantities from the given corners*/
90   void interpolate(double dx,double dy,
91                    CfdLoc* src,
92                    CfdLoc* tl,CfdLoc* tr,
93                    CfdLoc* bl,CfdLoc* br);
94   
95   /*Copy all quantities from the given*/
96   void copy(CfdLoc* src) {
97     Vx=src->Vx;
98     Vy=src->Vy;
99     P =src->P;
100     T =src->T;
101     //  q =src.q;
102   }
103
104   CfdLoc average(CfdLoc src) {
105     double vx,vy,p,t;
106     vx = (Vx+src.Vx)/2;
107     vy = (Vy+src.Vy)/2;
108     p = (P+src.P)/2;
109     t = (T+src.T)/2;
110     return CfdLoc(p,vx,vy,t);
111   }
112
113   void pup(PUP::er& p) {
114     p|Vx;
115     p|Vy;
116     p|P;
117     p|T;
118   }
119   //Print out debugging info.
120   void printLoc() {
121     printf("P=%d  kpa  Vx=%lf  m/s  Vy=%lf m/s    T=%d K per n \n",
122            (int)(P/1000),Vx,Vy,(int)T);
123   }
124 };
125
126 class CfdGrid
127 {
128  private:
129   //The grid size and storage
130   int gridW,gridH;
131   CfdLoc **cfdNow,**cfdNext, **cur;
132   Properties properties;
133  public:
134   int getWidth() {return gridW;}
135   int getHeight() {return gridH;}
136
137   CfdGrid(){}
138
139   CfdGrid(int w,int h,Properties* prop) {
140     gridW=w;gridH=h;
141     properties=*prop;
142     cfdNow= new CfdLoc* [gridW+GHOSTWIDTH];
143     cfdNext= new CfdLoc* [gridW+GHOSTWIDTH];
144     for(int x=0; x<gridW+GHOSTWIDTH;x++) {
145       cfdNow[x] = new CfdLoc [gridH +GHOSTWIDTH];
146       cfdNext[x] = new CfdLoc [gridH +GHOSTWIDTH];
147     }
148
149     cur=cfdNow;
150     //    delete prop;
151   }
152
153   Properties getProperties() {return properties;}
154   
155   //Get the primary grid
156   CfdLoc** getPrimary() {return cfdNow;}
157   CfdLoc** getSecondary() {return cfdNext;}
158
159   void setPrimary() {cur=cfdNow;}
160   void setSecondary() {cur=cfdNext;}
161   
162   //Swap the primary and secondary grids
163   void swap(){
164     CfdLoc** tmp=cfdNow;
165     cfdNow=cfdNext;
166     cfdNext=tmp; 
167   }
168
169   //Get location x,y from the primary grid
170   CfdLoc* at(int x,int y){
171     return &(cur[x][y]);
172   }
173   
174   //Set for initial conditions
175   void set(int x,int y,CfdLoc* n) {
176     cfdNow[x][y].copy(n);
177     cfdNext[x][y].copy(n);
178     delete n;
179   }
180
181   double averageP(){
182     /*double sumP =0.0;
183     int count;
184     for(int x=1;x<gridW+1;x++)
185       for(int y=1;y<gridH+1;y++)
186       sumP += cur[x][y].P;*/
187     int sumP =0;
188     int count;
189    
190     for(int y=1;y<gridH+1;y++)
191       for(int x=1;x<gridW+1;x++)
192         sumP += (int)cur[x][y].P;
193     count = gridH*gridW;
194     return (double)(sumP/count);
195   }
196
197   //Determine grid values for the next timestep
198   void update(CfdLoc **src,CfdLoc **dest);
199
200   //Resample src into dest based on src velocities
201   void resample(CfdLoc **src,CfdLoc **dest,double velScale);
202
203   void printGrid() {
204     for(int y=0; y<gridH; y++)
205       for(int x=0;x<gridW;x++)
206         cur[x][y].printLoc();
207   }
208   
209   void pup(PUP::er& p) {
210     p|gridW;
211     p|gridH;
212     properties.pup(p);
213     if(p.isUnpacking()) { 
214       cfdNow= new CfdLoc* [gridW+GHOSTWIDTH];
215       cfdNext= new CfdLoc* [gridW+GHOSTWIDTH];
216       for(int x=0; x<gridW+GHOSTWIDTH;x++) {
217         cfdNow[x] = new CfdLoc [gridH +GHOSTWIDTH];
218         cfdNext[x] = new CfdLoc [gridH +GHOSTWIDTH];
219       }
220       cur = cfdNow;
221     }
222    
223     for(int y=0; y<gridH+GHOSTWIDTH; y++)
224       for(int x=0; x<gridW+GHOSTWIDTH;x++){
225         cfdNow[x][y].pup(p);
226         cfdNext[x][y].pup(p);
227       }
228   }
229
230   ~CfdGrid() {
231      for(int x=0;x<gridW+GHOSTWIDTH;x++) {
232       delete []cfdNow[x];
233       delete []cfdNext[x];
234       }
235     delete []cfdNow;
236     delete []cfdNext;
237   }
238         
239 };