testFmmAlgorithmProc.cpp 44.3 KB
Newer Older
1
2
// /!\ Please, you must read the license at the bottom of this page

3
#include "../Src/Utils/FMpi.hpp"
4
#include "../Src/Utils/FTic.hpp"
5

6
#include "../Src/Containers/FOctree.hpp"
7
#include "../Src/Containers/FVector.hpp"
8
#include "../Src/Utils/FParameters.hpp"
9
#include "../Src/Utils/FGlobal.hpp"
10

11
#include "../Src/Components/FSimpleLeaf.hpp"
12

13
#include "../Src/Utils/F3DPosition.hpp"
14
#include "../Src/Utils/FAbstractSendable.hpp"
15
#include "../Src/Utils/FQuickSort.hpp"
16

17
18
19
20
#include "../Src/Components/FFmaParticle.hpp"
#include "../Src/Components/FTestParticle.hpp"
#include "../Src/Components/FTestCell.hpp"
#include "../Src/Components/FTestKernels.hpp"
21
#include "../Src/Extensions/FExtendPhysicalValue.hpp"
22

23
#include "../Src/Core/FFmmAlgorithmThreadProc.hpp"
berenger-bramas's avatar
berenger-bramas committed
24
25
#include "../Src/Core/FFmmAlgorithmThread.hpp"

26
27
#include "../Src/Files/FFmaBinLoader.hpp"
#include "../Src/Files/FProcFmaLoader.hpp"
28

29
#include "../Src/Components/FBasicKernels.hpp"
30

berenger-bramas's avatar
berenger-bramas committed
31
32
33
34
35
36
#include <iostream>

#include <stdio.h>
#include <stdlib.h>


37
// Compile by : g++ testFmmAlgorithmProc.cpp ../Src/Utils/FDebug.cpp ../Src/Utils/FTrace.cpp -lgomp -fopenmp -O2 -o testFmmAlgorithmProc.exe
38
39

/** This program show an example of use of
berenger-bramas's avatar
berenger-bramas committed
40
  * the fmm threaded + mpi algo
41
42
43
44
  * it also check that each particles is impacted each other particles
  */


45
/** Fmb class has to extend {FExtendForces,FExtendPotential,FExtendPhysicalValue}
46
  * Because we use fma loader it needs {FExtendPhysicalValue}
47
48
49
50
51
  */
class TestParticle : public FTestParticle, public FExtendPhysicalValue {
public:
};

berenger-bramas's avatar
berenger-bramas committed
52
class FTestCellPar : public FTestCell, public FAbstractSendable{
53
public :
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
    int bytesToSendUp() const{
        return sizeof(long);
    }
    int writeUp(void* const buffer, const int) const {
        const long tmpUp = getDataUp();
        memcpy(buffer,&tmpUp,bytesToSendUp());
        return bytesToSendUp();
    }
    int bytesToReceiveUp() const{
        return sizeof(long);
    }
    int readUp(void* const buffer, const int) {
        long tmpUp;
        memcpy(&tmpUp,buffer,bytesToSendUp());
        setDataUp(tmpUp);
        return bytesToReceiveUp();
    }

    int bytesToSendDown() const{
        return sizeof(long);
    }
    int writeDown(void* const buffer, const int) const {
        const long tmpDown = getDataDown();
        memcpy(buffer,&tmpDown,bytesToSendDown());
        return bytesToSendDown();
    }
    int bytesToReceiveDown() const{
        return sizeof(long);
    }
    int readDown(void* const buffer, const int) {
        long tmpDown;
        memcpy(&tmpDown,buffer,bytesToSendDown());
        setDataDown(tmpDown + getDataDown());
        return bytesToReceiveDown();
88
89
90
91
    }
};


92
93
94
95
96
97
98
/////////////////////////////////////////////////////////////////////////////
// Test function
/////////////////////////////////////////////////////////////////////////////

/** This function test the octree to be sure that the fmm algorithm
  * has worked completly.
  */
