RotationMPIFMM.cpp 10.7 KB
Newer Older
COULAUD Olivier's avatar
COULAUD Olivier committed
1
// ===================================================================================
2 3 4 5
// Copyright ScalFmm 2016 INRIA, Olivier Coulaud, Bérenger Bramas,
// Matthias Messner olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the
// FMM.
6
//
7
// This software is governed by the CeCILL-C and LGPL licenses and
COULAUD Olivier's avatar
COULAUD Olivier committed
8
// abiding by the rules of distribution of free software.
9 10 11
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
COULAUD Olivier's avatar
COULAUD Olivier committed
12 13 14 15
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 17 18
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
COULAUD Olivier's avatar
COULAUD Olivier committed
19 20 21 22 23 24 25 26 27 28 29 30 31
// ===================================================================================

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

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


#include "ScalFmmConfig.h"
32 33 34
#include "Containers/FOctree.hpp"
#include "Utils/FMpi.hpp"
#include "Core/FFmmAlgorithmThreadProc.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
35

36 37 38
#include "Files/FFmaGenericLoader.hpp"
#include "Files/FMpiFmaGenericLoader.hpp"
#include "Files/FMpiTreeBuilder.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
39

40
#include "BalanceTree/FLeafBalance.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
41

42 43
#include "Kernels/Rotation/FRotationKernel.hpp"
#include "Kernels/Rotation/FRotationCell.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
44

45 46
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
47

48
#include "Utils/FParameters.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
49

50
#include "Utils/FParameterNames.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

/// \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
//!
66
//!     \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
67 68 69
//!
//

70

COULAUD Olivier's avatar
COULAUD Olivier committed
71 72 73 74

// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
75 76 77 78 79 80 81 82 83 84 85 86 87 88
    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
89
#ifdef _OPENMP
90 91
    omp_set_num_threads(NbThreads);
    std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
92
#else
93
    std::cout << "\n>> Sequential version.\n" << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
94
#endif
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
95
    //
96 97 98 99 100 101 102 103
    std::cout <<	 "Parameters  "<< std::endl
                  <<     "      Octree Depth      "<< TreeHeight <<std::endl
                      <<	 "      SubOctree depth "<< SubTreeHeight <<std::endl
                          <<     "      Input file  name: " <<filename <<std::endl
                              <<     "      Thread number:  " << NbThreads <<std::endl
                                  <<std::endl;
    //init values for MPI
    FMpi app(argc,argv);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
104
    //
105 106 107
    // init timer
    FTic time;

108 109
    typedef double FReal;
    FMpiFmaGenericLoader<FReal>* loader = new FMpiFmaGenericLoader<FReal>(filename,app.global());
110 111 112 113 114 115 116 117 118 119

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



    // begin spherical kernel
    // accuracy
    const unsigned int P = 22;
    // typedefs
120 121 122 123
    typedef FP2PParticleContainerIndexed<FReal>                     ContainerClass;
    typedef FSimpleLeaf< FReal, ContainerClass >                       LeafClass;
    typedef FRotationCell<FReal, P>                                             CellClass;
    typedef FOctree<FReal, CellClass,ContainerClass,LeafClass>  OctreeClass;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
124
    //
125
    typedef FRotationKernel<FReal, CellClass, ContainerClass , P>   KernelClass;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
126
    //
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
    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{
143
            FSize indexInFile;
144
            FPoint<FReal> position;
145
            FReal physicalValue;
146
            const FPoint<FReal>& getPosition(){
147 148 149 150
                return position;
            }
        };
        TestParticle* particles = new TestParticle[loader->getMyNumberOfParticles()];
151
        memset(particles, 0, (sizeof(TestParticle) * loader->getMyNumberOfParticles()));
152 153

        //idx (in file) of the first part that will be used by this proc.
154
        FSize idxStart = loader->getStart();
155
        
156
        for(FSize idxPart = 0 ; idxPart < loader->getMyNumberOfParticles() ; ++idxPart){
157 158 159 160 161 162 163 164 165
            //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;
166
        FMpiTreeBuilder<FReal, TestParticle >::DistributeArrayToContainer(app.global(), particles, loader->getMyNumberOfParticles(),
167 168 169 170
                                                                    tree.getBoxCenter(),
                                                                    tree.getBoxWidth(),
                                                                    tree.getHeight(), &finalParticles,&balancer);

171 172
        std::cout << "Local nb particles after sort "
                  << finalParticles.getSize() << std::endl;
173
        for(FSize idx = 0 ; idx < finalParticles.getSize(); ++idx){
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
            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;
    }
    // -----------------------------------------------------
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
201
    //
202
    // Some output
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
203 204
    //
    //
205 206

    { // -----------------------------------------------------
207
        FSize N1=0, N2= loader->getNumberOfParticles()/2, N3= (loader->getNumberOfParticles()-1); ;
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
        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();
224
            const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
225 226
            const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();

227
            const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
228

229 230
            for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                const FSize indexPartOrig = indexes[idxPart];
231 232 233 234 235 236 237 238 239 240 241 242 243 244
                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;
        }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
245
    }
246 247 248 249
    delete loader;
    // -----------------------------------------------------

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