testFmmAlgorithmProc.cpp 20.1 KB
Newer Older
1
// See LICENCE file at project root
2

3 4 5 6
// ==== CMAKE =====
// @FUSE_MPI
// ================

7 8
#include "../../Src/Utils/FMpi.hpp"
#include "../../Src/Utils/FTic.hpp"
9

10 11 12 13
#include "../../Src/Containers/FOctree.hpp"
#include "../../Src/Containers/FVector.hpp"
#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Utils/FGlobal.hpp"
14

15
#include "../../Src/Components/FSimpleLeaf.hpp"
16

17
#include "../../Src/Utils/FPoint.hpp"
18

19 20
#include "../../Src/Components/FTestCell.hpp"
#include "../../Src/Components/FTestKernels.hpp"
21
#include "../../Src/Components/FTestParticleContainer.hpp"
22

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

COULAUD Olivier's avatar
COULAUD Olivier committed
27
#include "../../Src/Files/FMpiFmaGenericLoader.hpp"
28
#include "../../Src/Files/FMpiTreeBuilder.hpp"
29

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

32
#include "../../Src/Utils/FLeafBalance.hpp"
33

34 35
#include "../../Src/Utils/FParameterNames.hpp"

berenger-bramas's avatar
berenger-bramas committed
36
#include <iostream>
37 38
#include <cstdio>
#include <cstdlib>
berenger-bramas's avatar
berenger-bramas committed
39 40


41
/** This program show an example of use of the fmm threaded + mpi algo
42 43
 * it also check that each particles is impacted each other particles
 */
44

45 46 47
/////////////////////////////////////////////////////////////////////////////
// Test function
/////////////////////////////////////////////////////////////////////////////
48