berenger-bramas's avatar
berenger-bramas committed
99
100
101
102
template<class OctreeClass, class ContainerClass, class FmmClassProc>
void ValidateFMMAlgoProc(OctreeClass* const badTree,
                         OctreeClass* const valideTree,
                         FmmClassProc*const fmm){
103
    const int OctreeHeight = badTree->getHeight();
104
105
    std::cout << "Check Result\n";
    {
berenger-bramas's avatar
berenger-bramas committed
106
        typename OctreeClass::Iterator octreeIterator(badTree);
107
108
        octreeIterator.gotoBottomLeft();

berenger-bramas's avatar
berenger-bramas committed
109
        typename OctreeClass::Iterator octreeIteratorValide(valideTree);
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
        octreeIteratorValide.gotoBottomLeft();

        for(int level = OctreeHeight - 1 ; level > 0 ; --level){
            int NbLeafs = 0;
            do{
                ++NbLeafs;
            } while(octreeIterator.moveRight());
            octreeIterator.gotoLeft();

            const int startIdx = fmm->getLeft(NbLeafs);
            const int endIdx = fmm->getRight(NbLeafs);
            // Check that each particle has been summed with all other

            for(int idx = 0 ; idx < startIdx ; ++idx){
                octreeIterator.moveRight();
                octreeIteratorValide.moveRight();
            }

            for(int idx = startIdx ; idx < endIdx ; ++idx){
                if(octreeIterator.getCurrentGlobalIndex() != octreeIteratorValide.getCurrentGlobalIndex()){
                    std::cout << "Error index are not equal!" << std::endl;
                }
                else{
                    if(octreeIterator.getCurrentCell()->getDataUp() != octreeIteratorValide.getCurrentCell()->getDataUp()){
                        std::cout << "M2M error at level " << level << " up bad " << octreeIterator.getCurrentCell()->getDataUp()
                                << " good " << octreeIteratorValide.getCurrentCell()->getDataUp() << " idx " << idx << std::endl;
                    }
                    if(octreeIterator.getCurrentCell()->getDataDown() != octreeIteratorValide.getCurrentCell()->getDataDown()){
                        std::cout << "L2L error at level " << level << " down bad " << octreeIterator.getCurrentCell()->getDataDown()
                                << " good " << octreeIteratorValide.getCurrentCell()->getDataDown() << " idx " << idx << std::endl;
                    }
                }

                octreeIterator.moveRight();
                octreeIteratorValide.moveRight();
            }

            octreeIterator.moveUp();
            octreeIterator.gotoLeft();

            octreeIteratorValide.moveUp();
            octreeIteratorValide.gotoLeft();
        }
    }
    {
        int NbPart = 0;
        int NbLeafs = 0;
        { // Check that each particle has been summed with all other
berenger-bramas's avatar
berenger-bramas committed
158
            typename OctreeClass::Iterator octreeIterator(badTree);
159
160
161
162
163
164
165
166
167
168
169
            octreeIterator.gotoBottomLeft();
            do{
                NbPart += octreeIterator.getCurrentListSrc()->getSize();
                ++NbLeafs;
            } while(octreeIterator.moveRight());
            std::cout << "There is " << NbPart << " particles on " << NbLeafs << " Leafs" << std::endl;
        }
        {
            const int startIdx = fmm->getLeft(NbLeafs);
            const int endIdx = fmm->getRight(NbLeafs);
            // Check that each particle has been summed with all other
berenger-bramas's avatar
berenger-bramas committed
170
            typename OctreeClass::Iterator octreeIterator(badTree);
171
172
173
174
175
176
177
            octreeIterator.gotoBottomLeft();

            for(int idx = 0 ; idx < startIdx ; ++idx){
                octreeIterator.moveRight();
            }

            for(int idx = startIdx ; idx < endIdx ; ++idx){
berenger-bramas's avatar
berenger-bramas committed
178
                typename ContainerClass::BasicIterator iter(*octreeIterator.getCurrentListTargets());
179
180
181

                const bool isUsingTsm = (octreeIterator.getCurrentListTargets() != octreeIterator.getCurrentListSrc());

berenger-bramas's avatar
berenger-bramas committed
182
                while( iter.hasNotFinished() ){
183
184
                    // If a particles has been impacted by less than NbPart - 1 (the current particle)
                    // there is a problem
berenger-bramas's avatar
berenger-bramas committed
185
186
187
                    if( (!isUsingTsm && iter.data().getDataDown() != NbPart - 1) ||
                        (isUsingTsm && iter.data().getDataDown() != NbPart) ){
                        std::cout << "Problem L2P + P2P, value on particle is : " << iter.data().getDataDown() << "\n";
188
                    }
berenger-bramas's avatar
berenger-bramas committed
189
                    iter.gotoNext();
190
191
192
193
194
195
196
197
198
199
200
201
                }
                octreeIterator.moveRight();
            }
        }
    }

    std::cout << "Done\n";
}

/** To print an octree
  * used to debug and understand how the values were passed
  */
berenger-bramas's avatar
berenger-bramas committed
202
203
204
template<class OctreeClass>
void print(OctreeClass* const valideTree){
    typename OctreeClass::Iterator octreeIterator(valideTree);
205
    for(int idxLevel = valideTree->getHeight() - 1 ; idxLevel > 1 ; --idxLevel ){
206
207
208
209
210
211
212
213
214
        do{
            std::cout << "[" << octreeIterator.getCurrentGlobalIndex() << "] up:" << octreeIterator.getCurrentCell()->getDataUp() << " down:" << octreeIterator.getCurrentCell()->getDataDown() << "\t";
        } while(octreeIterator.moveRight());
        std::cout << "\n";
        octreeIterator.gotoLeft();
        octreeIterator.moveDown();
    }
}

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
class MyFackParticle : public FBasicParticle {
    int myIndex;
public:
    MyFackParticle() : myIndex(0) {
    }
    void setIndex(const int inIndex){
        this->myIndex = inIndex;
    }
    int getIndex() const{
        return this->myIndex;
    }
};

// My leaf store the indexes of the particles it receives
// in a vector
class MyLeaf : public FAbstractLeaf<MyFackParticle, FVector<int> > {
    FVector<int> indexes;
public:
    void push(const MyFackParticle& particle){
        indexes.push( particle.getIndex() );
    }
    FVector<int>* getSrc(){
        return &indexes;
    }
    FVector<int>* getTargets(){
        return &indexes;
    }
};

struct ParticlesGroup {
    int number;
    int positionInArray;
    MortonIndex index;
    ParticlesGroup(const int inNumber = 0 , const int inPositionInArray = 0, const MortonIndex inIndex = 0)
        : number(inNumber), positionInArray(inPositionInArray), index(inIndex) {
    }
};

253

254
255
256
257
258
259
260
typedef TestParticle               ParticleClass;
typedef FTestCellPar               CellClass;
typedef FVector<ParticleClass>     ContainerClass;

typedef FSimpleLeaf<ParticleClass, ContainerClass >                     LeafClass;
typedef FOctree<ParticleClass, CellClass, ContainerClass , LeafClass >  OctreeClass;
typedef FTestKernels<ParticleClass, CellClass, ContainerClass >         KernelClass;
berenger-bramas's avatar
berenger-bramas committed
261

262
263
typedef FFmmAlgorithmThread<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass >     FmmClass;
typedef FFmmAlgorithmThreadProc<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass >     FmmClassProc;
berenger-bramas's avatar
berenger-bramas committed
264

