-
Notifications
You must be signed in to change notification settings - Fork 30
/
Copy pathTreePiece.cpp
6748 lines (5930 loc) · 220 KB
/
TreePiece.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/** @file TreePiece.cpp
*/
#include <cstdio>
#include <algorithm>
#include <fstream>
#include <assert.h>
#include <float.h>
// jetley
#include "limits.h"
#include "ParallelGravity.h"
#include "DataManager.h"
#include "Reductions.h"
// jetley
#include "MultistepLB.h"
#include "MultistepLB_SFC.h"
#include "MultistepLB_notopo.h"
#include "MultistepNodeLB_notopo.h"
#include "Orb3dLB.h"
#include "Orb3dLB_notopo.h"
#include "HierarchOrbLB.h"
// jetley - refactoring
//#include "codes.h"
#include "Opt.h"
#include "Compute.h"
#include "TreeWalk.h"
//#include "State.h"
#include "Space.h"
#include "gravity.h"
#include "smooth.h"
#include "PETreeMerger.h"
#include "IntraNodeLBManager.h"
#include "CkLoopAPI.h"
#include "formatted_string.h"
#if !CMK_LB_USER_DATA
#error "Please recompile charm with --enable-lbuserdata"
#endif
#ifdef PUSH_GRAVITY
#include "ckmulticast.h"
#endif
using namespace std;
using namespace SFC;
using namespace TreeStuff;
using namespace TypeHandling;
int TreeStuff::maxBucketSize;
#ifdef PUSH_GRAVITY
extern CkGroupID ckMulticastGrpId;
#endif
CkpvExtern(int, _lb_obj_index);
//forward declaration
string getColor(GenericTreeNode*);
const char *typeString(NodeType type);
/**
* @brief glassDamping applies a damping force to a particle's velocity.
*
* This can be useful for generating glasses. It is mean to model a damping
* term vdot = -damping * v
* @param v Particle velocity (v or vPred), a Vector3D
* @param dDelta Time step to apply damping over
* @param damping Inverse timescale of the damping
*/
inline void glassDamping(Vector3D<cosmoType> &v, double dDelta, double damping) {
#ifdef DAMPING
v *= exp(-dDelta * damping);
#endif
}
/*
* set periodic information in all the TreePieces
*/
void TreePiece::setPeriodic(int nRepsPar, // Number of replicas in
// each direction
Vector3D<cosmoType> fPeriodPar, // Size of periodic box
int bEwaldPar, // Use Ewald summation
double fEwCutPar, // Cutoff on real summation
double dEwhCutPar, // Cutoff on Fourier summation
int bPeriodPar, // Periodic boundaries
int bComovePar, // Comoving coordinates
double dRhoFacPar // Background density
)
{
nReplicas = nRepsPar;
fPeriod = fPeriodPar;
bEwald = bEwaldPar;
fEwCut = fEwCutPar;
dEwhCut = dEwhCutPar;
bPeriodic = bPeriodPar;
bComove = bComovePar;
dRhoFac = dRhoFacPar;
if(ewt == NULL) {
ewt = new EWT[nMaxEwhLoop];
}
}
// Scale velocities (needed to convert to canonical momenta for
// comoving coordinates.
void TreePiece::velScale(double dScale, const CkCallback& cb)
{
for(unsigned int i = 0; i < myNumParticles; ++i)
{
myParticles[i+1].velocity *= dScale;
if(TYPETest(&myParticles[i+1], TYPE_GAS)
#ifdef COLLISION
|| TYPETest(&myParticles[i+1], TYPE_DARK)
#endif
) {
myParticles[i+1].vPred() *= dScale;
}
}
contribute(cb);
}
/// After the bounding box has been found, we can assign keys to the particles
void TreePiece::assignKeys(CkReductionMsg* m) {
if(m->getSize() != sizeof(OrientedBox<float>)) {
ckerr << thisIndex << ": TreePiece: Fatal: Wrong size reduction message received!" << endl;
CkAssert(0);
callback.send(0);
delete m;
return;
}
boundingBox = *static_cast<OrientedBox<float> *>(m->getData());
delete m;
if(thisIndex == 0 && verbosity > 1)
ckout << "TreePiece: Bounding box originally: "
<< boundingBox << endl;
//give particles keys, using bounding box to scale
if((domainDecomposition!=ORB_dec)
&& (domainDecomposition!=ORB_space_dec)){
// get longest axis
Vector3D<float> bsize = boundingBox.size();
float max = (bsize.x > bsize.y) ? bsize.x : bsize.y;
max = (max > bsize.z) ? max : bsize.z;
//
// Make the bounding box cubical.
//
Vector3D<float> bcenter = boundingBox.center();
// The magic number below is approximately 2^(-19)
const float fEps = 1.0 + 1.91e-6; // slop to ensure keys fall
// between 0 and 1.
bsize = Vector3D<float>(fEps*0.5*max);
boundingBox = OrientedBox<float>(bcenter-bsize, bcenter+bsize);
if(thisIndex == 0 && verbosity > 1)
ckout << "TreePiece: Bounding box now: "
<< boundingBox << endl;
if(myNumParticles > 0) {
myParticles[0].key = firstPossibleKey;
myParticles[myNumParticles+1].key = lastPossibleKey;
for(unsigned int i = 0; i < myNumParticles; ++i) {
myParticles[i+1].key = generateKey(myParticles[i+1].position,
boundingBox);
}
sort(&myParticles[1], &myParticles[myNumParticles+1]);
}
}
#if COSMO_DEBUG > 1
auto file_name = make_formatted_string("tree.%d.%d.before",thisIndex,iterationNo);
char const* fout = file_name.c_str();
ofstream ofs(fout);
for (int i=1; i<=myNumParticles; ++i)
ofs << keyBits(myParticles[i].key,KeyBits) << " " << myParticles[i].position[0] << " "
<< myParticles[i].position[1] << " " << myParticles[i].position[2] << endl;
ofs.close();
#endif
if(verbosity >= 5)
cout << thisIndex << ": TreePiece: Assigned keys to all my particles" << endl;
contribute(callback);
}
/**************ORB Decomposition***************/
/// Three comparison routines used in sort and upper_bound
/// to order particles in each of three dimensions, respectively
bool comp_dim0(GravityParticle p1, GravityParticle p2) {
return p1.position[0] < p2.position[0];
}
bool comp_dim1(GravityParticle p1, GravityParticle p2) {
return p1.position[1] < p2.position[1];
}
bool comp_dim2(GravityParticle p1, GravityParticle p2) {
return p1.position[2] < p2.position[2];
}
///Initialize stuff before doing ORB decomposition
void TreePiece::initORBPieces(const CkCallback& cb){
OrientedBox<float> box = boundingBox;
orbBoundaries.clear();
orbBoundaries.push_back(myParticles+1);
orbBoundaries.push_back(myParticles+myNumParticles+1);
firstTime=true;
phase=0;
//Initialize function pointers
compFuncPtr[0]= &comp_dim0;
compFuncPtr[1]= &comp_dim1;
compFuncPtr[2]= &comp_dim2;
myBinCountsORB.clear();
myBinCountsORB.push_back(myNumParticles);
//Find out how many levels will be there before we go into the tree owned
//completely by the TreePiece
chunkRootLevel=0;
unsigned int tmp = numTreePieces;
while(tmp){
tmp >>= 1;
chunkRootLevel++;
}
chunkRootLevel--;
boxes = new OrientedBox<float>[chunkRootLevel+1];
splitDims = new char[chunkRootLevel+1];
boxes[0] = boundingBox;
contribute(sizeof(OrientedBox<float>), &box, boxReduction, cb);
}
/// Allocate memory for sorted particles.
/// @param cb callback after everything is sorted.
/// @param cback callback to perform now.
void TreePiece::initBeforeORBSend(unsigned int myCount,
unsigned int myCountGas,
unsigned int myCountStar,
const CkCallback& cb,
const CkCallback& cback){
callback = cb;
CkCallback nextCallback = cback;
if(numTreePieces == 1) {
myCount = myNumParticles;
myCountGas = myNumSPH;
myCountStar = myNumStar;
}
myExpectedCount = myCount;
myExpectedCountSPH = myCountGas;
myExpectedCountStar = myCountStar;
mySortedParticles.clear();
mySortedParticles.reserve(myExpectedCount);
mySortedParticlesSPH.clear();
mySortedParticlesSPH.reserve(myExpectedCountSPH);
mySortedParticlesStar.clear();
mySortedParticlesStar.reserve(myExpectedCountStar);
contribute(nextCallback);
}
void TreePiece::sendORBParticles(){
std::list<GravityParticle *>::iterator iter;
std::list<GravityParticle *>::iterator iter2;
int i=0;
int nCounts = myBinCountsORB.size()/3; // Gas counts are in second
// half and star counts are
// in third half
for(iter=orbBoundaries.begin();iter!=orbBoundaries.end();iter++,i++){
iter2=iter;
iter2++;
if(iter2==orbBoundaries.end())
break;
if(myBinCountsORB[i]>0) {
extraSPHData *pGasOut = NULL;
if(myBinCountsORB[nCounts+i] > 0) {
pGasOut = new extraSPHData[myBinCountsORB[nCounts+i]];
if (verbosity>=3)
CkPrintf("me:%d to:%d nPart :%ld, nGas:%d\n", thisIndex, i,
*iter2 - *iter, myBinCountsORB[nCounts+i]);
int iGasOut = 0;
for(GravityParticle *pPart = *iter; pPart < *iter2; pPart++) {
if(TYPETest(pPart, TYPE_GAS)) {
pGasOut[iGasOut] = *(extraSPHData *)pPart->extraData;
iGasOut++;
}
}
}
extraStarData *pStarOut = NULL;
if(myBinCountsORB[2*nCounts+i] > 0) {
pStarOut = new extraStarData[myBinCountsORB[2*nCounts+i]];
int iStarOut = 0;
for(GravityParticle *pPart = *iter; pPart < *iter2; pPart++) {
if(pPart->isStar()) {
pStarOut[iStarOut] = *(extraStarData *)pPart->extraData;
iStarOut++;
}
}
}
if(i==thisIndex){
acceptORBParticles(*iter,myBinCountsORB[i], pGasOut,
myBinCountsORB[nCounts+i], pStarOut,
myBinCountsORB[2*nCounts+i]);
}
else{
pieces[i].acceptORBParticles(*iter,myBinCountsORB[i], pGasOut,
myBinCountsORB[nCounts+i],
pStarOut,
myBinCountsORB[2*nCounts+i]);
}
if(pGasOut != NULL)
delete[] pGasOut;
if(pStarOut != NULL)
delete[] pStarOut;
}
}
if(myExpectedCount > (int) myNumParticles){
delete [] myParticles;
nStore = (int)((myExpectedCount + 2)*(1.0 + dExtraStore));
myParticles = new GravityParticle[nStore];
}
myNumParticles = myExpectedCount;
if(myExpectedCountSPH > (int) myNumSPH){
if(nStoreSPH > 0) delete [] mySPHParticles;
nStoreSPH = (int)(myExpectedCountSPH*(1.0 + dExtraStore));
mySPHParticles = new extraSPHData[nStoreSPH];
}
myNumSPH = myExpectedCountSPH;
if(myExpectedCountStar > (int) myNumStar){
delete [] myStarParticles;
allocateStars();
}
myNumStar = myExpectedCountStar;
if(myExpectedCount == 0) // No particles. Make sure transfer is
// complete
acceptORBParticles(NULL, 0, NULL, 0, NULL, 0);
}
/// Accept particles from other TreePieces once the sorting has finished
void TreePiece::acceptORBParticles(const GravityParticle* particles,
const int n,
const extraSPHData *pGas,
const int nGasIn,
const extraStarData *pStar,
const int nStarIn) {
copy(particles, particles + n, back_inserter(mySortedParticles));
copy(pGas, pGas + nGasIn, back_inserter(mySortedParticlesSPH));
copy(pStar, pStar + nStarIn, back_inserter(mySortedParticlesStar));
if(myExpectedCount == mySortedParticles.size()) {
//I've got all my particles
//Assigning keys to particles
for(int i=0;i<myExpectedCount;i++){
mySortedParticles[i].key = thisIndex;
}
copy(mySortedParticles.begin(), mySortedParticles.end(), &myParticles[1]);
copy(mySortedParticlesSPH.begin(), mySortedParticlesSPH.end(),
&mySPHParticles[0]);
copy(mySortedParticlesStar.begin(), mySortedParticlesStar.end(),
&myStarParticles[0]);
// assign gas and star data pointers
int iGas = 0;
int iStar = 0;
for(int i=0;i<myExpectedCount;i++){
if(myParticles[i+1].isGas()) {
myParticles[i+1].extraData
= (extraSPHData *)&mySPHParticles[iGas];
iGas++;
}
else if(myParticles[i+1].isStar()) {
myParticles[i+1].extraData
= (extraStarData *)&myStarParticles[iStar];
iStar++;
}
else
myParticles[i+1].extraData = NULL;
}
//signify completion with a reduction
if(verbosity>1)
ckout << thisIndex <<" contributing to accept particles"<<endl;
deleteTree();
contribute(callback);
}
}
/// @brief Determine my boundaries at the end of ORB decomposition
void TreePiece::finalizeBoundaries(ORBSplittersMsg *splittersMsg){
CkCallback cback = splittersMsg->cb;
std::list<GravityParticle *>::iterator iter;
std::list<GravityParticle *>::iterator iter2;
iter = orbBoundaries.begin();
iter2 = orbBoundaries.begin();
iter2++;
phase++;
int index = thisIndex >> (chunkRootLevel-phase+1);
Key lastBit;
lastBit = thisIndex >> (chunkRootLevel-phase);
lastBit = lastBit & 0x1;
boxes[phase] = boxes[phase-1];
if(lastBit){
boxes[phase].lesser_corner[splittersMsg->dim[index]] = splittersMsg->pos[index];
}
else{
boxes[phase].greater_corner[splittersMsg->dim[index]] = splittersMsg->pos[index];
}
splitDims[phase-1]=splittersMsg->dim[index];
for(int i=0;i<splittersMsg->length;i++){
int dimen=(int)splittersMsg->dim[i];
//Current location of the division is stored in a variable
//Evaluate the number of particles in each division
GravityParticle dummy;
Vector3D<double> divide(0.0,0.0,0.0);
divide[dimen] = splittersMsg->pos[i];
dummy.position = divide;
GravityParticle* divEnd = upper_bound(*iter,*iter2,dummy,compFuncPtr[dimen]);
orbBoundaries.insert(iter2,divEnd);
iter = iter2;
iter2++;
}
firstTime = true;
// First part is total particles, second part is gas counts, third
// part is star counts
myBinCountsORB.assign(6*splittersMsg->length,0);
copy(tempBinCounts.begin(),tempBinCounts.end(),myBinCountsORB.begin());
delete splittersMsg;
contribute(cback);
}
/// @brief Evaluate particle counts for ORB decomposition
/// @param m A message containing splitting dimensions and splits, and
/// the callback to contribute
/// Counts the particles of this treepiece on each side of the
/// splits. These counts are summed in a contribution to the
/// specified callback.
///
void TreePiece::evaluateParticleCounts(ORBSplittersMsg *splittersMsg)
{
CkCallback& cback = splittersMsg->cb;
// For each split, BinCounts has total lower, total higher.
// The second half of the array has the counts for gas particles.
// The third half of the array has the counts for star particles.
tempBinCounts.assign(6*splittersMsg->length,0);
std::list<GravityParticle *>::iterator iter;
std::list<GravityParticle *>::iterator iter2;
iter = orbBoundaries.begin();
iter2 = orbBoundaries.begin();
iter2++;
for(int i=0;i<splittersMsg->length;i++){
int dimen = (int)splittersMsg->dim[i];
if(firstTime){
sort(*iter,*iter2,compFuncPtr[dimen]);
}
//Evaluate the number of particles in each division
GravityParticle dummy;
GravityParticle* divStart = *iter;
Vector3D<double> divide(0.0,0.0,0.0);
divide[dimen] = splittersMsg->pos[i];
dummy.position = divide;
GravityParticle* divEnd = upper_bound(*iter,*iter2,dummy,compFuncPtr[dimen]);
tempBinCounts[2*i] = divEnd - divStart;
tempBinCounts[2*i + 1] = myBinCountsORB[i] - (divEnd - divStart);
int nGasLow = 0;
int nGasHigh = 0;
int nStarLow = 0;
int nStarHigh = 0;
for(GravityParticle *pPart = divStart; pPart < divEnd; pPart++) {
// Count gas
if(TYPETest(pPart, TYPE_GAS))
nGasLow++;
// Count stars
if(TYPETest(pPart, TYPE_STAR))
nStarLow++;
}
for(GravityParticle *pPart = divEnd; pPart < *iter2; pPart++) {
// Count gas
if(TYPETest(pPart, TYPE_GAS))
nGasHigh++;
// Count stars
if(TYPETest(pPart, TYPE_STAR))
nStarHigh++;
}
tempBinCounts[2*splittersMsg->length + 2*i] = nGasLow;
tempBinCounts[2*splittersMsg->length + 2*i + 1] = nGasHigh;
tempBinCounts[4*splittersMsg->length + 2*i] = nStarLow;
tempBinCounts[4*splittersMsg->length + 2*i + 1] = nStarHigh;
iter++; iter2++;
}
if(firstTime)
firstTime=false;
contribute(6*splittersMsg->length*sizeof(int), &(*tempBinCounts.begin()), CkReduction::sum_int, cback);
delete splittersMsg;
}
#ifdef REDUCTION_HELPER
void ReductionHelper::evaluateBoundaries(SFC::Key* keys, const int n, int skipEvery, const CkCallback& cb){
splitters.assign(keys, keys + n);
if(localTreePieces.presentTreePieces.size() == 0){
int numBins = skipEvery ? n - (n-1)/(skipEvery+1) - 1 : n - 1;
int64_t *dummy = new int64_t[numBins];
for(int i = 0; i < numBins; i++) dummy[i] = 0;
contribute(sizeof(int64_t)*numBins, dummy, CkReduction::sum_long, cb);
delete [] dummy;
return;
}
for(int i = 0; i < localTreePieces.presentTreePieces.size(); i++){
localTreePieces.presentTreePieces[i]->evaluateBoundaries(keys, n, skipEvery, cb);
}
}
void ReductionHelper::evaluateBoundaries(const CkBitVector &binsToSplit, const CkCallback& cb) {
std::vector<SFC::Key> newSplitters;
SFC::Key leftBound, rightBound;
newSplitters.reserve(splitters.size() * 4);
newSplitters.push_back(SFC::firstPossibleKey);
for (int i = 0; i < binsToSplit.Length(); i++) {
if (binsToSplit.Test(i) == true) {
leftBound = splitters[i];
rightBound = splitters[i + 1];
if (newSplitters.back() != (rightBound | 7L) ) {
if (newSplitters.back() != (leftBound | 7L)) {
newSplitters.push_back(leftBound | 7L);
}
newSplitters.push_back((leftBound / 4 * 3 + rightBound / 4) | 7L);
newSplitters.push_back((leftBound / 2 + rightBound / 2) | 7L);
newSplitters.push_back((leftBound / 4 + rightBound / 4 * 3) | 7L);
newSplitters.push_back(rightBound | 7L);
}
}
}
if (newSplitters.back() != lastPossibleKey) {
newSplitters.push_back(lastPossibleKey);
}
evaluateBoundaries(&newSplitters[0], newSplitters.size(), 0, cb);
}
#endif
/// Determine my part of the sorting histograms by counting the number
/// of my particles in each bin.
/// This routine assumes the particles in key order.
/// The parameter skipEvery means that every "skipEvery" bins counted, one
/// must be skipped. When skipEvery is set, the keys are in groups of
/// "skipEvery" size, and only splits within each group need to be
/// evaluated. Hence the counts between the end of one group, and the
/// start of the next group are not evaluated. This feature is used
/// by the Oct decomposition.
void TreePiece::evaluateBoundaries(SFC::Key* keys, const int n, int skipEvery, const CkCallback& cb){
#ifdef COSMO_EVENT
double startTimer = CmiWallTimer();
#endif
int numBins = skipEvery ? n - (n-1)/(skipEvery+1) - 1 : n - 1;
//this array will contain the number of particles I own in each bin
int64_t *myCounts;
#ifdef REDUCTION_HELPER
myCounts = new int64_t[numBins];
#else
//myBinCounts.assign(numBins, 0);
myBinCounts.resize(numBins);
myCounts = myBinCounts.getVec();
#endif
memset(myCounts, 0, numBins*sizeof(int64_t));
if (myNumParticles > 0) {
Key* endKeys = keys+n;
GravityParticle *binBegin = &myParticles[1];
GravityParticle *binEnd;
GravityParticle dummy;
GravityParticle *interpolatedBound;
GravityParticle *refinedLowerBound;
GravityParticle *refinedUpperBound;
//int binIter = 0;
//vector<int>::iterator binIter = myBinCounts.begin();
//vector<Key>::iterator keyIter = dm->boundaryKeys.begin();
Key* keyIter = lower_bound(keys, keys+n, binBegin->key);
int binIter = skipEvery ? (keyIter-keys) - (keyIter-keys-1) / (skipEvery+1) - 1: keyIter - keys - 1;
int skip = skipEvery ? skipEvery - (keyIter-keys-1) % (skipEvery+1) : -1;
if (binIter == -1) {
dummy.key = keys[0];
binBegin = upper_bound(binBegin, &myParticles[myNumParticles+1], dummy);
keyIter++;
binIter++;
skip = skipEvery ? skipEvery : -1;
}
for( ; keyIter != endKeys; ++keyIter) {
dummy.key = *keyIter;
if (domainDecomposition == SFC_dec ||
domainDecomposition == SFC_peano_dec ||
domainDecomposition == SFC_peano_dec_3D ||
domainDecomposition == SFC_peano_dec_2D) {
// try to guess a better upper bound
ptrdiff_t remainingParticles = &myParticles[myNumParticles + 1] - binBegin;
ptrdiff_t remainingBins = endKeys - keyIter;
ptrdiff_t interpolationInterval = remainingParticles / remainingBins;
ptrdiff_t scaledInterval =
(ptrdiff_t) ( (double) interpolationInterval * 1.5);
if (remainingParticles > scaledInterval) {
interpolatedBound = binBegin + scaledInterval;
}
else {
interpolatedBound = binBegin + interpolationInterval;
}
if (interpolatedBound->key <= dummy.key) {
refinedLowerBound = interpolatedBound;
refinedUpperBound = &myParticles[myNumParticles + 1];
}
else {
refinedLowerBound = binBegin;
refinedUpperBound = interpolatedBound;
}
/// find the last place I could put this splitter key in
/// my array of particles
binEnd = upper_bound(refinedLowerBound, refinedUpperBound, dummy);
}
else {
binEnd = upper_bound(binBegin, &myParticles[myNumParticles+1], dummy);
}
/// this tells me the number of particles between the
/// last two splitter keys
if (skip != 0) {
myCounts[binIter] = ((int64_t)(binEnd - binBegin));
++binIter;
--skip;
} else {
skip = skipEvery;
}
if(&myParticles[myNumParticles+1] <= binEnd) break;
binBegin = binEnd;
}
#ifdef COSMO_EVENTS
traceUserBracketEvent(boundaryEvaluationUE, startTimer, CmiWallTimer());
#endif
}
//send my bin counts back in a reduction
#ifdef REDUCTION_HELPER
reductionHelperProxy.ckLocalBranch()->reduceBinCounts(numBins, myCounts, cb);
delete[] myCounts;
#else
contribute(numBins * sizeof(int64_t), myCounts, CkReduction::sum_long, cb);
#endif
}
void TreePiece:: resetObjectLoad(const CkCallback& cb) {
// We assume that a rung 0 step has just been calculated before
// the checkpoint. Note that the load data here is from the
// previous rung 0 calculation. This should not be a bad approximation.
CkAssert(iPrevRungLB == 0); // check that the above comment is correct
setObjTime(savedPhaseLoad[0]);
contribute(cb);
}
void TreePiece::unshuffleParticlesWoDD(const CkCallback& callback) {
double tpLoad;
myShuffleMsg = NULL;
after_dd_callback = callback;
if (dm == NULL) {
dm = (DataManager*)CkLocalNodeBranch(dataManagerID);
}
tpLoad = getObjTime();
populateSavedPhaseData(iPrevRungLB, tpLoad, nPrevActiveParts);
//find my responsibility
myPlace = find(dm->responsibleIndex.begin(), dm->responsibleIndex.end(), thisIndex) - dm->responsibleIndex.begin();
if (myPlace == dm->responsibleIndex.size()) {
myPlace = -2;
}
setNumExpectedNeighborMsgs();
if (myNumParticles == 0) {
incomingParticlesSelf = true;
if (thisIndex == 0) {
CkCallback cbqd(CkIndex_TreePiece::shuffleAfterQD(), thisProxy);
CkStartQD(cbqd);
}
return;
}
sendParticlesDuringDD(true);
if (thisIndex == 0) {
CkCallback cbqd(CkIndex_TreePiece::shuffleAfterQD(), thisProxy);
CkStartQD(cbqd);
}
}
/*
* Accepts sorted particles from external TreePieces
*/
void TreePiece::acceptSortedParticlesFromOther(ParticleShuffleMsg *shuffleMsg) {
if(shuffleMsg == NULL) {
return;
}
// Copy the particles from shuffleMsg to the tmpShuffle array
myTmpShuffleParticle.insert(myTmpShuffleParticle.end(), shuffleMsg->particles,
shuffleMsg->particles + shuffleMsg->n);
// Copy the SPH particles from shuffleMsg to the tmpShuffle array
myTmpShuffleSphParticle.insert(myTmpShuffleSphParticle.end(),
shuffleMsg->pGas, shuffleMsg->pGas + shuffleMsg->nSPH);
// Copy the Star particles from shuffleMsg to the tmpShuffle array
myTmpShuffleStarParticle.insert(myTmpShuffleStarParticle.end(),
shuffleMsg->pStar, shuffleMsg->pStar + shuffleMsg->nStar);
incomingParticlesArrived += shuffleMsg->n;
savePhaseData(savedPhaseLoadTmp, savedPhaseParticleTmp, shuffleMsg->loads,
shuffleMsg->parts_per_phase, shuffleMsg->nloads);
delete shuffleMsg;
}
/*
* When reusing the splitters, we perform the migration of particles and then
* wait for QD. This method is called once the quiescence is detected.
*/
void TreePiece::shuffleAfterQD() {
// myShuffleMsg holds the particles that came from within this TreePiece
// (internal transfer).
if (myShuffleMsg != NULL) {
incomingParticlesArrived += myShuffleMsg->n;
savePhaseData(savedPhaseLoadTmp, savedPhaseParticleTmp, myShuffleMsg->loads,
myShuffleMsg->parts_per_phase, myShuffleMsg->nloads);
}
// This function is called after QD which means all the particles have arrived
// therefore set the particleCounts in the DataManager based on total external
// particles that arrived and the ones that arrived within the TreePiece.
if (myPlace != -2) {
dm->particleCounts[myPlace] = incomingParticlesArrived;
}
if (myPlace == -2 || dm->particleCounts[myPlace] == 0) {
// Special case where no particle is assigned to this TreePiece
incomingParticlesSelf = false;
incomingParticlesMsg.clear();
savedPhaseLoad.swap(savedPhaseLoadTmp);
savedPhaseParticle.swap(savedPhaseParticleTmp);
savedPhaseLoadTmp.clear();
savedPhaseParticleTmp.clear();
if(verbosity>1) ckout << thisIndex <<" no particles assigned"<<endl;
deleteTree();
int isTPEmpty = 0;
if (myPlace != -2) {
isTPEmpty = 1;
}
contribute(sizeof(int), &isTPEmpty, CkReduction::logical_or, after_dd_callback);
if (myShuffleMsg != NULL) {
delete myShuffleMsg;
}
return;
}
nStore = (int)((dm->particleCounts[myPlace] + 2)*(1.0 + dExtraStore));
myParticles = new GravityParticle[nStore];
myNumParticles = dm->particleCounts[myPlace];
incomingParticlesArrived = 0;
incomingParticlesSelf = false;
savedPhaseLoad.swap(savedPhaseLoadTmp);
savedPhaseParticle.swap(savedPhaseParticleTmp);
savedPhaseLoadTmp.clear();
savedPhaseParticleTmp.clear();
// Merge all the particles which includes the ones received from within and
// from outside
mergeAllParticlesAndSaveCentroid();
//signify completion with a reduction
if(verbosity>1) ckout << thisIndex <<" contributing to accept particles"
<<endl;
if (myShuffleMsg != NULL) {
delete myShuffleMsg;
}
deleteTree();
int isTPEmpty = 0;
contribute(sizeof(int), &isTPEmpty, CkReduction::logical_or, after_dd_callback);
}
/*
* Merge the particles in sorted order. The particles include the ones that are
* within this TreePiece as well as received from outside.
* Since we are iterating over the particles, we might as well calculate the
* centroid of the TreePiece.
*/
void TreePiece::mergeAllParticlesAndSaveCentroid() {
int nSPH = 0;
int nStar = 0;
// myShuffleLoc keeps a count of number of external particles received.
nSPH = myTmpShuffleSphParticle.size();
nStar = myTmpShuffleStarParticle.size();
if (myShuffleMsg != NULL) {
nSPH += myShuffleMsg->nSPH;
nStar += myShuffleMsg->nStar;
}
myNumSPH = nSPH;
nStoreSPH = (int)(myNumSPH*(1.0 + dExtraStore));
if(nStoreSPH > 0) mySPHParticles = new extraSPHData[nStoreSPH];
else mySPHParticles = NULL;
myNumStar = nStar;
allocateStars();
nSPH = 0;
nStar = 0;
// myTmpShuffle__Particle holds the particles that came from external
// TreePieces
// Copy the gas and star data from the tmpshufflemsg array
memcpy(&mySPHParticles[nSPH], &myTmpShuffleSphParticle[0],
myTmpShuffleSphParticle.size()*sizeof(extraSPHData));
nSPH += myTmpShuffleSphParticle.size();
memcpy(&myStarParticles[nStar], &myTmpShuffleStarParticle[0],
myTmpShuffleStarParticle.size()*sizeof(extraStarData));
nStar += myTmpShuffleStarParticle.size();
int iGas = 0;
int iStar = 0;
// Set the location of the star data in the particle data for the ones that
// came in from an external TreePiece.
for (int i = 0; i < myTmpShuffleParticle.size(); i++) {
if (myTmpShuffleParticle[i].isGas()) {
myTmpShuffleParticle[i].extraData =
(extraSPHData *) &mySPHParticles[iGas];
iGas++;
} else if (myTmpShuffleParticle[i].isStar()) {
myTmpShuffleParticle[i].extraData =
(extraStarData *) &myStarParticles[iStar];
iStar++;
} else {
myTmpShuffleParticle[i].extraData = NULL;
}
}
// Now copy the SPH and Star for particles that moved within the TreePiece
if (myShuffleMsg != NULL) {
memcpy(&mySPHParticles[nSPH], myShuffleMsg->pGas,
(myShuffleMsg->nSPH)*sizeof(extraSPHData));
nSPH += (myShuffleMsg->nSPH);
memcpy(&myStarParticles[nStar], myShuffleMsg->pStar,
(myShuffleMsg->nStar)*sizeof(extraStarData));
nStar += (myShuffleMsg->nStar);
}
// sort is [first, last)
// Note that the particles that were received from outside were just
// appended to the myTmpShuffleParticle. Though within each shuffleMsg they
// will be sorted, when they are appended the myTmpShuffleParticle need not
// sorted.
// Sort the items that came from outside. Since we expect them to be a small
// number this sort is used
sort(myTmpShuffleParticle.begin(), myTmpShuffleParticle.end());
int left, right, tmp;
left = right = 0;
tmp = 1;
int leftend = 0;
if (myShuffleMsg != NULL) {
leftend = myShuffleMsg->n;
}
int rightend = myTmpShuffleParticle.size();
Vector3D<double> vCenter(0.0, 0.0, 0.0);
// merge sort. myShuffleMsg contains sorted particles and so does
// myTmpShuffleParticles. myTmpShuffleParticles contain particles that were
// transferred from outside this TreePiece. As they come in, they are merged
// in sorted order.
// For the particles that came from outside, the extraData field is already
// set above.
// when merging the pointer to the star data has to be set for the particles
// that moved within.
while (left < leftend && right < rightend) {
if (myShuffleMsg->particles[left] < myTmpShuffleParticle[right]) {
myParticles[tmp] = myShuffleMsg->particles[left++];
if (myParticles[tmp].isGas()) {
myParticles[tmp].extraData = (extraSPHData *) &mySPHParticles[iGas++];
} else if (myParticles[tmp].isStar()) {
myParticles[tmp].extraData = (extraStarData *) &myStarParticles[iStar++];
} else {
myParticles[tmp].extraData = NULL;
}
} else {
myParticles[tmp] = myTmpShuffleParticle[right++];
}
vCenter += myParticles[tmp].position;
tmp++;
}
while (left < leftend) {
myParticles[tmp] = myShuffleMsg->particles[left++];
if (myParticles[tmp].isGas()) {
myParticles[tmp].extraData = (extraSPHData *) &mySPHParticles[iGas++];
} else if (myParticles[tmp].isStar()) {
myParticles[tmp].extraData = (extraStarData *) &myStarParticles[iStar++];
} else {
myParticles[tmp].extraData = NULL;
}
vCenter += myParticles[tmp].position;
tmp++;
}
while (right < rightend) {
myParticles[tmp] = myTmpShuffleParticle[right++];
vCenter += myParticles[tmp].position;
tmp++;
}
// Clear all the tmp datastructures which were holding the migrated particles
myTmpShuffleParticle.clear();
myTmpShuffleSphParticle.clear();
myTmpShuffleStarParticle.clear();
// Save the centroid
savedCentroid = vCenter/(double)myNumParticles;
}
void TreePiece::setNumExpectedNeighborMsgs() {
nbor_msgs_count_ = 2;
// This TreePiece is out of the responsible index range
if (myPlace == -2) {
nbor_msgs_count_ = 0;
}
// This TreePiece is the first one so will get only from my right neighbor
if (myPlace == 0) {
nbor_msgs_count_--;
}
// This TreePiece is the last one so will get only from my left neighbor
if (myPlace == (dm->responsibleIndex.size()-1)) {
nbor_msgs_count_--;
}