Several features regarding Particle were included.
[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; /* CHECKME */
26 /* readonly */ CProxy_Interaction interactionArray; /* CHECKME */
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> incomingParts;
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> updates);
51     void force(CkVec<Particle> particles);
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         // initializing a number of particles
198         for(i = 0; i < numParts / (m * n); i++){
199                 particles.push_back(Particle());
200                 particles[i].x = drand48() * L + thisIndex.x * L;
201     particles[i].y = drand48() * L + thisIndex.y * L;
202         }       
203
204         /* Particle initialization */
205
206   forceCount = 0;
207   stepCount = 0;
208 }
209
210 // Constructor needed for chare object migration (ignore for now)
211 // NOTE: This constructor does not need to appear in the ".ci" file
212 Cell::Cell(CkMigrateMessage *msg) { }                                         
213 Cell::~Cell() {
214   /* FIXME */ // Deallocate particle lists
215 }
216
217
218 void Cell::start() {
219
220   int x = thisIndex.x;
221   int y = thisIndex.y;
222
223   int i, j, k, l;
224
225   #ifdef DEBUG
226     CkPrintf("START:( %d, %d) ( %d , %d )\n", x,y,x,y);
227   #endif
228   
229   // self interaction
230   interactionArray( x, y, x, y).interact(particles, x, y);
231
232   // interaction with (x-1, y-1)
233   (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);
234   interactionArray( i, j, k, l ).interact(particles, x, y);
235
236   // interaction with (x-1, y)
237   (x == 0) ? (i=x, k=(x-1+m)%m) : (i=x-1, k=x);
238   interactionArray( i, y, k, y).interact(particles, x, y);
239
240   // interaction with (x-1, y+1)
241   (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);
242   interactionArray( i, j, k, l ).interact(particles, x, y);
243
244
245   // interaction with (x, y-1)
246   (y == 0) ? (j=y, l=(y-1+n)%n) : (j=(y-1+n)%n, l=y);
247   interactionArray( x, j, x, l ).interact(particles, x, y);
248
249   // interaction with (x, y+1)
250   (y == n-1) ? (j=(y+1)%n, l=y) : (j=y, l=y+1);// compute
251   interactionArray( x, j, x, l ).interact(particles, x, y);
252
253
254   // interaction with (x+1, y-1)
255   (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 );
256   interactionArray( i, j, k, l ).interact(particles, x, y);
257
258   // interaction with (x+1, y)
259   (x == m-1) ? (i=(x+1)%m, k=x) : (i=x, k=x+1);
260   interactionArray( i, y, k, y).interact(particles, x, y);
261
262   // interaction with (x+1, y+1)
263   (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 );
264   interactionArray( i, j, k, l ).interact(particles, x, y);
265
266 }
267
268 void Cell::force(CkVec<Particle> particles) {
269   forceCount++;
270   if( forceCount >= 9) {
271     
272     // Received all it's forces from the interactions.
273     
274     stepCount++;
275     forceCount = 0;
276     
277     /* FIX ME*/ // Methods to migrate atoms.
278       
279     #ifdef DEBUG
280       CkPrintf("STEP: %d DONE:( %d , %d )\n", stepCount, thisIndex.x, thisIndex.y);
281     #endif
282
283     if(stepCount >= finalStepCount) {
284       mainProxy.checkIn();
285     } else {
286       thisProxy( thisIndex.x, thisIndex.y ).start();
287     }
288   }
289     
290 }
291
292 // Function that receives a set of particles and updates the forces of them into the local set
293 void Cell::updateParticles(CkVec<Particle> updates) {
294         updateCount++;
295
296         //CHECKincomingParts.append(updates);
297
298         if(updateCount >= 8 ) {
299         
300                 
301                 
302                 updateCount = 0;
303         }
304
305 }
306
307 // Default constructor
308 Interaction::Interaction() {
309   cellCount = 0;
310
311   /* FIXME */
312     bufferedX = 0;
313     bufferedY = 0;
314 }
315
316 Interaction::Interaction(CkMigrateMessage *msg) { }
317   
318
319 // Function to receive vector of particles
320 void Interaction::interact(CkVec<Particle> particles, int x, int y ) {
321
322   if(cellCount == 0) {
323
324     // self interaction check
325     if( thisIndex.x == thisIndex.z && thisIndex.y == thisIndex.w ) {
326       CkPrintf("SELF: ( %d , %d )\n", thisIndex.x, thisIndex.y );
327       cellCount = 0;
328                         interact(particles,particles);
329       cellArray( x, y).force(particles);
330     } else {
331                          bufferedX = x;
332          bufferedY = y;
333                          bufferedParticles = particles;
334                 }
335
336   }
337
338   cellCount++;
339
340         // if both particle sets are received, compute interaction
341   if(cellCount >= 2) {
342
343     CkPrintf("PAIR:( %d , %d )  ( %d , %d ) \n", bufferedX, bufferedY, x, y );
344     cellCount = 0;
345
346                 interact(bufferedParticles,particles);
347
348     cellArray(bufferedX, bufferedY).force(bufferedParticles);
349     cellArray(x, y).force(particles);
350
351   }
352 }
353
354 // Function to compute all the interactions between pairs of particles in two sets
355 void Interaction::interact(CkVec<Particle> &first, CkVec<Particle> &second){
356         int i, j;
357         for(i = 0; i < first.length(); i++)
358                 for(j = 0; j < second.length(); j++)
359                         interact(first[i], second[j]);
360 }
361
362 // Function for computing interaction among two particles
363 // There is an extra test for interaction of identical particles, in which case there is no effect
364 void Interaction::interact(Particle &first, Particle &second){
365         float rx,ry,rz,r,fx,fy,fz,f;
366
367         // computing base values
368         rx = first.x - second.x;
369         ry = first.y - second.y;
370         r = sqrt(rx*rx + ry*ry);
371
372         if(r == 0.0 || r >= DEFAULT_RADIUS)
373                 return;
374
375         f = A / pow(r,12) - B / pow(r,6);
376         fx = f * rx / r;
377         fy = f * ry / r;
378
379         // updating particle properties
380         second.fx += fx;
381         second.fy += fy;
382         first.fx += -fx;
383         first.fy += -fy;
384
385 }
386
387 #include "main.def.h"