265
266
267
268

struct IndexedParticle{
    MortonIndex index;
    ParticleClass particle;
269
270
271
272

    operator MortonIndex(){
        return this->index;
    }
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
};

long getTreeCoordinate(const FReal inRelativePosition, const FReal boxWidthAtLeafLevel) {
        const FReal indexFReal = inRelativePosition / boxWidthAtLeafLevel;
        const long index = FMath::dfloor(indexFReal);
        if( index && FMath::LookEqual(inRelativePosition, boxWidthAtLeafLevel * index ) ){
                return index - 1;
        }
        return index;
}



// Simply create particles and try the kernels
int main(int argc, char ** argv){
288
289
290
291
    ///////////////////////What we do/////////////////////////////
    std::cout << ">> This executable has to be used to test the FMM algorithm.\n";
    //////////////////////////////////////////////////////////////

292
293
    FMpi app( argc, argv);

294
295
    const int NbLevels = FParameters::getValue(argc,argv,"-h", 9);
    const int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 3);
296
    char defaultFilename[] = "testLoaderFMA.fma"; //../../Data/ "testLoaderFMA.fma" "testFMAlgorithm.fma" Sphere.fma
297
    char* filename;
298
299
    FTic counter;

300
301
302
303
304
305
306
307
308
    if(argc == 1){
        std::cout << "You have to give a .fma file in argument.\n";
        std::cout << "The program will try a default file : " << defaultFilename << "\n";
        filename = defaultFilename;
    }
    else{
        filename = argv[1];
        std::cout << "Opening : " << filename << "\n";
    }
309

310
    FProcFmaLoader<ParticleClass> loader(filename,app);
311
    if(!loader.isOpen()){
312
313
314
        std::cout << "Loader Error, " << filename << " is missing\n";
        return 1;
    }
315

