RotationMPIFMM.cpp 9.7 KB
Newer Older
1
// See LICENCE file at project root
COULAUD Olivier's avatar
COULAUD Olivier committed
2 3 4 5 6 7 8 9 10 11 12 13

// ==== CMAKE =====
// @FUSE_MPI
// ================

#include <iostream>
#include <stdexcept>
#include <cstdio>
#include <cstdlib>


#include "ScalFmmConfig.h"
14 15 16
#include "Containers/FOctree.hpp"
#include "Utils/FMpi.hpp"
#include "Core/FFmmAlgorithmThreadProc.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
17

18 19 20
#include "Files/FFmaGenericLoader.hpp"
#include "Files/FMpiFmaGenericLoader.hpp"
#include "Files/FMpiTreeBuilder.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
21

22
#include "Utils/FLeafBalance.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
23

24 25
#include "Kernels/Rotation/FRotationKernel.hpp"
#include "Kernels/Rotation/FRotationCell.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
26

27 28
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
29

30
#include "Utils/FParameters.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
31

32
#include "Utils/FParameterNames.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

/// \file  RotationFMMProc.cpp
//!
//! \brief This program runs the MPI FMM  with harmonic spherical approximation of 1/r kernel
//!  \authors B. Bramas, O. Coulaud
//!
//!  This code is a short example to use the rotation harmonic spherical approximation for the 1/r kernel
//!
//!
//!  <b> General arguments:</b>
//!     \param   -help(-h)      to see the parameters available in this driver
//!     \param   -depth          The depth of the octree
//!     \param   -subdepth     Specifies the size of the sub octree
//!     \param   -t                   The number of threads
//!
48
//!     \param   -f name          Name of the particles file with extension (.fma or .bfma). The data in  file have to be in our FMA format
COULAUD Olivier's avatar
COULAUD Olivier committed
49 50 51
//!
//

52

COULAUD Olivier's avatar
COULAUD Olivier committed
53 54 55 56

// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
57 58 59 60 61 62 63 64 65 66 67 68 69 70
    FHelpDescribeAndExit(argc, argv,
                         "Driver for Rotation Spherical kernel using MPI  (1/r kernel).\n"
                         "CMD >> mpirun -np nb_proc_needed ./RotationFMMProc [params]",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
                         FParameterDefinitions::NbThreads);

    const std::string defaultFile(/*SCALFMMDataPath+*/"../Data/test20k.fma");
    const std::string filename                = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFile.c_str());
    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);


COULAUD Olivier's avatar
COULAUD Olivier committed
71
#ifdef _OPENMP
72 73
    omp_set_num_threads(NbThreads);
    std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
74
#else
75
    std::cout << "\n>> Sequential version.\n" << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
76
#endif
77
    //
78
    std::cout <<         "Parameters  "<< std::endl
79
                  <<     "      Octree Depth      "<< TreeHeight <<std::endl
80
                      <<         "      SubOctree depth "<< SubTreeHeight <<std::endl
81 82 83 84 85
                          <<     "      Input file  name: " <<filename <<std::endl
                              <<     "      Thread number:  " << NbThreads <<std::endl
                                  <<std::endl;
    //init values for MPI
    FMpi app(argc,argv);
86
    //
87 88 89
    // init timer
    FTic time;

90 91
    typedef double FReal;
    FMpiFmaGenericLoader<FReal>* loader = new FMpiFmaGenericLoader<FReal>(filename,app.global());
92 93 94 95 96 97 98 99 100 101

    if(!loader->isOpen()) throw std::runtime_error("Particle file couldn't be opened!") ;
    ////////////////////////////////////////////////////////////////////



    // begin spherical kernel
    // accuracy
    const unsigned int P = 22;
    // typedefs
102 103 104 105
    typedef FP2PParticleContainerIndexed<FReal>                     ContainerClass;
    typedef FSimpleLeaf< FReal, ContainerClass >                       LeafClass;
    typedef FRotationCell<FReal, P>                                             CellClass;
    typedef FOctree<FReal, CellClass,ContainerClass,LeafClass>  OctreeClass;
106
    //
107
    typedef FRotationKernel<FReal, CellClass, ContainerClass , P>   KernelClass;
108
    //
