working version. needs small amount of placing tweaking to put
[charm.git] / examples / charm++ / Molecular / cell.C
1 /*
2  University of Illinois at Urbana-Champaign
3  Department of Computer Science
4  Parallel Programming Lab
5  2008
6 */
7
8 #include "liveViz.h"
9 #include "common.h"
10 #include "cell.decl.h"
11 #include "cell.h"
12 #include "main.decl.h"
13 #include "time.h"
14
15 extern /* readonly */ CProxy_Main mainProxy;
16 extern /* readonly */ CProxy_Cell cellArray;
17 extern /* readonly */ CProxy_Interaction interactionArray;
18
19 extern /* readonly */ int numParts;
20 extern /* readonly */ int m; // Number of Chare Rows
21 extern /* readonly */ int n; // Number of Chare Columns
22 extern /* readonly */ int L; 
23 extern /* readonly */ double radius;
24 extern /* readonly */ int finalStepCount; 
25
26 double A = 2.0; // Force Calculation parameter 1
27 double B = 1.0; // Force Calculation parameter 2
28
29 // Default constructor
30 Cell::Cell() {
31         int i;
32  
33         // starting random generator
34         srand48( thisIndex.x * 1000 + thisIndex.y +time(NULL));
35
36         /* Particle initialization */
37         for(i = 0; i < numParts / (m * n); i++){
38             particles.push_back(Particle());
39
40             particles[i].x = drand48() * L + thisIndex.x * L;
41             particles[i].y = drand48() * L + thisIndex.y * L;
42             particles[i].vx = (drand48() - 0.5) * .2 * MAX_VELOCITY;
43             particles[i].vy = (drand48() - 0.5) * .2 * MAX_VELOCITY;
44             particles[i].id = (thisIndex.x*m + thisIndex.y) * numParts / (m*n)  + i;
45 /*
46             particles[i].x = 0.0 + thisIndex.x * L;
47             particles[i].y = (float) i* 0.9 * L + thisIndex.y * L;
48             particles[i].vx = (drand48() - 0.5) * .2 * MAX_VELOCITY;
49             particles[i].vy = 0.0;
50             particles[i].id = (thisIndex.x*m + thisIndex.y) * numParts / (m*n)  + i;
51 */
52         }       
53
54   updateCount = 0;
55   forceCount = 0;
56   stepCount = 0;
57         updateFlag = false;
58         incomingFlag = false;
59   incomingParticles.resize(0);
60
61 }
62
63 // Constructor for chare object migration
64 Cell::Cell(CkMigrateMessage *msg) { }  
65                                        
66 Cell::~Cell() {}
67
68 // Function to start interaction among particles in neighboring cells as well as its own particles
69 void Cell::start() {
70
71   int x = thisIndex.x;
72   int y = thisIndex.y;
73   int i, j, k, l;
74   
75   // self interaction
76   interactionArray( x, y, x, y).interact(particles, x, y);
77
78   // interaction with (x-1, y-1)
79   (x == 0) ? ( i=x, k=(x-1+m)%m, j=y, l=(y-1+n)%n ) : (i=x-1, k=x, j=(y-1+n)%n, l=y);
80   interactionArray( i, j, k, l ).interact(particles, x, y);
81
82   // interaction with (x-1, y)
83   (x == 0) ? (i=x, k=(x-1+m)%m) : (i=x-1, k=x);
84   interactionArray( i, y, k, y).interact(particles, x, y);
85
86   // interaction with (x-1, y+1)
87   (x == 0) ? ( i=x, k=(x-1+m)%m, j=y, l=(y+1)%n ) : (i=x-1, k=x, j=(y+1)%n, l=y);
88   interactionArray( i, j, k, l ).interact(particles, x, y);
89
90   // interaction with (x, y-1)
91   (y == 0) ? (j=y, l=(y-1+n)%n) : (j=y-1, l=y);
92   interactionArray( x, j, x, l ).interact(particles, x, y);
93
94   // interaction with (x, y+1)
95   (y == n-1) ? (j=(y+1)%n, l=y) : (j=y, l=y+1);// compute
96   interactionArray( x, j, x, l ).interact(particles, x, y);
97
98   // interaction with (x+1, y-1)
99   (x == m-1) ? ( i=0, k=x, j=(y-1+n)%n, l=y ) : (i=x, k=x+1, j=y, l=(y-1+n)%n );
100   interactionArray( i, j, k, l ).interact(particles, x, y);
101
102   // interaction with (x+1, y)
103   (x == m-1) ? (i=0, k=x) : (i=x, k=x+1);
104   interactionArray( i, y, k, y).interact(particles, x, y);
105
106   // interaction with (x+1, y+1)
107   (x == m-1) ? ( i=0, k=x, j=(y+1)%n, l=y ) : (i=x, k=x+1, j=y, l=(y+1)%n );
108   interactionArray( i, j, k, l ).interact(particles, x, y);
109
110 }
111
112 // Function to update forces coming from a neighbor interaction chare
113 void Cell::updateForces(CkVec<Particle> &updates) {
114         int i, x ,y;
115         CkVec<Particle> outgoing;
116
117         // incrementing the counter for receiving updates
118   forceCount++;
119
120         // updating force information
121         for(i = 0; i < updates.length(); i++){
122                 particles[i].fx += updates[i].fx;
123                 particles[i].fy += updates[i].fy;
124         }
125
126         // if all forces are received, then it must recompute particles location
127   if( forceCount >= 9) {
128     // Received all it's forces from the interactions.
129     forceCount = 0;
130     
131     // Update properties on own particles
132                 updateProperties();
133
134                 // Sending particles to neighboring cells
135                 x = thisIndex.x;
136                 y = thisIndex.y;
137
138                 // particles sent to (x-1,y-1)          
139                 outgoing.removeAll();
140                 i = 0;
141                 while(i < particles.length()){
142                         if(particles[i].x < x*L && particles[i].y < y*L){
143                                 outgoing.push_back(wrapAround(particles[i]));
144                                 particles.remove(i);
145                         }else
146                                 i++;
147                 }
148                 cellArray((x-1+m)%m,(y-1+n)%n).updateParticles(outgoing);
149                       
150                 // particles sent to (x-1,y)            
151                 outgoing.removeAll();
152                 i = 0;
153                 while(i < particles.length()){
154                         if(particles[i].x < x*L && particles[i].y <= (y+1)*L){
155                                 outgoing.push_back(wrapAround(particles[i]));
156                                 particles.remove(i);
157                         }else
158                                 i++;
159                 }
160                 cellArray((x-1+m)%m,y).updateParticles(outgoing);
161
162                 // particles sent to (x-1,y+1)
163                 outgoing.removeAll();
164                 i = 0;
165                 while(i < particles.length()){
166                         if(particles[i].x < x*L && particles[i].y > (y+1)*L){
167                                 outgoing.push_back(wrapAround(particles[i]));
168                                 particles.remove(i);
169                         }else
170                                 i++;
171                 }
172                 cellArray((x-1+m)%m,(y+1)%n).updateParticles(outgoing);
173
174                 // particles sent to (x+1,y-1)
175                 outgoing.removeAll();
176                 i = 0;
177                 while(i < particles.length()){
178                         if(particles[i].x > (x+1)*L && particles[i].y < y*L){
179                                 outgoing.push_back(wrapAround(particles[i]));
180                                 particles.remove(i);
181                         }else
182                                 i++;
183                 }
184                 cellArray((x+1)%m,(y-1+n)%n).updateParticles(outgoing);
185
186                 // particles sent to (x+1,y)
187                 outgoing.removeAll();
188                 i = 0;
189                 while(i < particles.length()){
190                         if(particles[i].x > (x+1)*L && particles[i].y <= (y+1)*L){
191                                 outgoing.push_back(wrapAround(particles[i]));
192                                 particles.remove(i);
193                         }else
194                                 i++;
195                 }
196                 cellArray((x+1)%m,y).updateParticles(outgoing);
197
198                 // particles sent to (x+1,y+1)
199                 outgoing.removeAll();
200                 i = 0;
201                 while(i < particles.length()){
202                         if(particles[i].x > (x+1)*L && particles[i].y > (y+1)*L){
203                                 outgoing.push_back(wrapAround(particles[i]));
204                                 particles.remove(i);
205                         }else
206                                 i++;
207                 }
208                 cellArray((x+1)%m,(y+1)%n).updateParticles(outgoing);
209
210                 // particles sent to (x,y-1)
211                 outgoing.removeAll();
212                 i = 0;
213                 while(i < particles.length()){
214                         if(particles[i].y < y*L){
215                                 outgoing.push_back(wrapAround(particles[i]));
216                                 particles.remove(i);
217                         }else
218                                 i++;
219                 }
220                 cellArray(x,(y-1+n)%n).updateParticles(outgoing);
221
222                 // particles sent to (x,y+1)
223                 outgoing.removeAll();
224                 i = 0;
225                 while(i < particles.length()){
226                         if(particles[i].y > (y+1)*L){
227                                 outgoing.push_back(wrapAround(particles[i]));
228                                 particles.remove(i);
229                         }else
230                                 i++;
231                 }
232                 cellArray(x,(y+1)%n).updateParticles(outgoing);
233     outgoing.removeAll();
234
235                 updateFlag = true;
236                 
237     // checking whether to proceed with next step
238                 checkNextStep();
239
240   }
241     
242 }
243
244 // Function that checks whether it must start the following step or wait until other messages are received
245 void Cell::checkNextStep(){
246         int i;
247
248         if(updateFlag && incomingFlag){
249
250                 // resetting flags
251                 updateFlag = false;
252                 incomingFlag = false;
253     stepCount++;
254
255                 // adding new elements
256                 for(i = 0; i < incomingParticles.length(); i++){
257                         particles.push_back(incomingParticles[i]);
258                 }
259                 incomingParticles.removeAll();
260
261                 // checking for next step
262                 if(stepCount >= finalStepCount) {
263                         print();
264                         mainProxy.checkIn();
265           } else {
266             thisProxy( thisIndex.x, thisIndex.y ).start();
267           }
268         }
269
270 }
271
272 // Function that receives a set of particles and updates the forces of them into the local set
273 void Cell::updateParticles(CkVec<Particle> &updates) {
274
275   updateCount++;
276
277   for( int i=0; i < updates.length(); i++) {
278     incomingParticles.push_back(updates[i]);
279   }
280
281         // if all the incoming particle updates have been received, we must check whether to proceed with next step
282         if(updateCount >= 8 ) {
283                 updateCount = 0;
284                 incomingFlag = true;
285                 checkNextStep();
286         }
287
288 }
289
290 // Function to update properties (i.e. acceleration, velocity and position) in particles
291 void Cell::updateProperties(){
292         int i;
293   double xDisp, yDisp;
294         
295         for(i = 0; i < particles.length(); i++){
296
297                 // applying kinetic equations
298                 particles[i].ax = particles[i].fx / DEFAULT_MASS;
299                 particles[i].ay = particles[i].fy / DEFAULT_MASS;
300                 particles[i].vx = particles[i].vx + particles[i].ax * DEFAULT_DELTA;
301                 particles[i].vy = particles[i].vy + particles[i].ay * DEFAULT_DELTA;
302
303     limitVelocity( particles[i] );
304
305                 particles[i].x = particles[i].x + particles[i].vx * DEFAULT_DELTA;
306                 particles[i].y = particles[i].y + particles[i].vy * DEFAULT_DELTA;
307
308                 particles[i].fx = 0.0;
309                 particles[i].fy = 0.0;
310
311         }
312
313 }
314
315 void Cell::limitVelocity(Particle &p) {
316
317     //if( fabs(p.vx * DEFAULT_DELTA) > DEFAULT_RADIUS ) {
318     if( fabs( p.vx ) > MAX_VELOCITY ) {
319       //if( p.vx * DEFAULT_DELTA < 0.0 )
320       if( p.vx < 0.0 )
321         p.vx = -MAX_VELOCITY;
322       else
323         p.vx = MAX_VELOCITY;
324       
325     }
326
327     //if( fabs(p.vy * DEFAULT_DELTA) > DEFAULT_RADIUS ) {
328     if( fabs(p.vy) > MAX_VELOCITY ) {
329
330       //if( p.vy * DEFAULT_DELTA < 0.0 )
331       if( p.vy < 0.0 )
332         p.vy = -MAX_VELOCITY;
333       else
334         p.vy = MAX_VELOCITY;
335     }
336 }
337
338 Particle& Cell::wrapAround(Particle &p) {
339
340                 if(p.x < 0.0) p.x += L*m;
341                 if(p.y < 0.0) p.y += L*n;
342                 if(p.x > L*m) p.x -= L*m;
343                 if(p.y > L*n) p.y -= L*n;
344
345     return p;
346 }
347
348 // Helper function to help with LiveViz
349 void color_pixel(unsigned char*buf,int width, int height, int xpos,int ypos,unsigned char R,unsigned char G,unsigned char B){
350   if(xpos>=0 && xpos<width && ypos>=0 && ypos<height){
351     buf[3*(ypos*width+xpos)+0] = R; 
352     buf[3*(ypos*width+xpos)+1] = G; 
353     buf[3*(ypos*width+xpos)+2] = B; 
354   }
355 }
356     
357 // Each chare provides its particle data to LiveViz
358 void Cell::requestNextFrame(liveVizRequestMsg *lvmsg) {
359   // These specify the desired total image size requested by the client viewer
360   int wdes = lvmsg->req.wid;
361   int hdes = lvmsg->req.ht;
362    
363   int myWidthPx = wdes / m;
364   int myHeightPx = hdes / n;
365   int sx=thisIndex.x*myWidthPx;
366   int sy=thisIndex.y*myHeightPx; 
367
368   // set the output pixel values for rectangle
369   // Each component is a char which can have 256 possible values
370   unsigned char *intensity= new unsigned char[3*myWidthPx*myHeightPx];
371   for(int i=0;i<myHeightPx;++i){
372     for(int j=0;j<myWidthPx;++j){
373                         
374       // black background
375       color_pixel(intensity,myWidthPx,myHeightPx,j,i,0,0,0);
376
377     } 
378   }
379
380   for (int i=0; i < particles.length(); i++ ){
381     
382     int xpos = (int)((particles[i].x /(double) (L*m)) * wdes) - sx;
383     int ypos = (int)((particles[i].y /(double) (L*n)) * hdes) - sy;
384
385     Color c(particles[i].id);
386     color_pixel(intensity,myWidthPx,myHeightPx,xpos+1,ypos,c.R,c.B,c.G);
387     color_pixel(intensity,myWidthPx,myHeightPx,xpos-1,ypos,c.R,c.B,c.G);
388     color_pixel(intensity,myWidthPx,myHeightPx,xpos,ypos+1,c.R,c.B,c.G);
389     color_pixel(intensity,myWidthPx,myHeightPx,xpos,ypos-1,c.R,c.B,c.G);
390     color_pixel(intensity,myWidthPx,myHeightPx,xpos,ypos,c.R,c.B,c.G);
391   }
392         
393   for(int i=0;i<myHeightPx;++i){
394     for(int j=0;j<myWidthPx;++j){
395       
396     // Draw red lines
397     if(i==0 || j==0){
398       color_pixel(intensity,myWidthPx,myHeightPx,j,i,128,0,0);
399     }
400     }
401   }
402         
403   liveVizDeposit(lvmsg, sx,sy, myWidthPx,myHeightPx, intensity, this, max_image_data);
404   delete[] intensity;
405
406 }
407
408 // Prints all particles 
409 void Cell::print(){
410 #ifdef PRINT
411         int i;
412         CkPrintf("*****************************************************\n");
413         CkPrintf("Cell (%d,%d)\n",thisIndex.x,thisIndex.y);
414         //CkPrintf("Part     x     y\n");
415         for(i = 0; i < particles.length(); i++){
416                 CkPrintf("Cell (%d,%d) %-5d %7.4f %7.4f \n",thisIndex.x,thisIndex.y,i,particles[i].x,particles[i].y);
417         }
418         CkPrintf("*****************************************************\n");
419 #endif
420 }
421
422 /* ------------------ Interaction Methods --------------------- */
423
424 // Interaction - Default constructor
425 Interaction::Interaction() {
426   cellCount = 0;
427   bufferedX = 0;
428   bufferedY = 0;
429 }
430
431 Interaction::Interaction(CkMigrateMessage *msg) { }
432   
433 // Function to receive vector of particles
434 void Interaction::interact(CkVec<Particle> particles, int x, int y ) {
435
436   int i;
437
438   // self interaction check
439   if( thisIndex.x == thisIndex.z && thisIndex.y == thisIndex.w ) {
440                 interact(particles,particles);
441     cellArray( x, y).updateForces(particles);
442   } else {
443     if(cellCount == 0) {
444
445               bufferedX = x;
446           bufferedY = y;
447         bufferedParticles = particles;
448         cellCount++;
449
450                 } else if(cellCount == 1) {
451     
452             // if both particle sets are received, compute interaction
453       cellCount = 0;
454                   interact(bufferedParticles,particles);
455       cellArray(bufferedX, bufferedY).updateForces(bufferedParticles);
456       cellArray(x, y).updateForces(particles);
457
458     }
459
460   }
461 }
462
463 // Function to compute all the interactions between pairs of particles in two sets
464 void Interaction::interact(CkVec<Particle> &first, CkVec<Particle> &second){
465         int i, j;
466         for(i = 0; i < first.length(); i++)
467                 for(j = 0; j < second.length(); j++)
468                         interact(first[i], second[j]);
469 }
470
471 // Function for computing interaction among two particles
472 // There is an extra test for interaction of identical particles, in which case there is no effect
473 void Interaction::interact(Particle &first, Particle &second){
474         float rx,ry,rz,r,fx,fy,fz,f;
475
476         // computing base values
477         rx = first.x - second.x;
478         ry = first.y - second.y;
479         r = sqrt(rx*rx + ry*ry);
480
481   // We include 0.000001 to ensure that r doesn't tend to zero in the force calculation
482         //if(r < 0.000001 || r >= DEFAULT_RADIUS)
483         if(r < 0.000001 || r >= L)
484                 return;
485
486         f = A / pow(r,12) - B / pow(r,6);
487         fx = f * rx / r;
488         fy = f * ry / r;
489
490         // updating particle properties
491         second.fx -= fx;
492         second.fy -= fy;
493         first.fx += fx;
494         first.fy += fy;
495
496 }
497
498 #include "cell.def.h"