316
317
318
319
320
321
322
323
324
325
326

    OctreeClass realTree(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());

    {
        OctreeClass treeInterval(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
        int myNbParticlesCounter = 0;

        //////////////////////////////////////////////////////////////////////////////////
        // We sort our particles
        //////////////////////////////////////////////////////////////////////////////////
        {
327
328
            std::cout << "Inserting & Sorting my particles ..." << std::endl;
            counter.tic();
329
330

            //////////////////////////////////////////////////////////////////////////////////
331
            // Quick sort
332
            //////////////////////////////////////////////////////////////////////////////////
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
            if(true){
                IndexedParticle* outputArray = 0;
                long outputSize = 0;
                {
                    // create particles
                    IndexedParticle*const realParticlesIndexed = reinterpret_cast<IndexedParticle*>(new char[loader.getNumberOfParticles() * sizeof(IndexedParticle)]);
                    F3DPosition boxCorner(loader.getCenterOfBox() - (loader.getBoxWidth()/2));
                    FTreeCoordinate host;
                    const FReal boxWidthAtLeafLevel = loader.getBoxWidth() / (1 << (NbLevels - 1) );
                    for(long idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
                        loader.fillParticle(realParticlesIndexed[idxPart].particle);
                        host.setX( getTreeCoordinate( realParticlesIndexed[idxPart].particle.getPosition().getX() - boxCorner.getX(), boxWidthAtLeafLevel ));
                        host.setY( getTreeCoordinate( realParticlesIndexed[idxPart].particle.getPosition().getY() - boxCorner.getY(), boxWidthAtLeafLevel ));
                        host.setZ( getTreeCoordinate( realParticlesIndexed[idxPart].particle.getPosition().getZ() - boxCorner.getZ(), boxWidthAtLeafLevel ));
                        realParticlesIndexed[idxPart].index = host.getMortonIndex(NbLevels - 1);
                    }
349

350
                    std::cout << "Start to sort " << loader.getNumberOfParticles() << " particles..." << std::endl;
351

352
353
354
                    // sort particles
                    FQuickSort::QsMpi<IndexedParticle,MortonIndex>(realParticlesIndexed, loader.getNumberOfParticles(),outputArray,outputSize);
                    delete [] reinterpret_cast<char*>(realParticlesIndexed);
355

356
                    std::cout << "Sorted "<< outputSize <<  " particles..." << std::endl;
357
                }
358
359
360
361
362
363
364
365
366
367
368
369
370
371
                // be sure there is no splited leaves
                {
                    MortonIndex otherFirstIndex = -1;
                    {
                        FMpi::Request req[2];
                        int reqiter = 0;
                        if( 0 < app.processId() && outputSize){
                            std::cout << "I send to left that my first index is " << outputArray[0].index << std::endl;
                            MPI_Isend( &outputArray[0].index, 1, MPI_LONG_LONG, app.processId() - 1, 0, MPI_COMM_WORLD, &req[reqiter++]);
                        }
                        if( app.processId() != app.processCount() - 1){
                            MPI_Irecv(&otherFirstIndex, 1, MPI_LONG_LONG, app.processId() + 1, 0, MPI_COMM_WORLD, &req[reqiter++]);
                            std::cout << "I receive index from right " << otherFirstIndex << std::endl;
                        }
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
                        MPI_Waitall(reqiter,req,0);
                        if( 0 < app.processId() && !outputSize){
                            std::cout << "I send to left the first index is from right " << otherFirstIndex << std::endl;
                            MPI_Send( &otherFirstIndex, 1, MPI_LONG_LONG, app.processId() - 1, 0, MPI_COMM_WORLD);
                        }
                    }
                    // at this point every one know the first index of his right neighbors
                    const bool needToRecvBeforeSend = (app.processId() != 0 && ((outputSize && outputArray[0].index == otherFirstIndex) || !outputSize));
                    if( needToRecvBeforeSend || (app.processId() == app.processCount() - 1) ){
                        long sendByOther = 0;
                        MPI_Recv(&sendByOther, 1, MPI_LONG, app.processId() - 1, 0, MPI_COMM_WORLD,0);
                        std::cout << "I will receive " << sendByOther << std::endl;
                        if(sendByOther){
                            const IndexedParticle* const reallocOutputArray = outputArray;
                            const long reallocOutputSize = outputSize;

                            outputSize += sendByOther;
                            outputArray = new IndexedParticle[outputSize];
                            memcpy(&outputArray[sendByOther], reallocOutputArray, reallocOutputSize * sizeof(IndexedParticle));
                            delete[] reallocOutputArray;

                            MPI_Recv(outputArray, sizeof(IndexedParticle) * sendByOther, MPI_BYTE, app.processId() - 1, 0, MPI_COMM_WORLD, 0);
                        }
                    }
                    if(app.processId() != app.processCount() - 1){
                        long idxPart = outputSize - 1 ;
                        while(idxPart >= 0 && outputArray[idxPart].index == otherFirstIndex){
                            --idxPart;
                        }
                        long toSend = outputSize - 1 - idxPart;
                        MPI_Send( &toSend, 1, MPI_LONG, app.processId() + 1, 0, MPI_COMM_WORLD);
                        std::cout << "send to other " << toSend << std::endl;
                        if(toSend){
                            MPI_Send( &outputArray[idxPart + 1], toSend * sizeof(IndexedParticle), MPI_BYTE, app.processId() + 1, 0, MPI_COMM_WORLD);
                        }
                        if( app.processId() != 0 && !needToRecvBeforeSend ){
                            long sendByOther = 0;
                            MPI_Recv(&sendByOther, 1, MPI_LONG, app.processId() - 1, 0, MPI_COMM_WORLD, 0);
                            std::cout << "I will receive " << sendByOther << std::endl;
                            if(sendByOther){
                                const IndexedParticle* const reallocOutputArray = outputArray;
                                const long reallocOutputSize = outputSize;

                                outputSize += sendByOther;
                                outputArray = new IndexedParticle[outputSize];
                                memcpy(&outputArray[sendByOther], reallocOutputArray, reallocOutputSize * sizeof(IndexedParticle));
                                delete[] reallocOutputArray;

                                MPI_Recv(outputArray, sizeof(IndexedParticle) * sendByOther, MPI_BYTE, app.processId() - 1, 0, MPI_COMM_WORLD,0);
                            }
                        }
                    }
                }
426

427
428
429
430
431
                // we can insert into the tree
                for(int idxPart = 0 ; idxPart < outputSize ; ++idxPart){
                    treeInterval.insert(outputArray[idxPart].particle);
                }
                myNbParticlesCounter = outputSize;
432
            }
433
434
435
436
437
438
            else{
                //////////////////////////////////////////////////////////////////////////////////
                // My method
                //////////////////////////////////////////////////////////////////////////////////
                FVector<ParticlesGroup> groups;
                ParticleClass*const realParticles = reinterpret_cast<ParticleClass*>(new char[loader.getNumberOfParticles() * sizeof(ParticleClass)]);
439

440
441
                {
                    OctreeClass sortingTree(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
442

443
444
445
446
447
448
                    ParticleClass particle;
                    for(long idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
                        loader.fillParticle(particle);
                        //printf("%f %f %f\n",particle.getPosition().getX(),particle.getPosition().getY(),particle.getPosition().getZ());
                        sortingTree.insert(particle);
                    }
449

450
451
452
                    //////////////////////////////////////////////////////////////////////////////////
                    //////////////////////////////////////////////////////////////////////////////////
                    int indexPart = 0;
453

454
455
456
457
458
                    OctreeClass::Iterator octreeIterator(&sortingTree);
                    octreeIterator.gotoBottomLeft();
                    do{
                        ContainerClass::ConstBasicIterator iter(*octreeIterator.getCurrentListTargets());
                        const MortonIndex indexAtThisLeaf = octreeIterator.getCurrentGlobalIndex();
459

460
                        groups.push(ParticlesGroup(octreeIterator.getCurrentListTargets()->getSize(),indexPart, indexAtThisLeaf));
461

462
463
                        while( iter.hasNotFinished() ){
                            realParticles[indexPart] = iter.data();
464

465
                            //std::cout << "Particles with index " << indexPart << " has a morton index of " << indexAtThisLeaf << std::endl;
466

467
468
                            //const F3DPosition& particlePosition = realParticles[indexPart].getPosition();
                            //std::cout << "\t The real position of this particle is (" << particlePosition.getX() << ";" << particlePosition.getY() << ";" << particlePosition.getZ() << ")" << std::endl;
469

470
471
472
473
                            ++indexPart;
                            iter.gotoNext();
                        }
                    } while(octreeIterator.moveRight());
474

475
                }
476

477
478
                counter.tac();
                std::cout << "Done  " << "(@Creating and Inserting Temporary Particles = " << counter.elapsed() << "s)." << std::endl;
479
480


481
482
483
                //////////////////////////////////////////////////////////////////////////////////
                // We send the particle that do not belong to us
                //////////////////////////////////////////////////////////////////////////////////
484

485
486
                const MortonIndex min = app.broadcast( app.reduceMin(groups[0].index) );
                const MortonIndex max = app.broadcast( app.reduceMax(groups[groups.getSize() - 1].index) );
487

488
489
                const MortonIndex startIndex = app.getLeft(max - min + 1) + min;
                const MortonIndex endIndex = app.getRight(max - min + 1) + min;
490

491
492
493
                printf("%d we are going from %lld to %lld\n",app.processId(), min, max);
                printf("%d I will go from %lld to %lld\n",app.processId(), startIndex, endIndex);
                printf("There is actually %d leafs\n", groups.getSize());
494

495
496
                int*const needToReceive = new int[app.processCount() * app.processCount()];
                memset(needToReceive,0,app.processCount() * app.processCount() * sizeof(int));
497

498
499
500
501
                FMpi::Request requests[app.processCount()];
                {
                    int needToSend[app.processCount()];
                    memset(needToSend, 0, sizeof(int) * app.processCount());
502

503
504
505
506
507
                    MortonIndex rightMortonIndex = min;
                    int groudIndex = 0;
                    for(int idxProc = 0 ; idxProc < app.processCount() && groudIndex < groups.getSize() ; ++idxProc){
                        rightMortonIndex = app.getOtherRight(max - min + 1, idxProc) + min;
                        printf("Working with %d, he goes to %lld \n",idxProc,rightMortonIndex);
508

509
510
511
512
513
514
515
516
                        if(idxProc != app.processId()){
                            int size = 0;
                            int currentGroupIndex = groudIndex;
                            while(groudIndex < groups.getSize() && groups[groudIndex].index < rightMortonIndex){
                                size += groups[groudIndex].number;
                                ++groudIndex;
                            }
                            needToSend[idxProc] = size;
517

518
519
520
521
522
523
524
525
526
527
528
529
530
                            printf("%d Send %d to %d\n",app.processId(), size, idxProc);
                            app.isendData( idxProc, sizeof(ParticleClass) * size, &realParticles[groups[currentGroupIndex].positionInArray], 1, &requests[idxProc]);
                        }
                        else{
                            needToSend[idxProc] = 0;
                            while(groudIndex < groups.getSize() && groups[groudIndex].index < rightMortonIndex){
                                const int end = groups[groudIndex].positionInArray + groups[groudIndex].number;
                                for(int idxPart = groups[groudIndex].positionInArray ; idxPart < end ; ++idxPart){
                                    //std::cout << "\t I keep (" << realParticles[idxPart].getPosition().getX() << ";" << realParticles[idxPart].getPosition().getY() << ";" << realParticles[idxPart].getPosition().getZ() << ")" << std::endl;
                                    treeInterval.insert(realParticles[idxPart]);
                                    ++myNbParticlesCounter;
                                }
                                ++groudIndex;
531
532
                            }
                        }
533
                    }
534

535
536
537
538
539
                    app.allgather(needToSend, app.processCount(), needToReceive, app.processCount());
                    for(int idxSrc = 0 ; idxSrc < app.processCount() ; ++idxSrc){
                        for(int idxTest = 0 ; idxTest < app.processCount() ; ++idxTest){
                            printf("[%d][%d] = %d\n", idxSrc, idxTest, needToReceive[idxSrc * app.processCount() + idxTest]);
                        }
540
541
542
543
                    }
                }


544
545
546
547
548
549
550
551
552
553
554
555
                //////////////////////////////////////////////////////////////////////////////////
                // We receive others particles and insert them in the tree
                //////////////////////////////////////////////////////////////////////////////////
                int CounterProcToReceive(0);
                int maxPartToReceive(0);
                for(int idxProc = 0 ; idxProc < app.processCount() ; ++idxProc){
                    if(idxProc != app.processId() && needToReceive[app.processCount() * idxProc + app.processId()]){
                        ++CounterProcToReceive;
                        if(maxPartToReceive < needToReceive[app.processCount() * idxProc + app.processId()]){
                            maxPartToReceive = needToReceive[app.processCount() * idxProc + app.processId()];
                        }
                        printf("needToReceive[idxProc][app.processId()] %d",needToReceive[app.processCount() * idxProc + app.processId()]);
556
557
                    }
                }
558

559
                printf("maxPartToReceive %d \n",maxPartToReceive);
560

561
562
563
564
565
566
                ParticleClass*const iterParticles = reinterpret_cast<ParticleClass*>(new char[maxPartToReceive * sizeof(ParticleClass)]);
                // we receive message from nb proc - 1 (from every other procs
                for(int idxProc = 0 ; idxProc < CounterProcToReceive ; ++idxProc){
                    int source(0);
                    printf("Wait data to receive\n");
                    app.waitMessage(&source);
567

568
569
                    const int nbPartFromProc = needToReceive[app.processCount() * source + app.processId()];
                    int received(0);
570

571
572
                    printf("%d receive %d\n",source,nbPartFromProc);
                    app.receiveDataFromTagAndSource(sizeof(ParticleClass) * nbPartFromProc, 1, source, iterParticles,&received);
573

574
                    printf("Received %d part*partcileSize %d \n",received,sizeof(ParticleClass) * nbPartFromProc);
575

576
577
578
579
580
581
                    printf("Insert into tree\n");
                    for(int idxPart = 0 ; idxPart < nbPartFromProc ; ++idxPart){
                        //std::cout << "\t We receive a new particle (" << (*iterParticles).getPosition().getX() << ";" << (*iterParticles).getPosition().getY() << ";" << (*iterParticles).getPosition().getZ() << ")" << std::endl;
                        treeInterval.insert(iterParticles[idxPart]);
                        ++myNbParticlesCounter;
                    }
582
                }
583

584
585
586
587
588
                printf("Wait all send\n");
                for(int idxProc = 0 ; idxProc < app.processCount() ; ++idxProc){
                    if(idxProc != app.processId() && needToReceive[app.processCount() * app.processId() + idxProc ]){
                        app.iWait(&requests[idxProc]);
                    }
589
590
                }

591
592
593
594
                printf("Delete particle array\n");
                delete [] reinterpret_cast<char*>(realParticles);
                delete [] needToReceive;
            }
595
596
597
598
599
600
601
602
603
604
605
        }

        //////////////////////////////////////////////////////////////////////////////////
        // Now we can build the real tree
        //////////////////////////////////////////////////////////////////////////////////


        {
            //////////////////////////////////////////////////////////////////////////////////
            // We inform the master proc about the data we have
            //////////////////////////////////////////////////////////////////////////////////
606
            printf("Inform other about leaves we have\n");
607
608

            FVector<ParticlesGroup> groups;
609
            ParticleClass*const realParticles = myNbParticlesCounter?new ParticleClass[myNbParticlesCounter]:0;
610
611
612
613
614
615
616

            int nbLeafs = 0;

            // we might now have any particles in our interval
            if(myNbParticlesCounter){
                OctreeClass::Iterator octreeIterator(&treeInterval);
                octreeIterator.gotoBottomLeft();
617
                int indexPart = 0;
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
                do{
                    ContainerClass::ConstBasicIterator iter(*octreeIterator.getCurrentListTargets());
                    const MortonIndex indexAtThisLeaf = octreeIterator.getCurrentGlobalIndex();

                    groups.push(ParticlesGroup(octreeIterator.getCurrentListTargets()->getSize(),indexPart, indexAtThisLeaf));

                    while( iter.hasNotFinished() ){
                        realParticles[indexPart] = iter.data();
                        ++indexPart;
                        iter.gotoNext();
                    }

                    ++nbLeafs;
                } while(octreeIterator.moveRight());
            }

            printf("%d I have %d leafs\n",app.processId(), nbLeafs);

            // receive from left and right
            if(app.isMaster()){
                app.sendData( 1, sizeof(int), &nbLeafs, 0);
            }
            else if(app.processId() == app.processCount() - 1){
                app.sendData( app.processId() - 1, sizeof(int), &nbLeafs, 0);
            }
            // receive
            int leftLeafs = 0;
            int rightLeafs = 0;
            if(!app.isMaster() && app.processId() != app.processCount() - 1){
                for(int idxToReceive = 0 ; idxToReceive < 2 ; ++idxToReceive){
                    int source(0);
                    int temp = 0;
                    app.receiveDataFromTag(sizeof(int), 0, &temp, &source);
                    if(source < app.processId()){ // come from left
                        leftLeafs = temp;
                        temp += nbLeafs;
                        app.sendData( app.processId() + 1, sizeof(int), &temp, 0);
                    }
                    else { // come from right
                        rightLeafs = temp;
                        temp += nbLeafs;
                        app.sendData( app.processId() - 1, sizeof(int), &temp, 0);
                    }
                }
            }
            else {
                if(app.isMaster()){ // come from right
                    app.receiveDataFromTag(sizeof(int), 0, &rightLeafs);
                }
                else { // come from left
                    app.receiveDataFromTag(sizeof(int), 0, &leftLeafs);
                }
            }

672
            printf("%d I have %d leafs on the left and %d on the right\n",app.processId(), leftLeafs, rightLeafs);
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

            //////////////////////////////////////////////////////////////////////////////////
            // We balance the data
            //////////////////////////////////////////////////////////////////////////////////

            const int totalNbLeafs = (leftLeafs + nbLeafs + rightLeafs);
            const int myLeftLeaf = app.getLeft(totalNbLeafs);
            const int myRightLeaf = app.getRight(totalNbLeafs);

            const bool iNeedToSendToLeft = leftLeafs < myLeftLeaf;
            const bool iNeedToSendToRight = myRightLeaf < leftLeafs + nbLeafs;

            const bool iWillReceiveFromRight = leftLeafs + nbLeafs < myRightLeaf;
            const bool iWillReceiveFromLeft = leftLeafs > myLeftLeaf;

            const bool iDoNotHaveEnoughtToSendRight = myRightLeaf < leftLeafs;
            const bool iDoNotHaveEnoughtToSendLeft = leftLeafs + nbLeafs < myLeftLeaf;

            ParticleClass* rpart(0);
            int rpartSize(0);

            // Do I need to send to right?
            if(iNeedToSendToRight){
                int iNeedToSend = leftLeafs + nbLeafs - myRightLeaf;
                int iCanSend = nbLeafs;
                int idxSend (0);

                printf("%d I need to send %d to right\n",app.processId(), iNeedToSend);
                app.sendData( app.processId() + 1, sizeof(int), &iNeedToSend, 0);

                while(idxSend < iNeedToSend && idxSend < iCanSend){
                    app.sendData(app.processId() + 1, sizeof(int),&groups[nbLeafs - idxSend - 1].number,0);
                    app.sendData(app.processId() + 1, sizeof(ParticleClass) * groups[nbLeafs - idxSend - 1].number,
706
                                 &realParticles[groups[nbLeafs - idxSend - 1].positionInArray], 0);
707
708
709
710
711
                    ++idxSend;
                }
                // I need to wait (idxSend == iCanSend && idxSend < iNeedToSend)
                if( iDoNotHaveEnoughtToSendRight ){ // I need to wait
                    int nbLeafsToRead(0);
712
713
                    app.receiveDataFromTag(sizeof(int), 0, &nbLeafsToRead);
                    printf("Receive %d from left to send to right\n", nbLeafsToRead);
714
715
                    for(int idxToRead = 0 ; idxToRead < nbLeafsToRead ; ++idxToRead){
                        int nbPartToRead(0);
716
                        app.receiveDataFromTag(sizeof(int), 0, &nbPartToRead);
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

                        if(rpartSize < nbPartToRead){
                            rpartSize = nbPartToRead;
                            delete [] (reinterpret_cast<char*>(rpart));
                            rpart = reinterpret_cast<ParticleClass*>(new char[nbPartToRead*sizeof(ParticleClass)]);
                        }

                        app.receiveDataFromTag(nbPartToRead*sizeof(ParticleClass), 0, rpart);
                        if(idxSend < iNeedToSend){
                            app.sendData(app.processId() + 1, sizeof(int),&nbPartToRead,0);
                            app.sendData(app.processId() + 1, sizeof(ParticleClass) * nbPartToRead, rpart,0);
                            ++idxSend;
                        }
                        else{
                            //insert into tree
                            for(int idxPart = 0 ; idxPart < nbPartToRead ; ++idxPart){
                                realTree.insert(rpart[idxPart]);
                            }
                        }
                    }
                }
            }
            // will I receive from left
            if(iNeedToSendToLeft){
                int iNeedToSend = myLeftLeaf - leftLeafs;
                int iCanSend = nbLeafs;
                int idxSend (0);

745
                printf("%d I need to send %d to left\n",app.processId(), iNeedToSend);
746
747
748
749
750
751
752
753
754
755
756
757
                app.sendData( app.processId() - 1, sizeof(int), &iNeedToSend, 1);

                while(idxSend < iNeedToSend && idxSend < iCanSend){
                    app.sendData(app.processId() - 1, sizeof(int),&groups[idxSend].number,1);
                    app.sendData(app.processId() - 1, sizeof(ParticleClass) * groups[idxSend].number,
                                 &realParticles[groups[idxSend].positionInArray],1);
                    ++idxSend;
                }
                // Can I do it now?
                if( iDoNotHaveEnoughtToSendLeft ){
                    int nbLeafsToRead(0);
                    app.receiveDataFromTag(sizeof(int), 1, &nbLeafsToRead);
758
                    printf("Receive %d from right to send to left\n", nbLeafsToRead);
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
                    for(int idxToRead = 0 ; idxToRead < nbLeafsToRead ; ++idxToRead){
                        int nbPartToRead(0);
                        app.receiveDataFromTag(sizeof(int), 1, &nbPartToRead);

                        if(rpartSize < nbPartToRead){
                            rpartSize = nbPartToRead;
                            delete [] (reinterpret_cast<char*>(rpart));
                            rpart = reinterpret_cast<ParticleClass*>(new char[nbPartToRead*sizeof(ParticleClass)]);
                        }

                        app.receiveDataFromTag(nbPartToRead*sizeof(ParticleClass), 1, rpart);
                        if(idxSend < iNeedToSend){
                            app.sendData(app.processId() - 1, sizeof(int),&nbPartToRead,1);
                            app.sendData(app.processId() - 1, sizeof(ParticleClass) * nbPartToRead, rpart,1);
                            ++idxSend;
                        }
                        else{
                            for(int idxPart = 0 ; idxPart < nbPartToRead ; ++idxPart){
                                realTree.insert(rpart[idxPart]);
                            }
                        }
                    }
                }
            }

            // If i will receive from left and I did no already have
            if(!(iNeedToSendToRight && iDoNotHaveEnoughtToSendRight) && iWillReceiveFromLeft){
                printf("%d I need to receive from left\n",app.processId());
                int nbLeafsToRead(0);
                app.receiveDataFromTag(sizeof(int), 0, &nbLeafsToRead);
                printf("%d I will receive from left %d\n",app.processId(), nbLeafsToRead);
                for(int idxToRead = 0 ; idxToRead < nbLeafsToRead ; ++idxToRead){
                    int nbPartToRead(0);
                    app.receiveDataFromTag(sizeof(int), 0, &nbPartToRead);
793
                    //printf("%d I will receive %d particles\n",app.processId(), nbPartToRead);
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
                    if(rpartSize < nbPartToRead){
                        rpartSize = nbPartToRead;
                        delete [] (reinterpret_cast<char*>(rpart));
                        rpart = reinterpret_cast<ParticleClass*>(new char[nbPartToRead*sizeof(ParticleClass)]);
                    }

                    app.receiveDataFromTag(nbPartToRead*sizeof(ParticleClass), 0, rpart);
                    for(int idxPart = 0 ; idxPart < nbPartToRead ; ++idxPart){
                        realTree.insert(rpart[idxPart]);
                    }
                }
            }
            // If i will receive from right and I did no already have
            if(!(iNeedToSendToLeft && iDoNotHaveEnoughtToSendLeft) && iWillReceiveFromRight){
                printf("%d I need to receive from right\n",app.processId());
                int nbLeafsToRead(0);
                app.receiveDataFromTag(sizeof(int), 1, &nbLeafsToRead);
811
                printf("%d I will receive from right %d\n",app.processId(), nbLeafsToRead);
812
813
814
                for(int idxToRead = 0 ; idxToRead < nbLeafsToRead ; ++idxToRead){
                    int nbPartToRead(0);
                    app.receiveDataFromTag(sizeof(int), 1, &nbPartToRead);
815
                    //printf("%d I will receive %d particles\n",app.processId(), nbPartToRead);
816
817
818
819
820
821
822
823
824
825
826
827
828
                    if(rpartSize < nbPartToRead){
                        rpartSize = nbPartToRead;
                        delete [] (reinterpret_cast<char*>(rpart));
                        rpart = reinterpret_cast<ParticleClass*>(new char[nbPartToRead*sizeof(ParticleClass)]);
                    }

                    app.receiveDataFromTag(nbPartToRead*sizeof(ParticleClass), 1, rpart);
                    for(int idxPart = 0 ; idxPart < nbPartToRead ; ++idxPart){
                        realTree.insert(rpart[idxPart]);
                    }
                }
            }

829
830
            printf("Will now take my own particles from %d to %d\n",FMath::Max(myLeftLeaf-leftLeafs,0) , FMath::Max(myLeftLeaf-leftLeafs,0) + FMath::Min(myRightLeaf,totalNbLeafs - rightLeafs) - myLeftLeaf);
            printf("myLeftLeaf %d leftLeafs %d myRightLeaf %d rightLeafs %d totalNbLeafs %d\n",myLeftLeaf,leftLeafs, myRightLeaf, rightLeafs, totalNbLeafs);
831
            // insert the particles we already have
832
833
834
835
836
            if(leftLeafs != totalNbLeafs){
                for(int idxLeafInsert = FMath::Max(myLeftLeaf-leftLeafs,0) ; idxLeafInsert < FMath::Max(myLeftLeaf-leftLeafs,0) + FMath::Min(myRightLeaf,totalNbLeafs- rightLeafs) - myLeftLeaf ; ++idxLeafInsert){
                    for(int idxPart = 0 ; idxPart < groups[idxLeafInsert].number ; ++idxPart){
                        realTree.insert(realParticles[groups[idxLeafInsert].positionInArray + idxPart]);
                    }
837
838
839
840
841
842
843
844
845
846
847
848
                }
            }

            delete [] reinterpret_cast<char*>(rpart);
        }

        app.processBarrier();

        //////////////////////////////////////////////////////////////////////////////////
        //////////////////////////////////////////////////////////////////////////////////
    }

849
    //////////////////////////////////////////////////////////////////////////////////
850
    // Create real tree
851
852
    //////////////////////////////////////////////////////////////////////////////////

berenger-bramas's avatar
berenger-bramas committed
853
    OctreeClass treeValide(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
854
855
856
857
858
859
860
861
    {
        FFmaBinLoader<ParticleClass> loaderSeq(filename);
        ParticleClass partToInsert;
        for(int idxPart = 0 ; idxPart < loaderSeq.getNumberOfParticles() ; ++idxPart){
            loaderSeq.fillParticle(partToInsert);
            treeValide.insert(partToInsert);
        }
    }
862

863
    //////////////////////////////////////////////////////////////////////////////////
864
    // Check particles in tree
865
866
    //////////////////////////////////////////////////////////////////////////////////

867
    {
868
869
870
871
872
873
874
875
876
        int totalNbLeafs = 0;
        {

            OctreeClass::Iterator octreeIterator(&treeValide);
            octreeIterator.gotoBottomLeft();
            do{
                ++totalNbLeafs;
            }while(octreeIterator.moveRight());
        }
877

878
879
        const int myLeftLeaf = app.getLeft(totalNbLeafs);
        const int myRightLeaf = app.getRight(totalNbLeafs);
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
        printf("%d should go from %d to %d leaf (on %d total leafs)\n", app.processId(), myLeftLeaf, myRightLeaf, totalNbLeafs);

        OctreeClass::Iterator octreeIteratorValide(&treeValide);
        octreeIteratorValide.gotoBottomLeft();
        for(int idxLeaf = 0 ; idxLeaf < myLeftLeaf ; ++idxLeaf){
            if(!octreeIteratorValide.moveRight()){
                printf("Error cannot access to the left leaf %d in the valide tree\n", myLeftLeaf);
            }
        }

        OctreeClass::Iterator octreeIterator(&realTree);
        octreeIterator.gotoBottomLeft();

        for(int idxLeaf = myLeftLeaf ; idxLeaf < myRightLeaf ; ++idxLeaf){
            if(octreeIteratorValide.getCurrentGlobalIndex() != octreeIterator.getCurrentGlobalIndex()){
                printf("Error index are different, valide %lld invalid %lld\n",octreeIteratorValide.getCurrentGlobalIndex(),
                       octreeIterator.getCurrentGlobalIndex());
                break;
            }
            if(octreeIteratorValide.getCurrentListSrc()->getSize() != octreeIterator.getCurrentListSrc()->getSize()){
                printf("Error leafs do not have the same number of particles, valide %d, invalide %d\n",
                       octreeIteratorValide.getCurrentListSrc()->getSize(), octreeIterator.getCurrentListSrc()->getSize() );
            }

905
            //printf("index %lld with %d particles\n", octreeIteratorValide.getCurrentGlobalIndex(), octreeIteratorValide.getCurrentListSrc()->getSize());
906

907
908
            if(!octreeIteratorValide.moveRight() && idxLeaf != myRightLeaf - 1){
                printf("Error cannot valide tree end to early, idxLeaf %d myRightLeaf %d\n", idxLeaf, myRightLeaf);
909
910
                break;
            }
911
912
913

            if(!octreeIterator.moveRight() && idxLeaf != myRightLeaf - 1){
                printf("Error cannot test tree end to early, idxLeaf %d myRightLeaf %d\n", idxLeaf, myRightLeaf);
914
915
                break;
            }
916
        }
917
918
919



920
    }
921
    return 0;
922
923
924
925
926
927
928

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    std::cout << "Working on particles ..." << std::endl;
    counter.tic();

berenger-bramas's avatar
berenger-bramas committed
929
    KernelClass kernels;
berenger-bramas's avatar
berenger-bramas committed
930

931
    FmmClassProc algo(app,&realTree,&kernels);
932
933
    algo.execute();

berenger-bramas's avatar
berenger-bramas committed
934
    FmmClass algoValide(&treeValide,&kernels);
berenger-bramas's avatar
berenger-bramas committed
935
936
    algoValide.execute();

937
    counter.tac();
938
    std::cout << "Done  " << "(@Algorithm Particles = " << counter.elapsed() << "s)." << std::endl;
939
940
941
942

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

943
    ValidateFMMAlgoProc<OctreeClass,ContainerClass,FmmClassProc>(&realTree,&treeValide,&algo);
944
945
946
947
948
949
950
951
952

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    return 0;
}


// [--LICENSE--]