testFmmAlgorithmPeriodic.cpp 6.06 KB
Newer Older
1
// See LICENCE file at project root
2 3 4 5


#include <iostream>

6 7
#include <cstdio>
#include <cstdlib>
8

9 10
#include "Utils/FParameters.hpp"
#include "Utils/FTic.hpp"
11

12
#include "Files/FRandomLoader.hpp"
13

14
#include "Files/FPerLeafLoader.hpp"
15

16 17
#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"
18

19
#include "Components/FSimpleLeaf.hpp"
20

21
#include "Components/FTestParticleContainer.hpp"
22

23
#include "Utils/FPoint.hpp"
24

25 26
#include "Components/FTestCell.hpp"
#include "Components/FTestKernels.hpp"
27

28
#include "Core/FFmmAlgorithmPeriodic.hpp"
29

30
#include "Utils/FParameterNames.hpp"
31

32 33 34 35 36 37 38 39
/** This program show an example of use of
  * the fmm basic algo
  * it also check that each particles is impacted each other particles
  */


// Simply create particles and try the kernels
int main(int argc, char ** argv){
40 41 42 43 44
    FHelpDescribeAndExit(argc, argv,
                         "Test FMM periodic algorithm by counting the nb of interactions each particle receive.",
                         FParameterDefinitions::OctreeHeight, FParameterDefinitions::OctreeSubHeight,
                         FParameterDefinitions::NbParticles, FParameterDefinitions::PeriodicityNbLevels);

45
    typedef double FReal;
46
    typedef FTestCell                   CellClass;
47
    typedef FTestParticleContainer<FReal>      ContainerClass;
48

49 50
    typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
    typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
51
    typedef FTestKernels< CellClass, ContainerClass >         KernelClass;
52

53
    typedef FFmmAlgorithmPeriodic<FReal, OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass >     FmmClass;
54 55 56 57
    ///////////////////////What we do/////////////////////////////
    std::cout << ">> This executable has to be used to test the FMM algorithm.\n";
    //////////////////////////////////////////////////////////////

58 59 60 61
    const int NbLevels          = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 5);
    const int SizeSubLevels     = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 3);
    const long NbParticles      = FParameters::getValue(argc,argv,FParameterDefinitions::NbParticles.options, 1000);
    const int PeriodicDeep      = FParameters::getValue(argc,argv,FParameterDefinitions::PeriodicityNbLevels.options, 2);
62
    // choose in +x dir or -/+x dir or all dirs
63 64 65 66 67 68

    FTic counter;

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

69
    FRandomLoader<FReal> loader(NbParticles);
70 71 72
    //FPerLeafLoader loader(NbLevels);

    std::cout << "Creating & Inserting " << loader.getNumberOfParticles() << " particles ..." << std::endl;
73 74 75
    std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
    counter.tic();

76
    OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
77 78

    {
79
        FPoint<FReal> particlePosition;
80
        for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
81 82 83 84
            loader.fillParticle(&particlePosition);
            tree.insert(particlePosition);
        }
    }
85 86 87 88 89 90 91 92 93 94 95

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

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

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

    KernelClass kernels;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
96
    FmmClass algo( &tree, PeriodicDeep);
97
    algo.setKernel(&kernels);
98 99 100 101 102 103 104 105
    algo.execute();

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

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

106
    { // Check that each particle has been summed with all other
107
        long long int counterNbPart = 0;
108

109
        tree.forEachCellLeaf([&](CellClass* cell, LeafClass* leaf){
110 111
            if(cell->getMultipoleData().get() != leaf->getSrc()->getNbParticles() ){
                    std::cout << "Problem P2M Data up = " << cell->getMultipoleData().get() <<
112
                                 " Size = " << leaf->getSrc()->getNbParticles() << "\n";
113 114
            }
            // we also count the number of particles.
115 116
            counterNbPart += leaf->getSrc()->getNbParticles();
        });
117

118 119 120
        if( counterNbPart != loader.getNumberOfParticles()){
            std::cout << "Problem global nb part, counter = " << counterNbPart << " created = " <<
                         loader.getNumberOfParticles() << std::endl;
121 122
        }
    }
123
    {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
124 125 126
        const long repetitions = algo.theoricalRepetition();
        const long totalRepeatedBox = repetitions * repetitions * repetitions;
        std::cout << "The box is repeated " << repetitions << " there are " << totalRepeatedBox << " boxes in total\n";
127
        const long long NbParticlesEntireSystem = loader.getNumberOfParticles() * totalRepeatedBox;
128
        std::cout << "The total number of particles is "  << NbParticlesEntireSystem << "\n";
129

130
        FTreeCoordinate min, max;
131
        algo.repetitionsIntervals(&min, &max);
132
        std::cout << "Min is " << min << " Max is " << max << std::endl;
133

134
        tree.forEachLeaf([&](LeafClass* leaf){
135
            for(FSize idxPart = 0 ; idxPart < leaf->getSrc()->getNbParticles() ; ++idxPart ){
136
                if( NbParticlesEntireSystem - 1 != leaf->getSrc()->getDataDown()[idxPart]){
137
                    std::cout << "P2P probleme, should be " << NbParticlesEntireSystem - 1 <<
138
                                 " iter.data().getDataDown() "<< leaf->getSrc()->getDataDown()[idxPart] << std::endl;
139 140
                }
            }
141
        });
142 143
    }

144 145
    std::cout << "Test done..." << std::endl;

146 147 148 149 150
    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    return 0;
}