fixed bugs in cleanup, topleveltree construction and idepth parameter acceptance...
[charm.git] / examples / charm++ / barnes-charm / ParticleChunk.cpp
1 #include "barnes.h"
2
3 #include "barnes.decl.h"
4 ParticleChunk::ParticleChunk(int mleaf, int mcell){
5
6   usesAtSync = CmiTrue;
7
8   nstep = 0;
9   numMsgsToEachTp = new int[numTreePieces];
10   particlesToTps.resize(numTreePieces);
11   for(int i = 0; i < numTreePieces; i++){
12     numMsgsToEachTp[i] = 0;
13     particlesToTps[i].reserve(MAX_PARTICLES_PER_MSG);
14   }
15
16 };
17
18 ParticleChunk::~ParticleChunk(){
19   for(int i = 0; i < numTreePieces; i++){
20     particlesToTps[i].free();
21   }
22 }
23
24 void ParticleChunk::doAtSync(CkCallback &cb_){
25   mainCb = cb_;
26   AtSync();
27 }
28
29 void ParticleChunk::ResumeFromSync(){
30   contribute(0,0,CkReduction::concat,mainCb);
31 }
32
33 void ParticleChunk::pup(PUP::er &p){
34   CBase_ParticleChunk::pup(p);
35   if(p.isUnpacking()){
36     particlesToTps.resize(numTreePieces);
37     for(int i = 0; i < numTreePieces; i++){
38       particlesToTps[i].reserve(MAX_PARTICLES_PER_MSG);
39     }
40   }
41
42   p | mybodytab;
43   p | mynbody;
44   p | bodytab;
45   p | nstep;
46   PUParray(p, rmin, NDIM);
47   p | rsize;
48   p | mainCb;
49   p | numMsgsToEachTp;
50
51 }
52
53 void ParticleChunk::SlaveStart(CmiUInt8 bodyptrstart_, CmiUInt8 bodystart_, CkCallback &cb_){
54   unsigned int ProcessId;
55
56   /* Get unique ProcessId */
57   ProcessId = thisIndex;
58
59   bodyptr *bodyptrstart = (bodyptr *)bodyptrstart_;
60   bodyptr bodystart = (bodyptr)bodystart_;
61
62   /* POSSIBLE ENHANCEMENT:  Here is where one might pin processes to
63      processors to avoid migration */
64
65   /* initialize mybodytabs */
66   mybodytab = bodyptrstart + (maxmybody * ProcessId);
67   bodytab = bodystart;
68
69   /* note that every process has its own copy   */
70   /* of mybodytab, which was initialized to the */
71   /* beginning of the whole array by proc. 0    */
72   /* before create                              */
73   /* POSSIBLE ENHANCEMENT:  Here is where one might distribute the
74      data across physically distributed memories as desired. 
75
76      One way to do this is as follows:
77
78      int i;
79
80      if (ProcessId == 0) {
81      for (i=0;i<NPROC;i++) {
82      Place all addresses x such that 
83      &(Local[i]) <= x < &(Local[i])+
84      sizeof(struct local_memory) on node i
85      Place all addresses x such that 
86      &(Local[i].mybodytab) <= x < &(Local[i].mybodytab)+
87      maxmybody * sizeof(bodyptr) - 1 on node i
88      Place all addresses x such that 
89      &(Local[i].mycelltab) <= x < &(Local[i].mycelltab)+
90      maxmycell * sizeof(cellptr) - 1 on node i
91      Place all addresses x such that 
92      &(Local[i].myleaftab) <= x < &(Local[i].myleaftab)+
93      maxmyleaf * sizeof(leafptr) - 1 on node i
94      }
95      }
96
97      barrier(Global->Barstart,NPROC);
98
99 */
100
101   find_my_initial_bodies(bodytab, nbody, ProcessId);
102
103   contribute(0,0,CkReduction::concat,cb_);
104 }
105
106 /* 
107  * FIND_MY_INITIAL_BODIES: puts into mybodytab the initial list of bodies 
108  * assigned to the processor.  
109  */
110
111 void ParticleChunk::find_my_initial_bodies(bodyptr btab, int nbody, unsigned int ProcessId)
112 {
113   int Myindex;
114   int equalbodies;
115   int extra,offset,i;
116
117   int NPROC = numParticleChunks;
118
119   mynbody = nbody / NPROC;
120   extra = nbody % NPROC;
121   if (ProcessId < extra) {
122     mynbody++;    
123     offset = mynbody * ProcessId;
124   }
125   if (ProcessId >= extra) {
126     offset = (mynbody+1)*extra + (ProcessId - extra)*mynbody; 
127   }
128   for (i=0; i < mynbody; i++) {
129      mybodytab[i] = &(btab[offset+i]);
130   }
131 }
132
133 void ParticleChunk::acceptRoot(CmiUInt8 root_, real rminx, real rminy, real rminz, real rs, CkCallback &mainCb_){
134   mainCb = mainCb_;
135   G_root = (cellptr) root_;
136   rsize = rs;
137   rmin[0] = rminx;
138   rmin[1] = rminy;
139   rmin[2] = rminz;
140
141 #ifdef VERBOSE_CHUNKS
142   CkPrintf("[%d] acceptRoot 0x%x\n", thisIndex, G_root);
143 #endif
144   CkCallback cb(CkIndex_ParticleChunk::stepsystemPartII(0), thisProxy);
145   contribute(0,0,CkReduction::concat,cb);
146 }
147
148 void ParticleChunk::partition(CkCallback &cb){
149   int ProcessId = thisIndex;
150   HouseKeep(); 
151   real Cavg = (real) Cost(G_root) / (real)NPROC ;
152   workMin = (int) (Cavg * ProcessId);
153   workMax = (int) (Cavg * (ProcessId + 1) + (ProcessId == (NPROC - 1)));
154 #ifdef VERBOSE_CHUNKS
155   CkPrintf("[%d] cost(root) = %f, workMin: %f, workMax: %f\n", ProcessId, Cavg, workMin, workMax);
156 #endif
157
158   mynbody = 0;
159   find_my_bodies((nodeptr)G_root, 0, BRC_FUC, ProcessId);
160 #ifdef VERBOSE_CHUNKS
161   CkPrintf("[%d] mynbody: %d\n", thisIndex, mynbody);
162 #endif
163   contribute(0,0,CkReduction::concat,cb);
164 #ifdef MEMCHECK
165   CkPrintf("[%d] memcheck after partition\n", thisIndex);
166   CmiMemoryCheck();
167 #endif
168 }
169
170 void ParticleChunk::HouseKeep(){
171   myn2bcalc = mynbccalc = myselfint = 0;
172 }
173
174 void ParticleChunk::find_my_bodies(nodeptr mycell, int work, int direction, unsigned ProcessId){
175   int i;
176   leafptr l;
177   nodeptr qptr;
178
179   if (Type(mycell) == LEAF) {
180     l = (leafptr) mycell;      
181     for (i = 0; i < l->num_bodies; i++) {                                                            
182       if (work >= workMin-.1) {                                                  
183         if((mynbody) > maxmybody) {                                             
184           CkPrintf("[%d] find_my_bodies: needs more than %d bodies; increase fleaves (%f) . mynbody: %d\n",ProcessId, maxmybody, fleaves, mynbody); 
185           CkAbort("fleaves\n");
186         }    
187         mybodytab[mynbody++] = Bodyp(l)[i];
188       }                                                                                             
189       work += Cost(Bodyp(l)[i]);                                                                    
190       if (work >= workMax-.1) {                                                    
191         break;
192       }                                                                                             
193     }                                                                                                
194   }
195   else {
196     for(i = 0; (i < NSUB) && (work < (workMax - .1)); i++){                         
197       qptr = Subp(mycell)[Child_Sequence[direction][i]];                                            
198       if (qptr!=NULL) {                                                                             
199         if ((work+Cost(qptr)) >= (workMin -.1)) {                                 
200           find_my_bodies(qptr,work, Direction_Sequence[direction][i],                             
201               ProcessId);                                                              
202         }
203         work += Cost(qptr);                                                                        
204       } 
205     }
206   }  
207 }
208
209 void ParticleChunk::ComputeForces (CkCallback &cb)
210 {
211    bodyptr p,*pp;
212    vector acc1, dacc, dvel, vel1, dpos;
213    unsigned ProcessId = thisIndex;
214
215    for (pp = mybodytab; pp < mybodytab+mynbody; pp++) {  
216       p = *pp;
217       SETV(acc1, Acc(p));
218       Cost(p)=0;
219       //CkPrintf("forces for particle %d\n", p->num);
220       hackgrav(p,ProcessId);
221       myn2bcalc += myn2bterm; 
222       mynbccalc += mynbcterm;
223       if (skipself) {       /*   did we miss self-int?  */
224          myselfint++;        /*   count another goofup   */
225       }
226       if (nstep > 0) {
227          /*   use change in accel to make 2nd order correction to vel      */
228          SUBV(dacc, Acc(p), acc1);
229          MULVS(dvel, dacc, dthf);
230          ADDV(Vel(p), Vel(p), dvel);
231       }
232 #ifdef OUTPUT_ACC
233       p->n2b = myn2bterm;
234       p->nbc = mynbcterm;
235 #endif
236    }
237
238    contribute(0,0,CkReduction::concat,cb);
239 #ifdef MEMCHECK
240   CkPrintf("[%d] memcheck after calcforces\n", thisIndex);
241   CmiMemoryCheck();
242 #endif
243 }
244
245 void ParticleChunk::stepsystemPartII(CkReductionMsg *msg){
246
247     delete msg;
248
249     unsigned int ProcessId = thisIndex;
250
251     /* load bodies into tree   */
252     maketree(ProcessId);
253     flushParticles();
254     doneSendingParticles();
255 }
256
257 /*
258  * MAKETREE: initialize tree structure for hack force calculation.
259  */
260
261 void ParticleChunk::maketree(unsigned int ProcessId)
262 {
263    bodyptr p, *pp;
264
265    Current_Root = (nodeptr) G_root;
266    for (pp = mybodytab; pp < mybodytab+mynbody; pp++) {
267       p = *pp;
268       if (Mass(p) != 0.0) {
269          Current_Root = (nodeptr) loadtree(p, (cellptr) Current_Root, ProcessId);
270       }
271       else {
272          fprintf(stderr, "Process %d found body 0x%x to have zero mass\n",
273                  ProcessId, p); 
274       }
275    }
276 }
277
278 /*
279  * LOADTREE: descend tree and insert particle.
280  */
281
282 nodeptr
283 ParticleChunk::loadtree(bodyptr p, cellptr root, unsigned ProcessId)
284 {
285    int l, xq[NDIM], xp[NDIM], flag;
286    int i, j, root_level;
287    bool valid_root;
288    int kidIndex;
289    nodeptr *qptr, mynode;
290    cellptr c;
291    leafptr le;
292
293    intcoord(xp, Pos(p), rmin, rsize);
294    root = G_root;
295    mynode = (nodeptr) root;
296    kidIndex = subindex(xp, Level(mynode));
297
298    l = Level(mynode) >> 1;
299
300    int depth = log8floor(numTreePieces);
301    int lowestLevel = Level(mynode) >> depth;
302    int fact = NDIM;
303    int whichTp = 0;
304    int d = depth-1;
305
306    for(int level = Level(mynode); level > lowestLevel; level >>= 1){
307      kidIndex = subindex(xp, Level(mynode));
308      mynode = Subp(mynode)[kidIndex];
309      whichTp += kidIndex*(1<<(fact*(d)));
310      d--;     
311    }
312
313    int howMany = particlesToTps[whichTp].push_back_v(p); 
314    if(howMany == MAX_PARTICLES_PER_MSG-1){ // enough particles to send 
315      sendParticlesToTp(whichTp);
316    }
317 }
318
319 void ParticleChunk::flushParticles(){
320   for(int tp = 0; tp < numTreePieces; tp++){
321     sendParticlesToTp(tp); 
322   }
323 }
324
325 void ParticleChunk::sendParticlesToTp(int tp){
326   int len = particlesToTps[tp].length();
327   if(len > 0){
328 #ifdef VERBOSE_CHUNKS
329     CkPrintf("[%d] sending %d particles to piece %d\n", thisIndex, len, tp);
330 #endif
331     ParticleMsg *msg = new (len) ParticleMsg(); 
332     memcpy(msg->particles, particlesToTps[tp].getVec(), len*sizeof(bodyptr));
333     msg->num = len; 
334     particlesToTps[tp].length() = 0;
335     pieces[tp].recvParticles(msg);
336   }
337 }
338
339 void ParticleChunk::doneSendingParticles(){
340 }
341
342 void ParticleChunk::doneTreeBuild(){
343 #ifdef VERBOSE_CHUNKS
344   CkPrintf("[%d] all pieces have completed buildTree()\n", thisIndex);
345 #endif
346   mainCb.send();
347 }
348
349 void ParticleChunk::advance(CkCallback &cb_){
350   /* advance my bodies */
351
352   real minmax[NDIM*2];
353   int i;                                      
354   bodyptr p,*pp;                                                          
355   vector acc1, dacc, dvel, vel1, dpos;
356
357   mainCb = cb_;
358   SETVS(minmax,1E99);
359   SETVS((minmax+NDIM),-1E99);
360
361   for (pp = mybodytab; pp < mybodytab+mynbody; pp++) {
362     p = *pp;
363     MULVS(dvel, Acc(p), dthf);
364     ADDV(vel1, Vel(p), dvel);
365     MULVS(dpos, vel1, dtime);
366     ADDV(Pos(p), Pos(p), dpos);
367     ADDV(Vel(p), vel1, dvel);
368
369     for (i = 0; i < NDIM; i++) {
370       if (Pos(p)[i]<minmax[i]) {
371         minmax[i]=Pos(p)[i];
372       }
373       if (Pos(p)[i]>minmax[NDIM+i]) {
374         minmax[NDIM+i]=Pos(p)[i] ;
375       }
376     }
377   }
378
379   //CkPrintf("chunk %d minmax: (%f,%f,%f) (%f,%f,%f)\n", thisIndex, minmax[0], minmax[1], minmax[2], minmax[3], minmax[4], minmax[5]);
380   CkCallback cb(CkIndex_Main::recvGlobalSizes(NULL), mainChare);
381   contribute(NDIM*2*sizeof(real), minmax, minmax_RealVectorType, cb);
382 #ifdef MEMCHECK
383   CkPrintf("[%d] memcheck after advance\n", thisIndex);
384   CmiMemoryCheck();
385 #endif
386 }
387
388 void ParticleChunk::cleanup(){
389   nstep++;
390   contribute(0,0,CkReduction::concat,mainCb);
391 }
392
393 void ParticleChunk::outputAccelerations(CkCallback &cb_){
394   bodyptr *pp;
395   bodyptr p; 
396   int i;
397 #ifdef OUTPUT_ACC
398   for (i = 0, pp = mybodytab; pp < mybodytab+mynbody; pp++, i++) {  
399     p = *pp;
400     real *xp = Pos(p);
401     real *ap = Acc(p);
402     //CkPrintf("[%d] %d: pos: (%f,%f,%f), acc: (%f,%f,%f), nbc: %d, n2b: %d\n", thisIndex, p->num, xp[0], xp[1], xp[2], ap[0], ap[1], ap[2], p->nbc, p->n2b);
403     ckerr << p->num << ": pos: (" << xp[0] << "," << xp[1] << "," << xp[2] << "), acc: (" << ap[0] << "," << ap[1] << "," << ap[2] << "), nbc: " << p->nbc << ", n2b: " << p->n2b << endl;
404     //CkPrintf("%d: pos: (%f,%f,%f), acc: (%f,%f,%f), nbc: %d, n2b: %d\n", p->num, xp[0], xp[1], xp[2], ap[0], ap[1], ap[2], p->nbc, p->n2b);
405
406   }
407 #endif
408   //CkPrintf("[%d] bc: %d bb: %d\n", thisIndex, mynbccalc, myn2bcalc);
409   contribute(0,0,CkReduction::concat,cb_);
410 }
411