d45af8f275573395bd19947f02cad46bc2923e94
[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 1000
18 #define DEFAULT_M 4
19 #define DEFAULT_N 3
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
44   public:
45     Cell();
46     Cell(CkMigrateMessage *msg);
47     ~Cell();
48
49     void start();
50     void updateParticles(CkVec<Particle>&);
51     void updateForces(CkVec<Particle>&);
52     void stepDone();
53 };
54
55 // Class representing the interaction agents between a couple of cells
56 class Interaction : public CBase_Interaction {
57   private:
58     int cellCount;                                                                                                                                      // to count the number of interact() calls
59     CkVec<Particle> bufferedParticles;
60     int bufferedX;
61                 int bufferedY;
62
63                 void interact(CkVec<Particle> &first, CkVec<Particle> &second);
64                 void interact(Particle &first, Particle &second);
65
66   public:
67     Interaction();
68     Interaction(CkMigrateMessage *msg);
69
70     void interact(CkVec<Particle> particles, int i, int j);
71
72
73 };
74     
75 // Main class
76 class Main : public CBase_Main {
77
78   private:
79     int checkInCount;                                                                                                                           // Count to terminate
80
81     #ifdef DEBUG
82       int interactionCount;
83     #endif
84
85   public:
86
87     /// Constructors ///
88     Main(CkArgMsg* msg);
89     Main(CkMigrateMessage* msg);
90
91     void checkIn();
92 };
93
94 // Entry point of Charm++ application
95 Main::Main(CkArgMsg* msg) {
96         int i, j, k, l;  
97
98   numParts = DEFAULT_PARTICLES;
99   m = DEFAULT_M;
100   n = DEFAULT_N;
101         L = DEFAULT_L;
102   radius = DEFAULT_RADIUS;
103   finalStepCount = DEFAULT_FINALSTEPCOUNT;
104
105   delete msg;
106   checkInCount = 0;
107
108   #ifdef DEBUG
109     interactionCount=0;
110   #endif
111
112   mainProxy = thisProxy;
113
114         // initializing the cell 2D array
115   cellArray = CProxy_Cell::ckNew(m,n);
116
117         // initializing the interaction 4D array
118   interactionArray = CProxy_Interaction::ckNew();
119   
120   for (int x = 0; x < m ; x++ ) {
121     for (int y = 0; y < n; y++ ) {
122
123       //Processor Round Robin needed /* FIXME */
124  
125       #ifdef DEBUG
126         CkPrintf("INITIAL:( %d, %d) ( %d , %d )\n", x,y,x,y);
127         interactionCount++;
128       #endif
129
130       // self interaction
131       interactionArray( x, y, x, y ).insert( /* processor number */0 );
132
133       // (x,y) and (x+1,y) pair
134       (x == m-1) ? (i=(x+1)%m, k=x) : (i=x, k=x+1);
135       #ifdef DEBUG
136         CkPrintf("INITIAL:( %d, %d) ( %d , %d )\n", i,y,k,y);
137         interactionCount++;
138       #endif
139       interactionArray( i, y, k, y ).insert( /* processor number */0 );
140
141       // (x,y) and (x,y+1) pair
142       (y == n-1) ? (j=(y+1)%n, l=y) : (j=y, l=y+1);
143       #ifdef DEBUG
144         CkPrintf("INITIAL:( %d, %d) ( %d , %d )\n", x,j,x,l);
145         interactionCount++;
146       #endif
147
148                         
149       interactionArray( x, j, x, l ).insert( /* processor number */0 );
150
151       // (x,y) and (x+1,y+1) pair, Irrespective of y /* UNDERSTAND */
152       (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 );
153       #ifdef DEBUG
154         CkPrintf("INITIAL:( %d, %d) ( %d , %d )\n", i,j,k,l);
155         interactionCount++;
156       #endif
157       interactionArray( i, j, k, l ).insert( /* processor number */0 );
158
159       // (x,y) and (x-1,y+1) pair /* UNDERSTAND */
160       (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 );
161       #ifdef DEBUG
162         CkPrintf("INITIAL:( %d, %d) ( %d , %d )\n", i,j,k,l);
163         interactionCount++;
164       #endif
165       interactionArray( i, j, k, l ).insert( /* processor number */0 );
166
167     }
168   }
169
170   interactionArray.doneInserting();
171
172   #ifdef DEBUG
173     CkPrintf("Interaction Count: %d\n", interactionCount);
174   #endif
175
176   cellArray.start();
177 }
178
179 // Constructor needed for chare object migration
180 Main::Main(CkMigrateMessage* msg) { }
181
182 void Main::checkIn() {
183
184   checkInCount ++;
185   if( checkInCount >= m*n)
186     CkExit();
187
188 }
189
190 // Default constructor
191 Cell::Cell() {
192         int i;
193
194         // starting random generator
195         srand48(time(NULL));
196   
197         /* Particle initialization */
198         // initializing a number of particles
199         for(i = 0; i < numParts / (m * n); i++){
200                 particles.push_back(Particle());
201                 particles[i].x = drand48() * L + thisIndex.x * L;
202     particles[i].y = drand48() * L + thisIndex.y * L;
203         }       
204
205   forceCount = 0;
206   stepCount = 0;
207 }
208
209 // Constructor needed for chare object migration (ignore for now)
210 Cell::Cell(CkMigrateMessage *msg) { }                                         
211 Cell::~Cell() {
212   /* FIXME */ // Deallocate particle lists
213 }
214
215
216 void Cell::start() {
217
218   int x = thisIndex.x;
219   int y = thisIndex.y;
220
221   int i, j, k, l;
222
223   #ifdef DEBUG
224     CkPrintf("START:( %d, %d) ( %d , %d )\n", x,y,x,y);
225   #endif
226   
227   // self interaction
228   interactionArray( x, y, x, y).interact(particles, x, y);
229
230   // interaction with (x-1, y-1)
231   (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);
232   interactionArray( i, j, k, l ).interact(particles, x, y);
233
234   // interaction with (x-1, y)
235   (x == 0) ? (i=x, k=(x-1+m)%m) : (i=x-1, k=x);
236   interactionArray( i, y, k, y).interact(particles, x, y);
237
238   // interaction with (x-1, y+1)
239   (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);
240   interactionArray( i, j, k, l ).interact(particles, x, y);
241
242
243   // interaction with (x, y-1)
244   (y == 0) ? (j=y, l=(y-1+n)%n) : (j=(y-1+n)%n, l=y);
245   interactionArray( x, j, x, l ).interact(particles, x, y);
246
247   // interaction with (x, y+1)
248   (y == n-1) ? (j=(y+1)%n, l=y) : (j=y, l=y+1);// compute
249   interactionArray( x, j, x, l ).interact(particles, x, y);
250
251
252   // interaction with (x+1, y-1)
253   (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 );
254   interactionArray( i, j, k, l ).interact(particles, x, y);
255
256   // interaction with (x+1, y)
257   (x == m-1) ? (i=(x+1)%m, k=x) : (i=x, k=x+1);
258   interactionArray( i, y, k, y).interact(particles, x, y);
259
260   // interaction with (x+1, y+1)
261   (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 );
262   interactionArray( i, j, k, l ).interact(particles, x, y);
263
264 }
265
266 void Cell::updateForces(CkVec<Particle> &particles) {
267   forceCount++;
268   if( forceCount >= 9) {
269     
270     // Received all it's forces from the interactions.
271     stepCount++;
272     forceCount = 0;
273     
274     /* FIX ME*/ // Update forces on own particles and determine which of them are being displaced to neighbors
275     /* FIX ME*/ // Calls to updateParticles on neighbours.
276       
277     #ifdef DEBUG
278       CkPrintf("STEP: %d DONE:( %d , %d )\n", stepCount, thisIndex.x, thisIndex.y);
279     #endif
280
281     if(stepCount >= finalStepCount) {
282       mainProxy.checkIn();
283     } else {
284       thisProxy( thisIndex.x, thisIndex.y ).start();
285     }
286   }
287     
288 }
289
290 // Function that receives a set of particles and updates the forces of them into the local set
291 void Cell::updateParticles(CkVec<Particle> &updates) {
292
293   updateCount++;
294
295   for( int i=0; i < updates.length(); i++) {
296     incomingParticles.push_back(updates[i]);
297   }
298
299
300   
301         //CHECKincomingParticles.append(updates);
302
303   
304
305         if(updateCount >= 8 ) {
306                 
307                 /* FIXME Synchronisation?? */
308         
309   
310     updateCount = 0;
311     /* FIXME Empty incomingParticles vector after appending it with particles */
312         }
313
314 }
315
316 // Default constructor
317 Interaction::Interaction() {
318   cellCount = 0;
319
320   /* FIXME */
321     bufferedX = 0;
322     bufferedY = 0;
323 }
324
325 Interaction::Interaction(CkMigrateMessage *msg) { }
326   
327
328 // Function to receive vector of particles
329 void Interaction::interact(CkVec<Particle> particles, int x, int y ) {
330
331   if(cellCount == 0) {
332
333     // self interaction check
334     if( thisIndex.x == thisIndex.z && thisIndex.y == thisIndex.w ) {
335       CkPrintf("SELF: ( %d , %d )\n", thisIndex.x, thisIndex.y );
336       cellCount = 0;
337                         interact(particles,particles);
338       cellArray( x, y).updateForces(particles);
339     } else {
340                          bufferedX = x;
341          bufferedY = y;
342                          bufferedParticles = particles; /* CHECKME */
343                 }
344
345   }
346
347   cellCount++;
348
349         // if both particle sets are received, compute interaction
350   if(cellCount >= 2) {
351
352     CkPrintf("PAIR:( %d , %d )  ( %d , %d ) \n", bufferedX, bufferedY, x, y );
353     cellCount = 0;
354
355                 interact(bufferedParticles,particles);
356
357     cellArray(bufferedX, bufferedY).updateForces(bufferedParticles);
358     cellArray(x, y).updateForces(particles);
359
360   }
361 }
362
363 // Function to compute all the interactions between pairs of particles in two sets
364 void Interaction::interact(CkVec<Particle> &first, CkVec<Particle> &second){
365         int i, j;
366         for(i = 0; i < first.length(); i++)
367                 for(j = 0; j < second.length(); j++)
368                         interact(first[i], second[j]);
369 }
370
371 // Function for computing interaction among two particles
372 // There is an extra test for interaction of identical particles, in which case there is no effect
373 void Interaction::interact(Particle &first, Particle &second){
374         float rx,ry,rz,r,fx,fy,fz,f;
375
376         // computing base values
377         rx = first.x - second.x;
378         ry = first.y - second.y;
379         r = sqrt(rx*rx + ry*ry);
380
381         if(r == 0.0 || r >= DEFAULT_RADIUS)
382                 return;
383
384         f = A / pow(r,12) - B / pow(r,6);
385         fx = f * rx / r;
386         fy = f * ry / r;
387
388         // updating particle properties
389         second.fx += fx;
390         second.fy += fy;
391         first.fx += -fx;
392         first.fy += -fy;
393
394 }
395
396 #include "main.def.h"