49 50 51
// Check if tree is built correctly
template<class OctreeClass>
void ValidateTree(OctreeClass& realTree,
52
                  OctreeClass& treeValide){
53

54 55
    typename OctreeClass::Iterator octreeIteratorValide(&treeValide);
    octreeIteratorValide.gotoBottomLeft();
56

57 58
    typename OctreeClass::Iterator octreeIterator(&realTree);
    octreeIterator.gotoBottomLeft();
59

60 61 62 63 64 65 66
    while(octreeIteratorValide.getCurrentGlobalIndex() != octreeIterator.getCurrentGlobalIndex()){
        if(octreeIteratorValide.moveRight() == false){
            std::cout << "Error the real tree smaller than the parallel one\n";
            std::cout << "Valide tree stop at " << octreeIteratorValide.getCurrentGlobalIndex() << "\n";
            std::cout << "Other at " << octreeIterator.getCurrentGlobalIndex() << "\n";
            return;
        }
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
    std::cout << "The tree starts at " << octreeIteratorValide.getCurrentGlobalIndex() << "\n";

    while(true){
        if(octreeIteratorValide.getCurrentGlobalIndex() != octreeIterator.getCurrentGlobalIndex()){
            std::cout << "Error the trees do not have the same idx.\n";
            std::cout << "Correct one is " << octreeIteratorValide.getCurrentGlobalIndex() << "\n";
            std::cout << "Incorrect one is " << octreeIterator.getCurrentGlobalIndex() << "\n";
            return;
        }

        if(octreeIteratorValide.getCurrentListSrc()->getNbParticles() != octreeIterator.getCurrentListSrc()->getNbParticles()){
            std::cout << "Error the trees do not have the same number of particles at idx " << octreeIteratorValide.getCurrentGlobalIndex() << ".\n";
            std::cout << "Correct one is " << octreeIteratorValide.getCurrentListSrc()->getNbParticles() << "\n";
            std::cout << "Incorrect one is " << octreeIterator.getCurrentListSrc()->getNbParticles() << "\n";
            return;
        }

        if(octreeIterator.moveRight() == false){
            break;
        }

        if(octreeIteratorValide.moveRight() == false){
            std::cout << "Error the real tree smaller than the parallel one\n";
        }
93
    }
94

95
    std::cout << "The tree stops at " << octreeIteratorValide.getCurrentGlobalIndex() << "\n";
96
}
97 98


99

berenger-bramas's avatar
berenger-bramas committed
100
/** This function tests the octree to be sure that the fmm algorithm
101 102
 * has worked completly.
 */
103
template<class OctreeClass, class ContainerClass, class FmmClassProc>
berenger-bramas's avatar
berenger-bramas committed
104
void ValidateFMMAlgoProc(OctreeClass* const badTree,
105 106 107 108
                         OctreeClass* const valideTree,
                         FmmClassProc* const fmm){
    std::cout << "The working interval is from " << fmm->getWorkingInterval(badTree->getHeight()-1).leftIndex << "\n";
    std::cout << "The working interval is to " << fmm->getWorkingInterval(badTree->getHeight()-1).rightIndex << "\n";
109
    std::cout << "\tValidate L2L M2M...\n";
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
    const int OctreeHeight = badTree->getHeight();
    {
        typename OctreeClass::Iterator octreeIterator(badTree);
        octreeIterator.gotoBottomLeft();

        typename OctreeClass::Iterator octreeIteratorValide(valideTree);
        octreeIteratorValide.gotoBottomLeft();

        for(int level = OctreeHeight - 1 ; level > 0 && fmm->hasWorkAtLevel(level) ; --level){

            while(octreeIteratorValide.getCurrentGlobalIndex() != octreeIterator.getCurrentGlobalIndex()) {
                octreeIteratorValide.moveRight();
            }

            while(octreeIteratorValide.getCurrentGlobalIndex() != fmm->getWorkingInterval(level).leftIndex){
                octreeIteratorValide.moveRight();
                octreeIterator.moveRight();
            }

            FSize countCheck = 0;
            do{
                if(octreeIterator.getCurrentGlobalIndex() != octreeIteratorValide.getCurrentGlobalIndex()){
                    std::cout << "Problem Error index are not equal! " << octreeIterator.getCurrentGlobalIndex() << " " << octreeIteratorValide.getCurrentGlobalIndex() << std::endl;
                }
                else{
135
                    if(octreeIterator.getCurrentCell()->getDataUp().get() != octreeIteratorValide.getCurrentCell()->getDataUp().get()){
136 137 138
                        std::cout << "Problem M2M error at level " << level << " up bad " << octreeIterator.getCurrentCell()->getDataUp()
                                  << " good " << octreeIteratorValide.getCurrentCell()->getDataUp() << " index " << octreeIterator.getCurrentGlobalIndex() << std::endl;
                    }
139
                    if(octreeIterator.getCurrentCell()->getDataDown().get() != octreeIteratorValide.getCurrentCell()->getDataDown().get()){
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
                        std::cout << "Problem L2L error at level " << level << " down bad " << octreeIterator.getCurrentCell()->getDataDown()
                                  << " good " << octreeIteratorValide.getCurrentCell()->getDataDown() << " index " << octreeIterator.getCurrentGlobalIndex() << std::endl;
                    }
                }
                ++countCheck;
            } while(octreeIteratorValide.moveRight() && octreeIterator.moveRight());

            // Check that each particle has been summed with all other

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

            octreeIteratorValide.moveUp();
            octreeIteratorValide.gotoLeft();
        }
155
    }
156
    std::cout << "\tValidate L2P P2P...\n";
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
        FSize NbPart = 0;
        FSize NbLeafs = 0;
        { // Check that each particle has been summed with all other
            typename OctreeClass::Iterator octreeIterator(valideTree);
            octreeIterator.gotoBottomLeft();
            do{
                NbPart += octreeIterator.getCurrentListSrc()->getNbParticles();
                ++NbLeafs;
            } while(octreeIterator.moveRight());
        }
        {
            // Check that each particle has been summed with all other
            typename OctreeClass::Iterator octreeIterator(badTree);
            octreeIterator.gotoBottomLeft();

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

                ContainerClass* container = (octreeIterator.getCurrentListTargets());
                const long long int*const dataDown = container->getDataDown();

                for(FSize idxPart = 0 ; idxPart < container->getNbParticles() ; ++idxPart){
                    // If a particles has been impacted by less than NbPart - 1 (the current particle)
                    // there is a problem
                    if( (!isUsingTsm && dataDown[idxPart] != NbPart - 1) ||
                            (isUsingTsm && dataDown[idxPart] != NbPart) ){
                        std::cout << "Problem L2P + P2P, value on particle is : " << dataDown[idxPart] <<
                                     " at pos " << idxPart << " index is " << octreeIterator.getCurrentGlobalIndex() << "\n";
                    }
                }
            } while( octreeIterator.moveRight());
        }
190
    }
