ChebyshevInterpolationMPIFMM.cpp 11.1 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 14

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

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


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

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

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

25 26 27
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
#include "Kernels/Chebyshev/FChebCell.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
28

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

32 33
#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
34

COULAUD Olivier's avatar
COULAUD Olivier committed
35 36 37
#ifdef SCALFMM_USE_EZTRACE
#include "eztrace.h"
#endif
38
/// \file
COULAUD Olivier's avatar
COULAUD Olivier committed
39 40 41 42 43 44 45 46 47 48
//!
//! \brief This program runs the MPI FMM with Chebyshev interpolation of 1/r kernel
//!  \authors B. Bramas, O. Coulaud
//!
//!  This code is a short example to use the FMM Algorithm Proc with Chebyshev Interpolation for the 1/r kernel


// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
49
    ///////// PARAMETERS HANDLING //////////////////////////////////////
50
    FHelpDescribeAndExit(argc, argv,
51 52 53 54 55
                         "Driver for Chebyshev Interpolation kernel using MPI  (1/r kernel). "
                         "Usully run using : mpirun -np nb_proc_needed ./ChebyshevInterpolationAlgorithm [params].",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
                         FParameterDefinitions::NbThreads);
56

57 58 59 60 61
    const std::string defaultFile("../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);
62

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

    //
67
    std::cout << "Parameters"<< std::endl
68 69 70 71
              << "      Octree Depth      " << TreeHeight    << std::endl
              << "      SubOctree depth   " << SubTreeHeight << std::endl
              << "      Input file  name: " << filename      << std::endl
              << "      Thread count :    " << NbThreads     << std::endl
72 73 74 75 76 77
              << std::endl;


    ///////// VAR INIT /////////////////////////////////////////////////

    // Initialize values for MPI
COULAUD Olivier's avatar
COULAUD Olivier committed
78 79 80
#ifdef SCALFMM_USE_EZTRACE
   eztrace_start();
#endif
81
    FMpi app(argc,argv);
COULAUD Olivier's avatar
COULAUD Olivier committed
82 83
#ifdef SCALFMM_USE_EZTRACE
   eztrace_pause();
84
#endif
85
    //
86
    // Initialize timer
87
    FTic time;
COULAUD Olivier's avatar
COULAUD Olivier committed
88

89
    typedef double FReal;
COULAUD Olivier's avatar
COULAUD Olivier committed
90

91 92
    // Creation of the particle loader
    FMpiFmaGenericLoader<FReal> loader(filename,app.global());
93
    if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!") ;
COULAUD Olivier's avatar
COULAUD Olivier committed
94 95


96 97
    // Begin spherical kernel
    // Accuracy
98
    const unsigned int ORDER = 7;
99 100 101 102 103 104
    // Typedefs
    using ContainerClass = FP2PParticleContainerIndexed<FReal>;
    using LeafClass      = FSimpleLeaf<FReal, ContainerClass>;
    using CellClass      = FChebCell<FReal,ORDER>;
    using OctreeClass    = FOctree<FReal,CellClass,ContainerClass,LeafClass>;

105
    using MatrixKernelClass = FInterpMatrixKernelR<FReal>;
106
    const MatrixKernelClass MatrixKernel;
COULAUD Olivier's avatar
COULAUD Olivier committed
107

108 109
    using KernelClass    = FChebSymKernel<FReal, CellClass,ContainerClass,MatrixKernelClass,ORDER>;
    using FmmClassProc   = FFmmAlgorithmThreadProc<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass>;
110

111
    // Initialize empty oct-tree
112 113 114 115 116 117 118 119 120 121 122
    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();

123
        /* Mock particle structure to balance the tree over the processes. */
124
        struct TestParticle{
125 126 127 128
            FSize index;             // Index of the particle in the original file.
            FPoint<FReal> position;  // Spatial position of the particle.
            FReal physicalValue;     // Physical value of the particle.
            /* Returns the particle position. */
129 130 131 132 133
            const FPoint<FReal>& getPosition(){
                return position;
            }
        };

134
        // Temporary array of particles read by this process.
135
        TestParticle* particles = new TestParticle[loader.getMyNumberOfParticles()];
136
        memset(particles, 0, (sizeof(TestParticle) * loader.getMyNumberOfParticles()));
137

138
        // Index (in file) of the first particle that will be read by this process.
139
        FSize idxStart = loader.getStart();
140
        std::cout << "Proc:" << app.global().processId() << " start-index: " << idxStart << std::endl;
141

142
        // Read particles from parts.
143
        for(FSize idxPart = 0 ; idxPart < loader.getMyNumberOfParticles() ; ++idxPart){
144
            // Store the index (in the original file) the particle.
145
            particles[idxPart].index = idxPart + idxStart;
146 147 148
            // Read particle from file
            loader.fillParticle(&particles[idxPart].position,
                                &particles[idxPart].physicalValue);
149 150
        }

151
        // Final vector of particles
152 153 154
        FVector<TestParticle> finalParticles;
        FLeafBalance balancer;

155 156 157 158 159 160 161 162 163 164
        // Redistribute particules between processes
        FMpiTreeBuilder< FReal, TestParticle >::
            DistributeArrayToContainer(app.global(),
                                       particles,
                                       loader.getMyNumberOfParticles(),
                                       tree.getBoxCenter(),
                                       tree.getBoxWidth(),
                                       tree.getHeight(),
                                       &finalParticles,
                                       &balancer);
165

166 167 168 169
        // Free temporary array memory.
        delete[] particles;

        // Insert final particles into tree.
170
        for(FSize idx = 0 ; idx < finalParticles.getSize(); ++idx){
171 172 173
            tree.insert(finalParticles[idx].position,
                        finalParticles[idx].index,
                        finalParticles[idx].physicalValue);
174 175 176 177 178
        }

        time.tac();
        double timeUsed = time.elapsed();
        double minTime,maxTime;
179 180
        std::cout << "Proc:" << app.global().processId()
                  << " "     << finalParticles.getSize()
181
                  << "particles have been inserted in the tree. (@Reading and Inserting Particles = "
182 183 184
                  << time.elapsed() << " s)."
                  << std::endl;

185 186 187
        MPI_Reduce(&timeUsed,&minTime,1,MPI_DOUBLE,MPI_MIN,0,app.global().getComm());
        MPI_Reduce(&timeUsed,&maxTime,1,MPI_DOUBLE,MPI_MAX,0,app.global().getComm());
        if(app.global().processId() == 0){
188
            std::cout << "readinsert-time-min:" << minTime
189
                      << " readinsert-time-max:" << maxTime
190
                      << std::endl;
191 192 193 194 195 196 197
        }
    } // -----------------------------------------------------

    { // -----------------------------------------------------
        std::cout << "\nChebyshev Interpolation  FMM Proc (P="<< ORDER << ") ... " << std::endl;

        time.tic();
198 199

        // Kernels to use (pointer because of the limited size of the stack)
200
        KernelClass *kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
201
        // MPI FMM algorithm
202
        FmmClassProc algorithm(app.global(),&tree, kernels);
203 204 205
        // FMM exectution
        algorithm.execute();

206 207 208 209 210 211 212
        time.tac();
        double timeUsed = time.elapsed();
        double minTime,maxTime;
        std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << " s)." << std::endl;
        MPI_Reduce(&timeUsed,&minTime,1,MPI_DOUBLE,MPI_MIN,0,app.global().getComm());
        MPI_Reduce(&timeUsed,&maxTime,1,MPI_DOUBLE,MPI_MAX,0,app.global().getComm());
        if(app.global().processId() == 0){
213
            std::cout << "exec-time-min:" << minTime
214
                      << " exec-time-max:" << maxTime
215
                      << std::endl;
216 217
        }

218
        // Free kernels from memory
219 220 221 222 223 224 225 226
        delete kernels;
    }
    // -----------------------------------------------------
    //
    // Some output
    //
    //
    { // -----------------------------------------------------
227
        FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= (loader.getNumberOfParticles()-1); ;
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244
        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();
245
            const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
246 247
            const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();

248
            const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
249

250 251
            for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                const FSize indexPartOrig = indexes[idxPart];
252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268
                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;
        }
    }
    // -----------------------------------------------------

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