ChebyshevInterpolationFMM.cpp 10.4 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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
// ===================================================================================

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

#include <iostream>

#include <cstdio>
#include <cstdlib>
#include <string>

#include "ScalFmmConfig.h"

#include "Files/FFmaGenericLoader.hpp"

#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"

#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"

#include "Utils/FParameters.hpp"

#include "Containers/FOctree.hpp"

#ifdef _OPENMP
#include "Core/FFmmAlgorithmThread.hpp"
#else
49
#include "Core/FFmmAlgorithm.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
50 51
#endif

52
#include "Utils/FParameterNames.hpp"
53

COULAUD Olivier's avatar
COULAUD Olivier committed
54 55 56 57 58 59 60 61 62 63 64 65 66 67
/**
 * This program runs the FMM Algorithm with the Chebyshev kernel and compares the results with a direct computation.
 */
/// \file  ChebyshevInterpolationFMM.cpp
//!
//! \brief This program runs the FMM Algorithm with the interpolation kernel based on Chebyshev interpolation (1/r kernel)
//!  \authors B. Bramas, O. Coulaud
//!
//!  This code is a short example to use the Chebyshev Interpolation approach for the 1/r kernel


// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
68 69 70 71 72 73 74
    FHelpDescribeAndExit(argc, argv,
                         "Driver for Chebyshev interpolation kernel  (1/r kernel).",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
                         FParameterDefinitions::NbThreads);


75
    const std::string defaultFile(/*SCALFMMDataPath+*/"Data/unitCubeXYZQ100.bfma" );
76 77 78 79 80 81
    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
82
#ifdef _OPENMP
83 84
        omp_set_num_threads(NbThreads);
        std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
85
#else
86
        std::cout << "\n>> Sequential version.\n" << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
87
#endif
88 89 90 91 92 93 94 95 96 97 98 99 100
        //
        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 timer
        FTic time;

        // open particle file
        ////////////////////////////////////////////////////////////////////
101 102
        typedef double FReal;
        FFmaGenericLoader<FReal> loader(filename);
103 104 105 106 107 108
        //
        ////////////////////////////////////////////////////////////////////
        // begin Chebyshev kernel
        // accuracy
        const unsigned int ORDER = 7;
        // typedefs
109 110 111 112
        typedef FP2PParticleContainerIndexed<FReal>                     ContainerClass;
        typedef FSimpleLeaf<FReal, ContainerClass >                        LeafClass;
        typedef FChebCell<FReal,ORDER>                                         CellClass;
        typedef FOctree<FReal,CellClass,ContainerClass,LeafClass>  OctreeClass;
113
        //
114 115 116
        typedef FInterpMatrixKernelR<FReal>                                                                              MatrixKernelClass;
        const MatrixKernelClass MatrixKernel;
        typedef FChebSymKernel<FReal, CellClass,ContainerClass,MatrixKernelClass,ORDER>  KernelClass;
117
        //
COULAUD Olivier's avatar
COULAUD Olivier committed
118
#ifdef _OPENMP
119
        typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
120
#else
121
        typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
122
#endif
123 124 125 126 127 128 129 130 131 132
        // init oct-tree
        OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());


        { // -----------------------------------------------------
                std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
                                                                                                                                                        << " particles ..." << std::endl;
                std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
                time.tic();
                //
133
                FPoint<FReal> position;
134 135
                FReal physicalValue = 0.0;
                //
136
                for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
                        //
                        // Read particle per particle from file
                        loader.fillParticle(&position,&physicalValue);
                        //
                        // put particle in octree
                        tree.insert(position, idxPart, physicalValue);
                }

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

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

                time.tic();
                //
                KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
                //
157
	      	// false : dynamic schedule.
BRAMAS Berenger's avatar
BRAMAS Berenger committed
158
                int inUserChunckSize = 10; // To specify the chunck size in the loops (-1 is static, 0 is N/p^2, otherwise i)
159
                FmmClass algo(&tree, &kernels, inUserChunckSize);
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
                //
                algo.execute();   // Here the call of the FMM algorithm
                //
                time.tac();
                std::cout << "Timers Far Field \n"
                          << "P2M " << algo.getTime(FAlgorithmTimers::P2MTimer) << " seconds\n"
                          << "M2M " << algo.getTime(FAlgorithmTimers::M2MTimer) << " seconds\n"
                          << "M2L " << algo.getTime(FAlgorithmTimers::M2LTimer) << " seconds\n"
                          << "L2L " << algo.getTime(FAlgorithmTimers::L2LTimer) << " seconds\n"
                          << "P2P and L2P " << algo.getTime(FAlgorithmTimers::NearTimer) << " seconds\n"
                          << std::endl;


                std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << " s) ." << std::endl;
        }
        // -----------------------------------------------------
        //
        // Some output
        //
        //
        { // -----------------------------------------------------
181
                FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= loader.getNumberOfParticles() -1; ;
182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
                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();
199
                        const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
200 201
                        const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();

202
                        const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
203

204 205
                        for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                                const FSize indexPartOrig = indexes[idxPart];
206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
                                if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3)  ) {
                                        std::cout << "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] ;
                        }
                });
                std::cout <<std::endl<<"Energy: "<< energy<<std::endl;
                std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;

        }
        // -----------------------------------------------------


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