191
    std::cout << "\tValidate P2M...\n";
192
    {
193 194 195 196 197 198 199
        {
            // Check that each particle has been summed with all other
            typename OctreeClass::Iterator octreeIterator(badTree);
            octreeIterator.gotoBottomLeft();

            do {
                if(octreeIterator.getCurrentListSrc()->getNbParticles() != octreeIterator.getCurrentCell()->getDataUp()){
200
                    printf("P2M Problem nb part %lld data up %lld \n",
201 202 203 204
                           octreeIterator.getCurrentListSrc()->getNbParticles(), octreeIterator.getCurrentCell()->getDataUp());
                }
            } while( octreeIterator.moveRight() );
        }
205
    }
206
    std::cout << "\tValidate Particles...\n";
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
    {
        // Check that each particle has been summed with all other
        typename OctreeClass::Iterator octreeIterator(badTree);
        octreeIterator.gotoBottomLeft();

        typename OctreeClass::Iterator valideOctreeIterator(valideTree);
        valideOctreeIterator.gotoBottomLeft();
        while(valideOctreeIterator.getCurrentGlobalIndex() != octreeIterator.getCurrentGlobalIndex()){
            valideOctreeIterator.moveRight();
        }

        do {
            if(valideOctreeIterator.getCurrentGlobalIndex() != octreeIterator.getCurrentGlobalIndex()){
                printf("Problem Do not have the same index valide %lld invalide %lld \n",
                       valideOctreeIterator.getCurrentGlobalIndex(), octreeIterator.getCurrentGlobalIndex());
                break;
            }

            if(octreeIterator.getCurrentListTargets()->getNbParticles() != valideOctreeIterator.getCurrentListTargets()->getNbParticles()){
226
                printf("Problem Do not have the same number of particle at leaf id %lld, valide %lld invalide %lld \n",
227 228 229 230 231 232 233 234 235 236
                       octreeIterator.getCurrentGlobalIndex(), valideOctreeIterator.getCurrentListTargets()->getNbParticles(),
                       octreeIterator.getCurrentListTargets()->getNbParticles());
            }
            else {
                ContainerClass* container = (octreeIterator.getCurrentListTargets());
                const long long int*const dataDown = container->getDataDown();

                ContainerClass* containerValide = (valideOctreeIterator.getCurrentListTargets());
                const long long int*const dataDownValide = containerValide->getDataDown();

237
                for(FSize idxPart = 0 ; idxPart < container->getNbParticles() ; ++idxPart){
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
                    // If a particles has been impacted by less than NbPart - 1 (the current particle)
                    // there is a problem
                    if( dataDown[idxPart] != dataDownValide[idxPart]){
                        std::cout << "Problem on leaf " << octreeIterator.getCurrentGlobalIndex() <<
                                     " part " << idxPart << " valide data down " << dataDownValide[idxPart] <<
                                     " invalide " << dataDown[idxPart] << "\n";
                        std::cout << "Data down for leaf is: valide " << valideOctreeIterator.getCurrentCell()->getDataDown()
                                  << " invalide " << octreeIterator.getCurrentCell()->getDataDown()
                                  << " size is: valide " <<  valideOctreeIterator.getCurrentListTargets()->getNbParticles()
                                  << " invalide " << octreeIterator.getCurrentListTargets()->getNbParticles() << std::endl;
                    }
                }
            }

        }while( octreeIterator.moveRight() && valideOctreeIterator.moveRight());
253
    }
254
    std::cout << "\tDone!\n";
255 256
}

257

258
/** To print an octree
259 260
 * used to debug and understand how the values were passed
 */
