*** empty log message ***
[charm.git] / examples / charm++ / Molecular / main.C
1 /*
2  University of Illinois at Urbana-Champaign
3  Department of Computer Science
4  Parallel Programming Lab
5  2008
6  Authors: Kumaresh Pattabiraman and Esteban Meneses
7 */
8
9 #include "main.h"
10 #include "main.decl.h"
11
12 #define DEBUG 1
13
14 #define A 2.0
15 #define B 1.0
16
17 #define DEFAULT_PARTICLES 350
18 #define DEFAULT_M 5
19 #define DEFAULT_N 7
20 #define DEFAULT_L 10
21 #define DEFAULT_RADIUS 10
22 #define DEFAULT_FINALSTEPCOUNT 2
23
24 /* readonly */ CProxy_Main mainProxy;
25 /* readonly */ CProxy_Cell cellArray;
26 /* readonly */ CProxy_Interaction interactionArray;
27
28 /* readonly */ int numParts;
29 /* readonly */ int m; // Number of Chare Rows
30 /* readonly */ int n; // Number of Chare Columns
31 /* readonly */ int L; 
32 /* readonly */ double radius;
33 /* readonly */ int finalStepCount; 
34
35 // Class representing a cell in the grid. We consider each cell as a square of LxL units
36 class Cell : public CBase_Cell {
37   private:
38     CkVec<Particle> particles;
39                 CkVec<Particle> incomingParticles;
40     int forceCount;                                                                                                                                     // to count the returns from interactions
41     int stepCount;                                                                                                                                      // to count the number of steps, and decide when to stop
42                 int updateCount;
43                 bool updateFlag;
44                 bool incomingFlag;
45
46                 void updateProperties();                                                                                                        // updates properties after receiving forces from interactions
47                 void checkNextStep();                                                                                                                   // checks whether to continue with next step
48                 void print();                                                                                                                                                   // prints all its particles
49
50   public:
51     Cell();
52     Cell(CkMigrateMessage *msg);
53     ~Cell();
54
55     void start();
56     void updateParticles(CkVec<Particle>&);
57     void updateForces(CkVec<Particle>&);
58     void stepDone();
59 };
60
61 // Class representing the interaction agents between a couple of cells
62 class Interaction : public CBase_Interaction {
63   private:
64     int cellCount;                                                                                                                                      // to count the number of interact() calls
65     CkVec<Particle> bufferedParticles;
66     int bufferedX;
67                 int bufferedY;
68
69                 void interact(CkVec<Particle> &first, CkVec<Particle> &second);
70                 void interact(Particle &first, Particle &second);
71
72   public:
73     Interaction();
74     Interaction(CkMigrateMessage *msg);
75
76     void interact(CkVec<Particle> particles, int i, int j);
77
78
79 };
80     
81 // Main class
82 class Main : public CBase_Main {
83
84   private:
85     int checkInCount;                                                                                                                           // Count to terminate
86
87     #ifdef DEBUG
88       int interactionCount;
89     #endif
90
91   public:
92
93     /// Constructors ///
94     Main(CkArgMsg* msg);
95     Main(CkMigrateMessage* msg);
96
97     void checkIn();
98 };
99
100 // Entry point of Charm++ application
101 Main::Main(CkArgMsg* msg) {
102         int i, j, k, l;  
103
104   numParts = DEFAULT_PARTICLES;
105   m = DEFAULT_M;
106   n = DEFAULT_N;
107         L = DEFAULT_L;
108   radius = DEFAULT_RADIUS;
109   finalStepCount = DEFAULT_FINALSTEPCOUNT;
110
111   delete msg;
112   checkInCount = 0;
113
114   #ifdef DEBUG
115     interactionCount=0;
116   #endif
117
118   mainProxy = thisProxy;
119
120         // initializing the cell 2D array
121   cellArray = CProxy_Cell::ckNew(m,n);
122
123         // initializing the interaction 4D array
124   interactionArray = CProxy_Interaction::ckNew();
125   
126   for (int x = 0; x < m ; x++ ) {
127     for (int y = 0; y < n; y++ ) {
128
129       //Processor Round Robin needed /* FIXME */
130  
131       #ifdef DEBUG
132         //CkPrintf("INITIAL:( %d, %d) ( %d , %d )\n", x,y,x,y);
133         interactionCount++;
134       #endif
135
136       // self interaction
137       interactionArray( x, y, x, y ).insert( /* processor number */0 );
138
139       // (x,y) and (x+1,y) pair
140       (x == m-1) ? (i=(x+1)%m, k=x) : (i=x, k=x+1);
141       #ifdef DEBUG
142         //CkPrintf("INITIAL:( %d, %d) ( %d , %d )\n", i,y,k,y);
143         interactionCount++;
144       #endif
145       interactionArray( i, y, k, y ).insert( /* processor number */0 );
146
147       // (x,y) and (x,y+1) pair
148       (y == n-1) ? (j=(y+1)%n, l=y) : (j=y, l=y+1);
149       #ifdef DEBUG
150         //CkPrintf("INITIAL:( %d, %d) ( %d , %d )\n", x,j,x,l);
151         interactionCount++;
152       #endif
153
154                         
155       interactionArray( x, j, x, l ).insert( /* processor number */0 );
156
157       // (x,y) and (x+1,y+1) pair, Irrespective of y /* UNDERSTAND */
158       (x == m-1) ? ( i=(x+1)%m, k=x, j=(y+1)%n, l=y ) : (i=x, k=x+1, j=y, l=(y+1)%n );
159       #ifdef DEBUG
160         //CkPrintf("INITIAL:( %d, %d) ( %d , %d )\n", i,j,k,l);
161         interactionCount++;
162       #endif
163       interactionArray( i, j, k, l ).insert( /* processor number */0 );
164
165       // (x,y) and (x-1,y+1) pair /* UNDERSTAND */
166       (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 );
167       #ifdef DEBUG
168         //CkPrintf("INITIAL:( %d, %d) ( %d , %d )\n", i,j,k,l);
169         interactionCount++;
170       #endif
171       interactionArray( i, j, k, l ).insert( /* processor number */0 );
172
173     }
174   }
175
176   interactionArray.doneInserting();
177
178   #ifdef DEBUG
179     CkPrintf("Interaction Count: %d\n", interactionCount);
180   #endif
181
182   cellArray.start();
183 }
184
185 // Constructor needed for chare object migration
186 Main::Main(CkMigrateMessage* msg) { }
187
188 void Main::checkIn() {
189
190   checkInCount ++;
191   if( checkInCount >= m*n)
192     CkExit();
193
194 }
195
196 // Default constructor
197 Cell::Cell() {
198         int i;
199
200         // starting random generator
201         srand48(time(NULL));
202   
203         /* Particle initialization */
204         // initializing a number of particles
205         for(i = 0; i < numParts / (m * n); i++){
206                 particles.push_back(Particle());
207                 particles[i].x = drand48() * L + thisIndex.x * L;
208     particles[i].y = drand48() * L + thisIndex.y * L;
209         }       
210
211   forceCount = 0;
212   stepCount = 0;
213         updateFlag = false;
214         incomingFlag = false;
215
216 }
217
218 // Constructor needed for chare object migration (ignore for now)
219 Cell::Cell(CkMigrateMessage *msg) { }                                         
220 Cell::~Cell() {
221   /* FIXME */ // Deallocate particle lists
222 }
223
224 // Function to start interaction among particles in neighboring cells as well as its own particles
225 void Cell::start() {
226
227   int x = thisIndex.x;
228   int y = thisIndex.y;
229
230   int i, j, k, l;
231
232   #ifdef DEBUG
233     //CkPrintf("START:( %d, %d) ( %d , %d )\n", x,y,x,y);
234   #endif
235   
236   // self interaction
237   interactionArray( x, y, x, y).interact(particles, x, y);
238
239   // interaction with (x-1, y-1)
240   (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);
241   interactionArray( i, j, k, l ).interact(particles, x, y);
242
243   // interaction with (x-1, y)
244   (x == 0) ? (i=x, k=(x-1+m)%m) : (i=x-1, k=x);
245   interactionArray( i, y, k, y).interact(particles, x, y);
246
247   // interaction with (x-1, y+1)
248   (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);
249   interactionArray( i, j, k, l ).interact(particles, x, y);
250
251
252   // interaction with (x, y-1)
253   (y == 0) ? (j=y, l=(y-1+n)%n) : (j=(y-1+n)%n, l=y);
254   interactionArray( x, j, x, l ).interact(particles, x, y);
255
256   // interaction with (x, y+1)
257   (y == n-1) ? (j=(y+1)%n, l=y) : (j=y, l=y+1);// compute
258   interactionArray( x, j, x, l ).interact(particles, x, y);
259
260
261   // interaction with (x+1, y-1)
262   (x == m-1) ? ( i=(x+1)%m, k=x, j=(y-1+n)%n, l=y ) : (i=x, k=x+1, j=y, l=(y-1+n)%n );
263   interactionArray( i, j, k, l ).interact(particles, x, y);
264
265   // interaction with (x+1, y)
266   (x == m-1) ? (i=(x+1)%m, k=x) : (i=x, k=x+1);
267   interactionArray( i, y, k, y).interact(particles, x, y);
268
269   // interaction with (x+1, y+1)
270   (x == m-1) ? ( i=(x+1)%m, k=x, j=(y+1)%n, l=y ) : (i=x, k=x+1, j=y, l=(y+1)%n );
271   interactionArray( i, j, k, l ).interact(particles, x, y);
272
273 }
274
275 // Function to update forces coming from a neighbor interaction chare
276 void Cell::updateForces(CkVec<Particle> &updates) {
277         int i, x ,y;
278         CkVec<Particle> outgoing;
279
280         // incrementing the counter for receiving updates
281   forceCount++;
282
283         // updating force information
284         for(i = 0; i < updates.length(); i++){
285                 particles[i].fx += updates[i].fx;
286                 particles[i].fy += updates[i].fy;
287         }
288
289         // if all forces are received, then it must recompute particles location
290   if( forceCount >= 9) {
291     
292                 CkPrintf("Forces done!\n");
293
294     // Received all it's forces from the interactions.
295     stepCount++;
296     forceCount = 0;
297     
298     // Update properties on own particles
299                 updateProperties();
300                 updateFlag = true;
301
302                 // Sending particles to neighboring cells
303                 x = thisIndex.x;
304                 y = thisIndex.y;
305                 // particles sent to (x-1,y-1)          
306                 outgoing.removeAll();
307                 i = 0;
308                 while(i < particles.length()){
309                         if(particles[i].x < x*L && particles[i].y < y*L){
310                                 outgoing.push_back(particles[i]);
311                                 particles.remove(i);
312                         }else
313                                 i++;
314                 }
315                 cellArray((x-1+m)%m,(y-1+n)%n).updateParticles(outgoing);
316                       
317                 // particles sent to (x-1,y)            
318                 outgoing.removeAll();
319                 i = 0;
320                 while(i < particles.length()){
321                         if(particles[i].x < x*L && particles[i].y <= y*(L+1)){
322                                 outgoing.push_back(particles[i]);
323                                 particles.remove(i);
324                         }else
325                                 i++;
326                 }
327                 cellArray((x-1+m)%m,y).updateParticles(outgoing);
328
329                 // particles sent to (x-1,y+1)
330                 outgoing.removeAll();
331                 i = 0;
332                 while(i < particles.length()){
333                         if(particles[i].x < x*L && particles[i].y > y*(L+1)){
334                                 outgoing.push_back(particles[i]);
335                                 particles.remove(i);
336                         }else
337                                 i++;
338                 }
339                 cellArray((x-1+m)%m,(y+1)%n).updateParticles(outgoing);
340
341                 // particles sent to (x+1,y-1)
342                 outgoing.removeAll();
343                 i = 0;
344                 while(i < particles.length()){
345                         if(particles[i].x > x*(L+1) && particles[i].y < y*L){
346                                 outgoing.push_back(particles[i]);
347                                 particles.remove(i);
348                         }else
349                                 i++;
350                 }
351                 cellArray((x+1)%m,(y-1+n)%n).updateParticles(outgoing);
352
353                 // particles sent to (x+1,y)
354                 outgoing.removeAll();
355                 i = 0;
356                 while(i < particles.length()){
357                         if(particles[i].x > x*(L+1) && particles[i].y <= y*(L+1)){
358                                 outgoing.push_back(particles[i]);
359                                 particles.remove(i);
360                         }else
361                                 i++;
362                 }
363                 cellArray((x+1)%m,y).updateParticles(outgoing);
364
365                 // particles sent to (x+1,y+1)
366                 outgoing.removeAll();
367                 i = 0;
368                 while(i < particles.length()){
369                         if(particles[i].x > x*(L+1) && particles[i].y > y*(L+1)){
370                                 outgoing.push_back(particles[i]);
371                                 particles.remove(i);
372                         }else
373                                 i++;
374                 }
375                 cellArray((x+1)%m,(y+1)%n).updateParticles(outgoing);
376
377                 // particles sent to (x,y-1)
378                 outgoing.removeAll();
379                 i = 0;
380                 while(i < particles.length()){
381                         if(particles[i].y < y*L){
382                                 outgoing.push_back(particles[i]);
383                                 particles.remove(i);
384                         }else
385                                 i++;
386                 }
387                 cellArray(x,(y-1+n)%n).updateParticles(outgoing);
388
389                 // particles sent to (x,y+1)
390                 outgoing.removeAll();
391                 i = 0;
392                 while(i < particles.length()){
393                         if(particles[i].y > (y+1)*L){
394                                 outgoing.push_back(particles[i]);
395                                 particles.remove(i);
396                         }else
397                                 i++;
398                 }
399                 cellArray(x,(y+1)%n).updateParticles(outgoing);
400
401     #ifdef DEBUG
402       CkPrintf("STEP: %d DONE:( %d , %d )\n", stepCount, thisIndex.x, thisIndex.y);
403     #endif
404
405                 // checking whether to proceed with next step
406                 checkNextStep();
407
408   }
409     
410 }
411
412 // Prints all particle set
413 void Cell::print(){
414         int i;
415         CkPrintf("*******************************************************************\n");
416         CkPrintf("Cell (%d,%d)\n",thisIndex.x,thisIndex.y);
417         CkPrintf("Part     x     y\n");
418         for(i = 0; i < particles.length(); i++){
419                 CkPrintf("%-5d %7.4f %7.4f \n",i,particles[i].x,particles[i].y);
420         }
421         CkPrintf("*******************************************************************\n");
422 }
423
424 // Function that checks whether it must start the following step or wait until other messages are received
425 void Cell::checkNextStep(){
426         int i;
427
428         if(updateFlag && incomingFlag){
429
430                 // resetting flags
431                 updateFlag = false;
432                 incomingFlag = false;
433
434                 // adding new elements
435                 for(i = 0; i < incomingParticles.length(); i++){
436                         particles.push_back(incomingParticles[i]);
437                 }
438                 incomingParticles.removeAll();
439
440                 // checking for next step
441                 if(stepCount >= finalStepCount) {
442                         print();
443                         mainProxy.checkIn();
444           } else {
445             thisProxy( thisIndex.x, thisIndex.y ).start();
446           }
447         }
448
449 }
450
451 // Function that receives a set of particles and updates the forces of them into the local set
452 void Cell::updateParticles(CkVec<Particle> &updates) {
453   updateCount++;
454
455   for( int i=0; i < updates.length(); i++) {
456     incomingParticles.push_back(updates[i]);
457   }
458
459         // if all the incoming particle updates have been received, we must check whether to proceed with next step
460         if(updateCount >= 8 ) {
461                 updateCount = 0;
462                 incomingFlag = true;
463                 checkNextStep();
464         }
465
466 }
467
468 // Function to update properties (i.e. acceleration, velocity and position) in particles
469 void Cell::updateProperties(){
470         int i;
471         
472         for(i = 0; i < particles.length(); i++){
473
474                 // applying kinetic equations
475                 particles[i].ax = particles[i].fx / DEFAULT_MASS;
476                 particles[i].ay = particles[i].fy / DEFAULT_MASS;
477                 particles[i].vx = particles[i].vx + particles[i].ax * DEFAULT_DELTA;
478                 particles[i].vy = particles[i].vy + particles[i].ay * DEFAULT_DELTA;
479                 particles[i].x = particles[i].x + particles[i].vx * DEFAULT_DELTA;
480                 particles[i].y = particles[i].y + particles[i].vy * DEFAULT_DELTA;
481                 particles[i].fx = 0.0;
482                 particles[i].fy = 0.0;
483
484                 // checking boundary conditions
485                 if(particles[i].x < 0.0) particles[i].x = 0.0;
486                 if(particles[i].y < 0.0) particles[i].y = 0.0;
487                 if(particles[i].x > L*m) particles[i].x = L*m;
488                 if(particles[i].y > L*n) particles[i].y = L*n;
489
490         }
491
492 }
493
494
495 // Default constructor
496 Interaction::Interaction() {
497   cellCount = 0;
498
499   /* FIXME */
500     bufferedX = 0;
501     bufferedY = 0;
502 }
503
504 Interaction::Interaction(CkMigrateMessage *msg) { }
505   
506
507 // Function to receive vector of particles
508 void Interaction::interact(CkVec<Particle> particles, int x, int y ) {
509
510   if(cellCount == 0) {
511
512     // self interaction check
513     if( thisIndex.x == thisIndex.z && thisIndex.y == thisIndex.w ) {
514       //CkPrintf("SELF: ( %d , %d )\n", thisIndex.x, thisIndex.y );
515       cellCount = 0;
516                         interact(particles,particles);
517       cellArray( x, y).updateForces(particles);
518     } else {
519                          bufferedX = x;
520          bufferedY = y;
521                          bufferedParticles = particles; /* CHECKME */
522                 }
523
524   }
525
526   cellCount++;
527
528         // if both particle sets are received, compute interaction
529   if(cellCount >= 2) {
530
531     //CkPrintf("PAIR:( %d , %d )  ( %d , %d ) \n", bufferedX, bufferedY, x, y );
532     cellCount = 0;
533
534                 interact(bufferedParticles,particles);
535
536     cellArray(bufferedX, bufferedY).updateForces(bufferedParticles);
537     cellArray(x, y).updateForces(particles);
538
539   }
540 }
541
542 // Function to compute all the interactions between pairs of particles in two sets
543 void Interaction::interact(CkVec<Particle> &first, CkVec<Particle> &second){
544         int i, j;
545         for(i = 0; i < first.length(); i++)
546                 for(j = 0; j < second.length(); j++)
547                         interact(first[i], second[j]);
548 }
549
550 // Function for computing interaction among two particles
551 // There is an extra test for interaction of identical particles, in which case there is no effect
552 void Interaction::interact(Particle &first, Particle &second){
553         float rx,ry,rz,r,fx,fy,fz,f;
554
555         // computing base values
556         rx = first.x - second.x;
557         ry = first.y - second.y;
558         r = sqrt(rx*rx + ry*ry);
559
560         if(r == 0.0 || r >= DEFAULT_RADIUS)
561                 return;
562
563         f = A / pow(r,12) - B / pow(r,6);
564         fx = f * rx / r;
565         fy = f * ry / r;
566
567         // updating particle properties
568         second.fx += fx;
569         second.fy += fy;
570         first.fx += -fx;
571         first.fy += -fy;
572
573 }
574
575 #include "main.def.h"