f6b8382f475fe59df8fd0a55dbd4359dcd4d9dfe
[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 }
165
166 void ParticleChunk::HouseKeep(){
167   myn2bcalc = mynbccalc = myselfint = 0;
168 }
169
170 void ParticleChunk::find_my_bodies(nodeptr mycell, int work, int direction, unsigned ProcessId){
171   int i;
172   leafptr l;
173   nodeptr qptr;
174
175   if (Type(mycell) == LEAF) {
176     l = (leafptr) mycell;      
177     for (i = 0; i < l->num_bodies; i++) {                                                            
178       if (work >= workMin-.1) {                                                  
179         if((mynbody) > maxmybody) {                                             
180           CkPrintf("[%d] find_my_bodies: needs more than %d bodies; increase fleaves (%f) . mynbody: %d\n",ProcessId, maxmybody, fleaves, mynbody); 
181           CkAbort("fleaves\n");
182         }    
183         mybodytab[mynbody++] = Bodyp(l)[i];
184       }                                                                                             
185       work += Cost(Bodyp(l)[i]);                                                                    
186       if (work >= workMax-.1) {                                                    
187         break;
188       }                                                                                             
189     }                                                                                                
190   }
191   else {
192     for(i = 0; (i < NSUB) && (work < (workMax - .1)); i++){                         
193       qptr = Subp(mycell)[Child_Sequence[direction][i]];                                            
194       if (qptr!=NULL) {                                                                             
195         if ((work+Cost(qptr)) >= (workMin -.1)) {                                 
196           find_my_bodies(qptr,work, Direction_Sequence[direction][i],                             
197               ProcessId);                                                              
198         }
199         work += Cost(qptr);                                                                        
200       } 
201     }
202   }  
203 }
204
205 void ParticleChunk::ComputeForces (CkCallback &cb)
206 {
207    bodyptr p,*pp;
208    vector acc1, dacc, dvel, vel1, dpos;
209    unsigned ProcessId = thisIndex;
210
211    for (pp = mybodytab; pp < mybodytab+mynbody; pp++) {  
212       p = *pp;
213       SETV(acc1, Acc(p));
214       Cost(p)=0;
215       //CkPrintf("forces for particle %d\n", p->num);
216       hackgrav(p,ProcessId);
217       myn2bcalc += myn2bterm; 
218       mynbccalc += mynbcterm;
219       if (skipself) {       /*   did we miss self-int?  */
220          myselfint++;        /*   count another goofup   */
221       }
222       if (nstep > 0) {
223          /*   use change in accel to make 2nd order correction to vel      */
224          SUBV(dacc, Acc(p), acc1);
225          MULVS(dvel, dacc, dthf);
226          ADDV(Vel(p), Vel(p), dvel);
227       }
228 #ifdef OUTPUT_ACC
229       p->n2b = myn2bterm;
230       p->nbc = mynbcterm;
231 #endif
232    }
233
234    contribute(0,0,CkReduction::concat,cb);
235 }
236
237 void ParticleChunk::stepsystemPartII(CkReductionMsg *msg){
238
239     delete msg;
240
241     unsigned int ProcessId = thisIndex;
242
243     /* load bodies into tree   */
244     maketree(ProcessId);
245     flushParticles();
246     doneSendingParticles();
247 }
248
249 /*
250  * MAKETREE: initialize tree structure for hack force calculation.
251  */
252
253 void ParticleChunk::maketree(unsigned int ProcessId)
254 {
255    bodyptr p, *pp;
256
257    Current_Root = (nodeptr) G_root;
258    for (pp = mybodytab; pp < mybodytab+mynbody; pp++) {
259       p = *pp;
260       if (Mass(p) != 0.0) {
261          Current_Root = (nodeptr) loadtree(p, (cellptr) Current_Root, ProcessId);
262       }
263       else {
264          fprintf(stderr, "Process %d found body 0x%x to have zero mass\n",
265                  ProcessId, p); 
266       }
267    }
268 }
269
270 /*
271  * LOADTREE: descend tree and insert particle.
272  */
273
274 nodeptr
275 ParticleChunk::loadtree(bodyptr p, cellptr root, unsigned ProcessId)
276 {
277    int l, xq[NDIM], xp[NDIM], flag;
278    int i, j, root_level;
279    bool valid_root;
280    int kidIndex;
281    nodeptr *qptr, mynode;
282    cellptr c;
283    leafptr le;
284
285    intcoord(xp, Pos(p), rmin, rsize);
286    root = G_root;
287    mynode = (nodeptr) root;
288    kidIndex = subindex(xp, Level(mynode));
289
290    l = Level(mynode) >> 1;
291
292    int depth = log8floor(numTreePieces);
293    int lowestLevel = Level(mynode) >> depth;
294    int fact = NSUB/2;
295    int whichTp = 0;
296    int d = depth;
297
298    for(int level = Level(mynode); level > lowestLevel; level >>= 1){
299      kidIndex = subindex(xp, Level(mynode));
300      mynode = Subp(mynode)[kidIndex];
301      whichTp += kidIndex*(1<<(fact*(d-1)));
302      d--;     
303    }
304
305    int howMany = particlesToTps[whichTp].push_back_v(p); 
306    if(howMany == MAX_PARTICLES_PER_MSG-1){ // enough particles to send 
307      sendParticlesToTp(whichTp);
308    }
309 }
310
311 void ParticleChunk::flushParticles(){
312   for(int tp = 0; tp < numTreePieces; tp++){
313     sendParticlesToTp(tp); 
314   }
315 }
316
317 void ParticleChunk::sendParticlesToTp(int tp){
318   int len = particlesToTps[tp].length();
319   if(len > 0){
320 #ifdef VERBOSE_CHUNKS
321     CkPrintf("[%d] sending %d particles to piece %d\n", thisIndex, len, tp);
322 #endif
323     ParticleMsg *msg = new (len) ParticleMsg(); 
324     memcpy(msg->particles, particlesToTps[tp].getVec(), len*sizeof(bodyptr));
325     msg->num = len; 
326     particlesToTps[tp].length() = 0;
327     pieces[tp].recvParticles(msg);
328   }
329 }
330
331 void ParticleChunk::doneSendingParticles(){
332 }
333
334 void ParticleChunk::doneTreeBuild(){
335 #ifdef VERBOSE_CHUNKS
336   CkPrintf("[%d] all pieces have completed buildTree()\n", thisIndex);
337 #endif
338   mainCb.send();
339 }
340
341 void ParticleChunk::advance(CkCallback &cb_){
342   /* advance my bodies */
343
344   real minmax[NDIM*2];
345   int i;                                      
346   bodyptr p,*pp;                                                          
347   vector acc1, dacc, dvel, vel1, dpos;
348
349   mainCb = cb_;
350   SETVS(minmax,1E99);
351   SETVS((minmax+NDIM),-1E99);
352
353   for (pp = mybodytab; pp < mybodytab+mynbody; pp++) {
354     p = *pp;
355     MULVS(dvel, Acc(p), dthf);
356     ADDV(vel1, Vel(p), dvel);
357     MULVS(dpos, vel1, dtime);
358     ADDV(Pos(p), Pos(p), dpos);
359     ADDV(Vel(p), vel1, dvel);
360
361     for (i = 0; i < NDIM; i++) {
362       if (Pos(p)[i]<minmax[i]) {
363         minmax[i]=Pos(p)[i];
364       }
365       if (Pos(p)[i]>minmax[NDIM+i]) {
366         minmax[NDIM+i]=Pos(p)[i] ;
367       }
368     }
369   }
370
371   CkCallback cb(CkIndex_Main::recvGlobalSizes(NULL), mainChare);
372   contribute(NDIM*2*sizeof(real), minmax, minmax_RealVectorType, cb);
373 }
374
375 void ParticleChunk::cleanup(){
376   nstep++;
377   contribute(0,0,CkReduction::concat,mainCb);
378 }
379
380 void ParticleChunk::outputAccelerations(CkCallback &cb_){
381   bodyptr *pp;
382   bodyptr p; 
383   int i;
384 #ifdef OUTPUT_ACC
385   for (i = 0, pp = mybodytab; pp < mybodytab+mynbody; pp++, i++) {  
386     p = *pp;
387     real *xp = Pos(p);
388     real *ap = Acc(p);
389     //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);
390     ckerr << p->num << ": pos: (" << xp[0] << "," << xp[1] << "," << xp[2] << "), acc: (" << ap[0] << "," << ap[1] << "," << ap[2] << "), nbc: " << p->nbc << ", n2b: " << p->n2b << endl;
391     //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);
392
393   }
394 #endif
395   //CkPrintf("[%d] bc: %d bb: %d\n", thisIndex, mynbccalc, myn2bcalc);
396   contribute(0,0,CkReduction::concat,cb_);
397 }
398