testChebAlgorithmProc.cpp 6.89 KB
Newer Older
1
// See LICENCE file at project root
2 3

// ==== CMAKE =====
4
// @FUSE_MPI
COULAUD Olivier's avatar
COULAUD Olivier committed
5
// @FUSE_BLAS
6 7 8 9 10 11 12 13
// ================

#include <iostream>

#include <cstdio>
#include <cstdlib>


14 15 16 17
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
18

19 20
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
21

22 23
#include "Utils/FParameters.hpp"
#include "Utils/FMemUtils.hpp"
24

25 26
#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"
27

28 29 30
#include "Files/FFmaGenericLoader.hpp"
#include "Files/FMpiFmaGenericLoader.hpp"
#include "Files/FMpiTreeBuilder.hpp"
31

32 33 34
#include "Core/FFmmAlgorithm.hpp"
#include "Core/FFmmAlgorithmThread.hpp"
#include "Core/FFmmAlgorithmThreadProc.hpp"
35

36
#include "Utils/FLeafBalance.hpp"
37

38
#include "Utils/FParameterNames.hpp"
39

40 41 42 43 44 45 46
/**
 * This program runs the FMM Algorithm Distributed with the Chebyshev kernel
 */

// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
47 48 49 50
    FHelpDescribeAndExit(argc, argv,
                         "Test with MPI the chebyshev FMM and compare it to the direct computation for debugging purpose.",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::NbThreads);
51

52
    typedef double FReal;
53 54
    const unsigned int ORDER = 7;
    const FReal epsilon = FReal(1e-7);
55

56 57
    typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
    typedef FSimpleLeaf<FReal, ContainerClass >  LeafClass;
58

59 60 61
    typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
    typedef FChebCell<FReal,ORDER> CellClass;
    typedef FOctree<FReal,CellClass,ContainerClass,LeafClass> OctreeClass;
62

63
    typedef FChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
64 65 66 67 68 69 70 71 72
    typedef FFmmAlgorithmThreadProc<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;

    FMpi app(argc,argv);

    const char* const filename       = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma");
    const unsigned int TreeHeight    = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5);
    const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
    const unsigned int NbThreads     = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
    FTic time;
73

74 75 76 77 78 79
    std::cout << ">> This executable has to be used to test Proc Chebyshev Algorithm. \n";

    omp_set_num_threads(NbThreads);
    std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;

    // Create Matrix Kernel
80 81
    const MatrixKernelClass MatrixKernel;

82 83
    // init particles position and physical value
    struct TestParticle{
84
        FPoint<FReal> position;
85
        FReal physicalValue;
86
        const FPoint<FReal>& getPosition(){
87 88 89
            return position;
        }
    };
90

91 92
    // open particle file
    std::cout << "Opening : " <<filename << "\n" << std::endl;
93
    FMpiFmaGenericLoader<FReal> loader(filename,app.global());
94
    if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");
95

96 97 98
    OctreeClass tree(TreeHeight, SubTreeHeight,loader.getBoxWidth(),loader.getCenterOfBox());

    time.tic();
99 100
    TestParticle* particles = new TestParticle[loader.getMyNumberOfParticles()];
    memset(particles,0,(unsigned int) (sizeof(TestParticle)* loader.getMyNumberOfParticles()));
101
    for(FSize idxPart = 0 ; idxPart < loader.getMyNumberOfParticles() ; ++idxPart){
102
        loader.fillParticle(&particles[idxPart].position,&particles[idxPart].physicalValue);
103 104
    }

105 106
    FVector<TestParticle> finalParticles;
    FLeafBalance balancer;
107
    FMpiTreeBuilder< FReal,TestParticle >::DistributeArrayToContainer(app.global(),particles,
108 109 110 111 112
                                                                loader.getMyNumberOfParticles(),
                                                                tree.getBoxCenter(),
                                                                tree.getBoxWidth(),tree.getHeight(),
                                                                &finalParticles, &balancer);
    { // -----------------------------------------------------
113 114
        std::cout << "Creating & Inserting " << loader.getMyNumberOfParticles() << " particles ..." << std::endl;
        std::cout << "For a total of " << loader.getNumberOfParticles() << " particles ..." << std::endl;
115 116 117
        std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
        time.tic();

118
        for(FSize idxPart = 0 ; idxPart < finalParticles.getSize() ; ++idxPart){
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
            // put in tree
            tree.insert(finalParticles[idxPart].position, idxPart, finalParticles[idxPart].physicalValue);
        }

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

    { // -----------------------------------------------------
        std::cout << "\nChebyshev FMM (ORDER="<< ORDER << ",EPS="<< epsilon <<") ... " << std::endl;
        time.tic();
        KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(), &MatrixKernel, epsilon);
        FmmClass algorithm(app.global(),&tree, &kernels);
        algorithm.execute();
        time.tac();
        std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
    } // -----------------------------------------------------


    { // -----------------------------------------------------
        std::cout << "\nError computation ... " << std::endl;
141 142
        FReal potential{};
        FReal fx{}, fy{}, fz{};
143 144 145 146 147 148 149
        { // Check that each particle has been summed with all other

            tree.forEachLeaf([&](LeafClass* leaf){
                const FReal*const potentials = leaf->getTargets()->getPotentials();
                const FReal*const forcesX = leaf->getTargets()->getForcesX();
                const FReal*const forcesY = leaf->getTargets()->getForcesY();
                const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
150
                const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
151

152
                for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169
                    potential += potentials[idxPart];
                    fx += forcesX[idxPart];
                    fy += forcesY[idxPart];
                    fz += forcesZ[idxPart];
                }
            });
        }

        // Print for information
        std::cout << "Potential " << potential << std::endl;
        std::cout << "Fx " << fx << std::endl;
        std::cout << "Fy " << fy << std::endl;
        std::cout << "Fz " << fz << std::endl;

    } // end Chebyshev kernel

    return 0;
170
}