261 262
template<class OctreeClass>
void print(OctreeClass* const valideTree){
263 264 265 266 267 268 269 270 271
    typename OctreeClass::Iterator octreeIterator(valideTree);
    for(int idxLevel = valideTree->getHeight() - 1 ; idxLevel > 1 ; --idxLevel ){
        do{
            std::cout << "[" << octreeIterator.getCurrentGlobalIndex() << "] up:" << octreeIterator.getCurrentCell()->getDataUp() << " down:" << octreeIterator.getCurrentCell()->getDataDown() << "\t";
        } while(octreeIterator.moveRight());
        std::cout << "\n";
        octreeIterator.gotoLeft();
        octreeIterator.moveDown();
    }
272 273
}

274

berenger-bramas's avatar
berenger-bramas committed
275 276 277 278
/////////////////////////////////////////////////////////////////////
// Define the classes to use
/////////////////////////////////////////////////////////////////////

279 280
typedef double FReal;

281
typedef FTestCell                  CellClass;
282
typedef FTestParticleContainer<FReal>     ContainerClass;
283

284 285
typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
286
typedef FTestKernels< CellClass, ContainerClass >         KernelClass;
berenger-bramas's avatar
berenger-bramas committed
287

288 289
typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass >     FmmClass;
typedef FFmmAlgorithmThreadProc<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass >     FmmClassProc;
berenger-bramas's avatar
berenger-bramas committed
290

291 292 293
/////////////////////////////////////////////////////////////////////
// Main
/////////////////////////////////////////////////////////////////////
294 295 296

// Simply create particles and try the kernels
int main(int argc, char ** argv){
297 298 299 300
    FHelpDescribeAndExit(argc, argv,
                         "Test FMM distributed algorithm by counting the nb of interactions each particle receive.",
                         FParameterDefinitions::OctreeHeight, FParameterDefinitions::OctreeSubHeight,
                         FParameterDefinitions::InputFile);
301 302 303 304 305 306
    ///////////////////////What we do/////////////////////////////
    std::cout << ">> This executable has to be used to test the FMM algorithm.\n";
    //////////////////////////////////////////////////////////////

    FMpi app( argc, argv);

307 308
    const int NbLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 5);
    const int SizeSubLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 3);
309 310 311 312
    FTic counter;
    const char* const defaultFilename = (sizeof(FReal) == sizeof(float))?
                "../../Data/test20k.bin.fma.single":
                "../../Data/test20k.bin.fma.double";
313
    const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFilename);
314 315
    std::cout << "Opening : " << filename << "\n";

316
    FMpiFmaGenericLoader<FReal> loader(filename,app.global());
317
    if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");
318

319
    const FSize nbPartsForMe = loader.getMyNumberOfParticles();
320
    const FReal boxWidth = loader.getBoxWidth();
321
    const FPoint<FReal> centerOfBox = loader.getCenterOfBox();
322

323
    std::cout << "Simulation properties :\n";
324
    std::cout << "Nb Particles For me " << nbPartsForMe << "\n";
325 326
    std::cout << "Box Width : " << boxWidth << "\n";
    std::cout << "Box Center : " << centerOfBox << "\n";
327 328

    // The real tree to work on
329
    OctreeClass realTree(NbLevels, SizeSubLevels,boxWidth,centerOfBox);
