examples/charm++: Always use CkWallTimer, not CmiWallTimer (tested)
[charm.git] / examples / charm++ / Molecular2D / Patch.C
1 /** \file Patch.C
2  *  Author: Abhinav S Bhatele
3  *  Date Created: July 1st, 2008
4  */
5
6 #include "time.h"
7 #include "common.h"
8 #ifdef RUN_LIVEVIZ
9   #include "liveViz.h"
10 #endif
11 #include "Patch.decl.h"
12 #include "Patch.h"
13 #include "Compute.h"
14
15 /* readonly */ CProxy_Main mainProxy;
16 /* readonly */ CProxy_Patch patchArray;
17 /* readonly */ CProxy_Compute computeArray;
18
19 /* readonly */ int numParts;
20 /* readonly */ int patchArrayDimX;      // Number of Chare Rows
21 /* readonly */ int patchArrayDimY;      // Number of Chare Columns
22 /* readonly */ int patchSize; 
23 /* readonly */ double radius;
24 /* readonly */ int finalStepCount; 
25 /* readonly */ double stepTime; 
26
27 double A = 2.0;                 // Force Calculation parameter 1
28 double B = 1.0;                 // Force Calculation parameter 2
29
30 // Entry point of Charm++ application
31 Main::Main(CkArgMsg* msg) {
32   stepTime = CkWallTimer();
33   CkPrintf("\nLENNARD JONES MOLECULAR DYNAMICS RUNNING ...\n");
34
35   numParts = DEFAULT_PARTICLES;
36   patchArrayDimX = PATCHARRAY_DIM_X;
37   patchArrayDimY = PATCHARRAY_DIM_Y;
38   patchSize = PATCH_SIZE;
39   radius = DEFAULT_RADIUS;
40   finalStepCount = DEFAULT_FINALSTEPCOUNT;
41
42   delete msg;
43   mainProxy = thisProxy;
44
45   // initializing the cell 2D array
46   patchArray = CProxy_Patch::ckNew(patchArrayDimX, patchArrayDimY);
47   CkPrintf("%d PATCHES CREATED\n", patchArrayDimX*patchArrayDimY);
48
49   // initializing the interaction 4D array
50   computeArray = CProxy_Compute::ckNew();
51  
52   for (int x=0; x<patchArrayDimX; x++)
53     for (int y=0; y<patchArrayDimY; y++)
54       patchArray(x, y).createComputes();
55
56 }
57
58 // Constructor for chare object migration
59 Main::Main(CkMigrateMessage* msg) { }
60
61 void Main::allDone() {
62   CkPrintf("SIMULATION COMPLETE.\n\n");
63   CkExit();
64 }
65
66 void Main::computeCreationDone() {
67   computeArray.doneInserting();
68   CkPrintf("%d COMPUTES CREATED\n", 5*patchArrayDimX*patchArrayDimY);
69
70 #ifdef RUN_LIVEVIZ
71   // setup liveviz
72   CkCallback c(CkIndex_Patch::requestNextFrame(0), patchArray);
73   liveVizConfig cfg(liveVizConfig::pix_color,true);
74   liveVizInit(cfg,patchArray,c);
75 #endif
76
77   patchArray.start();
78 }
79
80 // Default constructor
81 Patch::Patch() {
82   int i;
83
84   // starting random generator
85   srand48( thisIndex.x * 1000 + thisIndex.y +time(NULL));
86
87   // Particle initialization
88   for(i=0; i < numParts/(patchArrayDimX*patchArrayDimY); i++) {
89     particles.push_back(Particle());
90
91     particles[i].x = drand48() * patchSize + thisIndex.x * patchSize;
92     particles[i].y = drand48() * patchSize + thisIndex.y * patchSize;
93     particles[i].vx = (drand48() - 0.5) * .2 * MAX_VELOCITY;
94     particles[i].vy = (drand48() - 0.5) * .2 * MAX_VELOCITY;
95     particles[i].id = (thisIndex.x*patchArrayDimX + thisIndex.y) * numParts / (patchArrayDimX*patchArrayDimY)  + i;
96   }     
97
98   updateCount = 0;
99   forceCount = 0;
100   stepCount = 0;
101   updateFlag = false;
102   incomingFlag = false;
103   incomingParticles.resize(0);
104 }
105
106 // Constructor for chare object migration
107 Patch::Patch(CkMigrateMessage *msg) { }  
108                                        
109 Patch::~Patch() {}
110
111 void Patch::createComputes() {
112   int i, j, k, l, num;  
113   
114   int x = thisIndex.x;
115   int y = thisIndex.y;
116   int px1, py1, dx, dy, px2, py2;
117
118   // For Round Robin insertion
119   int numPes = CkNumPes();
120   int currPE = CkMyPe();
121  
122   /*  The computes X are inserted by a given patch:
123    *
124    *    ^  X  X  X
125    *    |  0  X  X
126    *    y  0  0  0
127    *       x ---->
128    */
129
130   // these computes will be created by other patches
131   for (num=0; num<NUM_NEIGHBORS; num++) {
132     dx = num/NBRS_Y - NBRS_X/2;
133     dy = num%NBRS_Y - NBRS_Y/2;
134     if (dx == 0) {
135       px1 = px2 = x;
136       if (dy == 0) { py1 = py2 = y; }
137       if (dy > 0) { (y >= patchArrayDimY - NBRS_Y/2) ? ( py1 = WRAP_Y(y+dy), py2 = y ) : ( py1 = y, py2 = y+dy ); }
138       if (dy < 0) { (y < NBRS_Y/2) ? ( py1 = y, py2 = WRAP_Y(y+dy) ) : ( py1 = y+dy, py2 = y ); }
139     }
140
141     if (dx > 0) {
142       (x >= patchArrayDimX - NBRS_X/2) ? 
143       ( px1 = WRAP_X(x+dx), py1 = WRAP_Y(y+dy), px2 = x, py2 = y ) : 
144       ( px1 = x, py1 = y, px2 = WRAP_X(x+dx), py2 = WRAP_Y(y+dy) ) ;
145     }
146
147     if (dx < 0) {
148       (x < NBRS_X/2) ? 
149       ( px1 = x, py1 = y, px2 = WRAP_X(x+dx), py2 = WRAP_Y(y+dy) ) :
150       ( px1 = WRAP_X(x+dx), py1 = WRAP_Y(y+dy), px2 = x, py2 = y ) ;
151     }
152
153     computesList[num][0] = px1; computesList[num][1] = py1; computesList[num][2] = px2; computesList[num][3] = py2;
154
155     //insert only the upper right half computes
156     if (num >= NUM_NEIGHBORS/2)
157       computeArray(px1, py1, px2, py2).insert((currPE++) % numPes);
158   } // end of for loop
159
160   contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::computeCreationDone(), mainProxy));
161 }
162
163 // Function to start interaction among particles in neighboring cells as well as its own particles
164 void Patch::start() {
165   int x = thisIndex.x;
166   int y = thisIndex.y;
167   int i, j, k, l;
168  
169   for(int num=0; num<NUM_NEIGHBORS; num++) {
170     i = computesList[num][0];
171     j = computesList[num][1];
172     k = computesList[num][2];
173     l = computesList[num][3];
174     computeArray(i, j, k, l).interact(particles, x, y);
175   }
176 }
177
178 // Function to update forces coming from a compute
179 void Patch::updateForces(CkVec<Particle> &updates) {
180   int i, x, y, x1, y1;
181   CkVec<Particle> outgoing[NUM_NEIGHBORS];
182
183   // incrementing the counter for receiving updates
184   forceCount++;
185
186   // updating force information
187   for(i = 0; i < updates.length(); i++){
188     particles[i].fx += updates[i].fx;
189     particles[i].fy += updates[i].fy;
190   }
191
192   // if all forces are received, then it must recompute particles location
193   if( forceCount == NUM_NEIGHBORS) {
194     // Received all it's forces from the interactions.
195     forceCount = 0;
196   
197     // Update properties on own particles
198     updateProperties();
199
200     // sending particles to neighboring cells
201     x = thisIndex.x;
202     y = thisIndex.y;
203
204     for(i=0; i<particles.length(); i++) {
205       migrateToPatch(particles[i], x1, y1);
206       if(x1 !=0 || y1!=0) {
207         outgoing[(x1+1)*NBRS_Y + (y1+1)].push_back(wrapAround(particles[i]));
208         particles.remove(i);
209       }
210     }
211    
212     for(i=0; i<NUM_NEIGHBORS; i++)
213       patchArray(WRAP_X(x + i/NBRS_Y - NBRS_X/2), WRAP_Y(y + i%NBRS_Y - NBRS_Y/2)).updateParticles(outgoing[i]);
214
215     updateFlag = true;
216               
217     // checking whether to proceed with next step
218     checkNextStep();
219   }
220   
221 }
222
223 void Patch::migrateToPatch(Particle p, int &px, int &py) {
224   // currently this is assuming that particles are
225   // migrating only to the immediate neighbors
226   int x = thisIndex.x * patchSize;
227   int y = thisIndex.y * patchSize;
228
229   if (p.x < x) px = -1;
230   else if (p.x > x+patchSize) px = 1;
231   else px = 0;
232
233   if (p.y < y) py = -1;
234   else if (p.y > y+patchSize) py = 1;
235   else py = 0;
236 }
237
238 // Function that checks whether it must start the following step or wait until other messages are received
239 void Patch::checkNextStep(){
240   int i;
241   double timer;
242
243   if (updateFlag && incomingFlag) {
244     // resetting flags
245     updateFlag = false;
246     incomingFlag = false;
247     stepCount++;
248
249     // adding new elements
250     for (i = 0; i < incomingParticles.length(); i++)
251       particles.push_back(incomingParticles[i]);
252     incomingParticles.removeAll();
253
254     if (thisIndex.x==0 && thisIndex.y==0 && stepCount%10==0) {
255       timer = CkWallTimer();
256       CkPrintf("Step %d Benchmark Time %f ms/step\n", stepCount, ((timer - stepTime)/10)*1000);
257       stepTime = timer;
258     }
259
260     // checking for next step
261     if (stepCount >= finalStepCount) {
262       print();
263       contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::allDone(), mainProxy)); 
264     } else {
265       thisProxy(thisIndex.x, thisIndex.y).start();
266     }
267   }
268 }
269
270 // Function that receives a set of particles and updates the 
271 // forces of them into the local set
272 void Patch::updateParticles(CkVec<Particle> &updates) {
273   updateCount++;
274
275   for( int i=0; i < updates.length(); i++) {
276     incomingParticles.push_back(updates[i]);
277   }
278
279   // if all the incoming particle updates have been received, we must check 
280   // whether to proceed with next step
281   if(updateCount == NUM_NEIGHBORS-1 ) {
282     updateCount = 0;
283     incomingFlag = true;
284     checkNextStep();
285   }
286 }
287
288 // Function to update properties (i.e. acceleration, velocity and position) in particles
289 void Patch::updateProperties() {
290   int i;
291   double xDisp, yDisp;
292         
293   for(i = 0; i < particles.length(); i++) {
294     // applying kinetic equations
295     particles[i].ax = particles[i].fx / DEFAULT_MASS;
296     particles[i].ay = particles[i].fy / DEFAULT_MASS;
297     particles[i].vx = particles[i].vx + particles[i].ax * DEFAULT_DELTA;
298     particles[i].vy = particles[i].vy + particles[i].ay * DEFAULT_DELTA;
299
300     limitVelocity( particles[i] );
301
302     particles[i].x = particles[i].x + particles[i].vx * DEFAULT_DELTA;
303     particles[i].y = particles[i].y + particles[i].vy * DEFAULT_DELTA;
304
305     particles[i].fx = 0.0;
306     particles[i].fy = 0.0;
307   }
308 }
309
310 void Patch::limitVelocity(Particle &p) {
311   if( fabs( p.vx ) > MAX_VELOCITY ) {
312     if( p.vx < 0.0 )
313       p.vx = -MAX_VELOCITY;
314     else
315       p.vx = MAX_VELOCITY;
316   }
317
318   if( fabs(p.vy) > MAX_VELOCITY ) {
319     if( p.vy < 0.0 )
320       p.vy = -MAX_VELOCITY;
321     else
322       p.vy = MAX_VELOCITY;
323   }
324 }
325
326 Particle& Patch::wrapAround(Particle &p) {
327   if(p.x < 0.0) p.x += patchSize*patchArrayDimX;
328   if(p.y < 0.0) p.y += patchSize*patchArrayDimY;
329   if(p.x > patchSize*patchArrayDimX) p.x -= patchSize*patchArrayDimX;
330   if(p.y > patchSize*patchArrayDimY) p.y -= patchSize*patchArrayDimY;
331
332   return p;
333 }
334
335 // Helper function to help with LiveViz
336 void color_pixel(unsigned char*buf,int width, int height, int xpos,int ypos,
337                              unsigned char R,unsigned char G,unsigned char B) {
338   if(xpos>=0 && xpos<width && ypos>=0 && ypos<height) {
339     buf[3*(ypos*width+xpos)+0] = R; 
340     buf[3*(ypos*width+xpos)+1] = G; 
341     buf[3*(ypos*width+xpos)+2] = B; 
342   }
343 }
344     
345 #ifdef RUN_LIVEVIZ
346 // Each chare provides its particle data to LiveViz
347 void Patch::requestNextFrame(liveVizRequestMsg *lvmsg) {
348   // These specify the desired total image size requested by the client viewer
349   int wdes = lvmsg->req.wid;
350   int hdes = lvmsg->req.ht;
351    
352   int myWidthPx = wdes / patchArrayDimX;
353   int myHeightPx = hdes / patchArrayDimY;
354   int sx=thisIndex.x*myWidthPx;
355   int sy=thisIndex.y*myHeightPx; 
356
357   // set the output pixel values for rectangle
358   // Each component is a char which can have 256 possible values
359   unsigned char *intensity= new unsigned char[3*myWidthPx*myHeightPx];
360   for(int i=0; i<myHeightPx; ++i)
361     for(int j=0; j<myWidthPx; ++j)
362       color_pixel(intensity,myWidthPx,myHeightPx,j,i,0,0,0);    // black background
363
364   for (int i=0; i < particles.length(); i++ ) {
365     int xpos = (int)((particles[i].x /(double) (patchSize*patchArrayDimX)) * wdes) - sx;
366     int ypos = (int)((particles[i].y /(double) (patchSize*patchArrayDimY)) * hdes) - sy;
367
368     Color c(particles[i].id);
369     color_pixel(intensity,myWidthPx,myHeightPx,xpos+1,ypos,c.R,c.B,c.G);
370     color_pixel(intensity,myWidthPx,myHeightPx,xpos-1,ypos,c.R,c.B,c.G);
371     color_pixel(intensity,myWidthPx,myHeightPx,xpos,ypos+1,c.R,c.B,c.G);
372     color_pixel(intensity,myWidthPx,myHeightPx,xpos,ypos-1,c.R,c.B,c.G);
373     color_pixel(intensity,myWidthPx,myHeightPx,xpos,ypos,c.R,c.B,c.G);
374   }
375         
376   for(int i=0; i<myHeightPx; ++i)
377     for(int j=0; j<myWidthPx; ++j) {
378       // Draw red lines
379       if(i==0 || j==0) {
380         color_pixel(intensity,myWidthPx,myHeightPx,j,i,128,0,0);
381       }
382     }
383         
384   liveVizDeposit(lvmsg, sx,sy, myWidthPx,myHeightPx, intensity, this, max_image_data);
385   delete[] intensity;
386 }
387 #endif
388
389 // Prints all particles 
390 void Patch::print(){
391 #ifdef PRINT
392   int i;
393   CkPrintf("*****************************************************\n");
394   CkPrintf("Patch (%d, %d)\n", thisIndex.x, thisIndex.y);
395
396   for(i=0; i < particles.length(); i++)
397     CkPrintf("Patch (%d,%d) %-5d %7.4f %7.4f \n", thisIndex.x, thisIndex.y, i, particles[i].x, particles[i].y);
398   CkPrintf("*****************************************************\n");
399 #endif
400 }
401
402 #include "Patch.def.h"