109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
    typedef FFmmAlgorithmThreadProc<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClassProc;

    // init oct-tree
    OctreeClass tree(TreeHeight, SubTreeHeight, loader->getBoxWidth(), loader->getCenterOfBox());


    { // -----------------------------------------------------
        if(app.global().processId() == 0){
            std::cout << "Creating & Inserting " << loader->getNumberOfParticles()
                      << " particles ..." << std::endl;
            std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
        }
        time.tic();
        //

        struct TestParticle{
125
            FSize indexInFile;
126
            FPoint<FReal> position;
127
            FReal physicalValue;
128
            const FPoint<FReal>& getPosition(){
129 130 131 132
                return position;
            }
        };
        TestParticle* particles = new TestParticle[loader->getMyNumberOfParticles()];
133
        memset(particles, 0, (sizeof(TestParticle) * loader->getMyNumberOfParticles()));
134 135

        //idx (in file) of the first part that will be used by this proc.
136
        FSize idxStart = loader->getStart();
137

138
        for(FSize idxPart = 0 ; idxPart < loader->getMyNumberOfParticles() ; ++idxPart){
139 140 141 142 143 144 145 146 147
            //Storage of the index (in the original file) of each part.
            particles[idxPart].indexInFile = idxPart + idxStart;

            // Read particles from file
            loader->fillParticle(&particles[idxPart].position,&particles[idxPart].physicalValue);
        }

        FVector<TestParticle> finalParticles;
        FLeafBalance balancer;
148
        FMpiTreeBuilder<FReal, TestParticle >::DistributeArrayToContainer(app.global(), particles, loader->getMyNumberOfParticles(),
149 150 151 152
                                                                    tree.getBoxCenter(),
                                                                    tree.getBoxWidth(),
                                                                    tree.getHeight(), &finalParticles,&balancer);

153 154
        std::cout << "Local nb particles after sort "
                  << finalParticles.getSize() << std::endl;
155
        for(FSize idx = 0 ; idx < finalParticles.getSize(); ++idx){
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
            tree.insert(finalParticles[idx].position,finalParticles[idx].indexInFile,finalParticles[idx].physicalValue);
        }

        delete[] particles;

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

    { // -----------------------------------------------------
        std::cout << "\nRotation harmonic Spherical FMM Proc (P="<< P << ") ... " << std::endl;

        time.tic();
        //
        // Here we use a pointer due to the limited size of the stack
        //
        KernelClass *kernels = new KernelClass(TreeHeight, loader->getBoxWidth(), loader->getCenterOfBox());
        //
        FmmClassProc algorithm(app.global(),&tree, kernels);
        //
        algorithm.execute();   // Here the call of the FMM algorithm
        //
        time.tac();
        std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
    }
    // -----------------------------------------------------
183
    //
184
    // Some output
185 186
    //
    //
187 188

    { // -----------------------------------------------------
189
        FSize N1=0, N2= loader->getNumberOfParticles()/2, N3= (loader->getNumberOfParticles()-1); ;
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
        FReal energy =0.0 ;
        //
        //   Loop over all leaves
        //
        std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
        std::cout << std::scientific;
        std::cout.precision(10) ;

        tree.forEachLeaf([&](LeafClass* leaf){
            const FReal*const posX = leaf->getTargets()->getPositions()[0];
            const FReal*const posY = leaf->getTargets()->getPositions()[1];
            const FReal*const posZ = leaf->getTargets()->getPositions()[2];
            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();
206
            const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
207 208
            const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();

209
            const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
210

211 212
            for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                const FSize indexPartOrig = indexes[idxPart];
213 214 215 216 217 218 219 220 221 222 223 224 225 226
                if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3)  ) {
                    std::cout << "Proc "<< app.global().processId() << " Index "<< indexPartOrig << "  potential  " << potentials[idxPart]
                                 << " Pos "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
                                    << "   Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
                }
                energy += potentials[idxPart]*physicalValues[idxPart] ;
            }
        });

        FReal gloEnergy = app.global().reduceSum(energy);
        if(0 == app.global().processId()){
            std::cout <<std::endl << "Proc "<< app.global().processId() << " Energy: "<< gloEnergy <<std::endl;
            std::cout <<std::endl <<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
        }
227
    }
228 229 230 231
    delete loader;
    // -----------------------------------------------------

    return 0;
COULAUD Olivier's avatar
COULAUD Olivier committed
232
}