330 331 332 333 334 335 336 337 338

    if( app.global().processCount() != 1){
        //////////////////////////////////////////////////////////////////////////////////
        // Build tree from mpi loader
        //////////////////////////////////////////////////////////////////////////////////
        std::cout << "Build Tree ..." << std::endl;
        counter.tic();

        struct TestParticle{
339 340
            FPoint<FReal> position;
            const FPoint<FReal>& getPosition(){
341 342 343 344
                return position;
            }
        };

345 346
        TestParticle* particles = new TestParticle[nbPartsForMe];
        memset(particles, 0, sizeof(TestParticle) * nbPartsForMe);
347
        for(FSize idxPart = 0 ; idxPart < nbPartsForMe ; ++idxPart){
348
            FPoint<FReal> position;
349 350 351
            FReal physicalValue;
            loader.fillParticle(&position, &physicalValue);
            particles[idxPart].position = position;
352
        }
353

354 355
        FVector<TestParticle> finalParticles;
        FLeafBalance balancer;
356
        FMpiTreeBuilder< FReal,TestParticle >::DistributeArrayToContainer(app.global(),particles,
357
                                                                    nbPartsForMe,
358 359 360 361
                                                                    realTree.getBoxCenter(),
                                                                    realTree.getBoxWidth(),realTree.getHeight(),
                                                                    &finalParticles, &balancer);
        std::cout << "I have now " << finalParticles.getSize() << " particles\n";
362

363 364 365
        for(int idx = 0 ; idx < finalParticles.getSize(); ++idx){
            realTree.insert(finalParticles[idx].position);
        }
366

367
        delete[] particles;
berenger-bramas's avatar
berenger-bramas committed
368

369 370 371 372 373 374
        counter.tac();
        std::cout << "Done  " << "(" << counter.elapsed() << "s)." << std::endl;

        //////////////////////////////////////////////////////////////////////////////////
    }
    else{
375
        for(FSize idxPart = 0 ; idxPart < nbPartsForMe ; ++idxPart){
376
            FPoint<FReal> position;
377 378 379
            FReal physicalValue;
            loader.fillParticle(&position, &physicalValue);
            realTree.insert(position);
380
        }
381
    }
382

383 384 385
    //////////////////////////////////////////////////////////////////////////////////
    // Create real tree
    //////////////////////////////////////////////////////////////////////////////////
386

387
    OctreeClass treeValide(NbLevels, SizeSubLevels,boxWidth,centerOfBox);
388
    {
389
        FFmaGenericLoader<FReal> loaderValide(filename);
390 391
        if(!loaderValide.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");

392
        const FSize nbPartsValide = loaderValide.getNumberOfParticles();
393
        const FReal boxWidthValide = loaderValide.getBoxWidth();
394
        const FPoint<FReal> centerOfBoxValide = loaderValide.getCenterOfBox();
395 396

        std::cout << "Simulation properties :\n";
397 398 399
        std::cout << "Nb Particles " << nbPartsValide << "\n";
        std::cout << "Box Width : " << boxWidthValide << "\n";
        std::cout << "Box Center : " << centerOfBoxValide << "\n";
400

401
        for(FSize idxPart = 0 ; idxPart < nbPartsValide ; ++idxPart){
402
            FPoint<FReal> position;
403 404 405
            FReal physicalValue;
            loaderValide.fillParticle(&position, &physicalValue);
            treeValide.insert(position);
406 407
        }
    }
408

409 410 411 412 413
    //////////////////////////////////////////////////////////////////////////////////
    // Check particles in tree
    //////////////////////////////////////////////////////////////////////////////////
    std::cout << "Validate tree ..." << std::endl;
    counter.tic();
414

415
    ValidateTree(realTree, treeValide);
berenger-bramas's avatar
berenger-bramas committed
416

417 418
    counter.tac();
    std::cout << "Done  " << "(" << counter.elapsed() << "s)." << std::endl;
419

420
    //////////////////////////////////////////////////////////////////////////////////
421

422 423
    std::cout << "Working parallel particles ..." << std::endl;
    counter.tic();
424

425
    KernelClass kernels;
426

427 428
    FmmClassProc algo(app.global(),&realTree,&kernels);
    algo.execute();
429

430 431
    counter.tac();
    std::cout << "Done  " << "(@Algorithm Particles = " << counter.elapsed() << "s)." << std::endl;
432

433
    //////////////////////////////////////////////////////////////////////////////////
434

435 436
    std::cout << "Working sequential particles ..." << std::endl;
    counter.tic();
437

438 439
    FmmClass algoValide(&treeValide,&kernels);
    algoValide.execute();
440

441 442
    counter.tac();
    std::cout << "Done  " << "(@Algorithm Particles = " << counter.elapsed() << "s)." << std::endl;
443

444 445 446 447 448 449 450
    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    std::cout << "Checking data ..." << std::endl;
    counter.tic();

    ValidateFMMAlgoProc<OctreeClass,ContainerClass, FmmClassProc>(&realTree,&treeValide,&algo);
451

452 453
    counter.tac();
    std::cout << "Done  " << "(" << counter.elapsed() << "s)." << std::endl;
454

455 456
    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////
457

458
    return 0;
459
}