compareAllPoissonKernels.cpp 40 KB
Newer Older
1
// See LICENCE file at project root
2 3 4 5 6
//
// ==== CMAKE =====
//
// ================

7
// Keep in private GIT
8
//
9

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

#include <iostream>
#include <stdexcept>
#include <string>
#include <cstdlib>
#include <cstdio>
//
#include "ScalFmmConfig.h"
#include "Utils/FTic.hpp"
#include "Utils/FParameters.hpp"

#include "Files/FFmaGenericLoader.hpp"

#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"

#include "Core/FFmmAlgorithmThread.hpp"

28
#ifdef SCALFMM_USE_BLAS
29 30 31 32 33 34 35 36 37 38
// chebyshev kernel

#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
#endif
//
// spherical kernel
#include "Kernels/Spherical/FSphericalCell.hpp"
39
#ifdef SCALFMM_USE_BLAS
40 41 42 43 44 45 46 47 48 49 50 51
#include "Kernels/Spherical/FSphericalBlasKernel.hpp"
#include "Kernels/Spherical/FSphericalBlockBlasKernel.hpp"
#endif
//
// taylor kernel
#include "Kernels/Taylor/FTaylorCell.hpp"
#include "Kernels/Taylor/FTaylorKernel.hpp"
//
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"

//Classical Spherical kernel
52 53
#include "Kernels/Spherical/FSphericalCell.hpp"
#include "Kernels/Spherical/FSphericalKernel.hpp"
54 55 56 57 58

//Rotation kernel
#include "Kernels/Rotation/FRotationKernel.hpp"
#include "Kernels/Rotation/FRotationCell.hpp"

59
#ifdef SCALFMM_USE_FFT
60 61 62 63 64 65 66
// Uniform grid kernel
#include "Kernels/Uniform/FUnifCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"

#include "Kernels/Uniform/FUnifKernel.hpp"
#endif

67
#include "Utils/FParameterNames.hpp"
68 69 70 71 72 73 74 75 76 77 78 79

/**
 * This program compares two different kernels, eg., the Chebyshev kernel with
 * the SphericalBlas kernel.
 */
//
/// \file  compareAllPoissonKernels.cpp
//!
//! \brief compareAllPoissonKernels: Driver to compare all different implementations for 1/r kernel.
//!
//! \details  This driver  compare all different implementations provided by our library for the classical Poisson kernel (1/r)<br>
//!     We check the following kernels <br>
80 81
//!      - Spherical expansion; Spherical with rotation optimization
//!      - Taylor expansion; Spherical with rotation optimization
82
//!      - if BLAS is activated
83 84
//!               - Blas And Block Blas optiization of the M2L operator
//!               - Chebychev and symetric Chebychev interpolation
85 86
//!      - if FFT is activated:  interpolation on uniform grid
//!<br>
87

88 89 90 91

// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
92 93 94 95 96
    FHelpDescribeAndExit(argc, argv,
                         "Driver for testing different approximations  for the  1/r kernel.",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
                         FParameterDefinitions::NbThreads, FParameterDefinitions::SHDevelopment);
97

98
    typedef double FReal;
99
        // get info from commande line
100 101 102 103 104
    const std::string  filename(FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/UTest/unitCubeRef20kDouble.bfma"));
    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, omp_get_max_threads());
    const int DevP                   = FParameters::getValue(argc, argv, FParameterDefinitions::SHDevelopment.options, 11);
105

106
        //
107
#ifdef _OPENMP
108
        omp_set_num_threads(NbThreads);
109
#else
110
    std::cout << "\n>> Sequential version.\n" << std::endl;
111 112
#endif

113 114 115 116 117 118
        std::cout <<     "Parameters  "<< std::endl
                        <<     "      Octree Depth      \t"<< TreeHeight <<std::endl
                        <<        "      SubOctree depth \t"<< SubTreeHeight <<std::endl
                        <<     "      Input file  name: \t" <<filename <<std::endl
                        <<     "      Thread number:  \t" << NbThreads <<std::endl
                        <<std::endl;
119

120 121
        // init timer
        FTic time;
122

123
    FFmaGenericLoader<FReal> loader(filename);
124 125
        //  if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");
        //
126

127
        FSize nbParticles = loader.getNumberOfParticles() ;
128
    FmaRWParticle<FReal, 8,8>* const particles = new FmaRWParticle<FReal, 8,8>[nbParticles];
129 130 131 132 133 134
        //
        loader.fillParticle(particles,nbParticles);
        //
        ////////////////////////////////////////////////////////////////////
        //  Compute direct energy
        FReal energyD =0.0, totPhysicalValue =0.0;
135 136

#pragma omp parallel for reduction(+:energyD,totPhysicalValue)
137
    for(FSize idx = 0 ; idx < loader.getNumberOfParticles()  ; ++idx){
138 139 140 141 142 143
                energyD             +=  particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
                totPhysicalValue += particles[idx].getPhysicalValue() ;
        }
        std::cout << " Total Physical values: "<< totPhysicalValue <<std::endl;
        std::cout << " Energy of the system: "<< energyD <<std::endl;
        ////////////////////////////////////////////////////////////////////
144

145
#ifdef  SCALFMM_USE_BLAS
146
        {	// begin Chebyshev kernel
147

148 149 150
                // accuracy
                const unsigned int ORDER = 7;
                std::cout << "\nFChebKernel FMM ... ORDER: " << ORDER <<std::endl;
151

152
                // typedefs
153 154 155 156 157
        typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
        typedef FSimpleLeaf<FReal, ContainerClass> LeafClass;
        typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
        typedef FChebCell<FReal, ORDER> CellClass;
        typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
158

159
        typedef FChebKernel<FReal, CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
160
                typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
161 162


163 164
                // init oct-tree
                OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
165 166 167 168

    // Create Matrix Kernel
    const MatrixKernelClass MatrixKernel;

169 170
                { // -----------------------------------------------------
                        time.tic();
171
            for(FSize idxPart = 0 ; idxPart < nbParticles; ++idxPart){
172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187
                                // put in tree
                                tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
                        }
                        time.tac();
                        std::cout << "(FChebKernel @ Inserting Particles = "<< time.elapsed() << " s)." << std::endl;
                } // -----------------------------------------------------

                { // -----------------------------------------------------
                        time.tic();
                        KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
                        FmmClass algorithm(&tree, &kernels);
                        algorithm.execute();
                        time.tac();
                        std::cout <<"(FChebKernel @Algorithm = " << time.elapsed() << " s)." << std::endl;
                } // -----------------------------------------------------
                FReal energy = 0.0;
188 189
        FMath::FAccurater<FReal> potentialDiff;
        FMath::FAccurater<FReal> fx, fy, fz;
190 191 192 193 194 195 196 197
                { // Check that each particle has been summed with all other

                        tree.forEachLeaf([&](LeafClass* leaf){
                                const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
                                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();
198 199
                const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
                const FVector<FSize>& indexes       = leaf->getTargets()->getIndexes();
200

201 202
                for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const FSize indexPartOrig = indexes[idxPart];
203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
                                        potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
                                        fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
                                        fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
                                        fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
                                        energy += potentials[idxPart]*physicalValues[idxPart] ;
                                }
                        });
                }

                // Print for information
                std::cout << "FChebKernel Energy "  << FMath::Abs(energy-energyD) /energyD << std::endl;
                std::cout << "FChebKernel Potential " << potentialDiff << std::endl;
                std::cout << "FChebKernel Fx " << fx << std::endl;
                std::cout << "FChebKernel Fy " << fy << std::endl;
                std::cout << "FChebKernel Fz " << fz << std::endl;

        } // end Chebyshev kernel
        {	// begin ChebSymKernel kernel

                // accuracy
                const unsigned int ORDER = 7;
                std::cout << "\nFChebSymKernel FMM ... ORDER: " << ORDER <<std::endl;

                // typedefs
227 228 229 230 231
        typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
        typedef FSimpleLeaf<FReal, ContainerClass> LeafClass;
        typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
        typedef FChebCell<FReal, ORDER> CellClass;
        typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
232

233
        typedef FChebSymKernel<FReal, CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
234
                typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
235 236


237 238
                // init oct-tree
                OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
239 240 241 242

    // Create Matrix Kernel
    const MatrixKernelClass MatrixKernel;

243 244
                { // -----------------------------------------------------
                        time.tic();
245

246
            for(FSize idxPart = 0 ; idxPart < nbParticles; ++idxPart){
247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263
                                // put in tree
                                tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
                        }

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

                { // -----------------------------------------------------
                        time.tic();
                        KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
                        FmmClass algorithm(&tree, &kernels);
                        algorithm.execute();
                        time.tac();
                        std::cout << "(FChebSymKernel @Algorithm = " << time.elapsed() << " s)." << std::endl;
                } // -----------------------------------------------------
                FReal energy = 0.0;
264 265
        FMath::FAccurater<FReal> potentialDiff;
        FMath::FAccurater<FReal> fx, fy, fz;
266 267 268 269 270 271 272 273
                { // Check that each particle has been summed with all other

                        tree.forEachLeaf([&](LeafClass* leaf){
                                const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
                                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();
274 275
                const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
                const FVector<FSize>& indexes       = leaf->getTargets()->getIndexes();
276

277 278
                for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const FSize indexPartOrig = indexes[idxPart];
279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302
                                        potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
                                        fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
                                        fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
                                        fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
                                        energy += potentials[idxPart]*physicalValues[idxPart] ;
                                }
                        });
                }

                // Print for information
                std::cout << "FChebSymKernel Energy "  << FMath::Abs(energy-energyD) /energyD << std::endl;
                std::cout << "FChebSymKernel Potential " << potentialDiff << std::endl;
                std::cout << "FChebSymKernel Fx " << fx << std::endl;
                std::cout << "FChebSymKernel Fy " << fy << std::endl;
                std::cout << "FChebSymKernel Fz " << fz << std::endl;

        } // end Chebyshev kernel
        //
        ////////////////////////////////////////////////////////////////////
        //
        {	// begin FFmaBlas kernel FSphericalBlockBlasKernel
                        std::cout << "\nFFmaBlas FMM ... P: " <<DevP << std::endl;

                        // typedefs
303 304 305 306 307
            typedef FSphericalCell<FReal>                 CellClass;
            typedef FP2PParticleContainerIndexed<FReal>         ContainerClass;
            typedef FSimpleLeaf<FReal,  ContainerClass >                     LeafClass;
            typedef FOctree<FReal,  CellClass, ContainerClass , LeafClass >  OctreeClass;
            typedef FSphericalBlasKernel<FReal,  CellClass, ContainerClass > KernelClass;
308
                        typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
309

310 311 312
                        // init cell class and oct-tree
                        CellClass::Init(DevP, true); // only for blas
                        OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
313

314 315
                        { // -----------------------------------------------------
                                time.tic();
316

317
                for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335
                                        // put in tree
                                        tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
                                }

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

                        // -----------------------------------------------------
                        time.tic();
                        KernelClass kernels(DevP, TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
                        FmmClass algorithm(&tree, &kernels);
                        algorithm.execute();
                        time.tac();
                        std::cout << "(FFmaBlas @Algorithm = " << time.elapsed() << " s)." << std::endl;
                        // -----------------------------------------------------

                        FReal energy = 0.0;
336 337
            FMath::FAccurater<FReal> potentialDiff;
            FMath::FAccurater<FReal> fx, fy, fz;
338 339 340 341 342 343 344 345
                        { // Check that each particle has been summed with all other

                                tree.forEachLeaf([&](LeafClass* leaf){
                                        const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
                                        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();
346 347
                    const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
                    const FVector<FSize>& indexes       = leaf->getTargets()->getIndexes();
348

349 350
                    for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                        const FSize indexPartOrig = indexes[idxPart];
351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373
                                                potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
                                                fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
                                                fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
                                                fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
                                                energy += potentials[idxPart]*physicalValues[idxPart] ;
                                        }
                                });
                        }

                        // Print for information
                        std::cout << "FFmaBlas Energy "  << FMath::Abs(energy-energyD) /energyD << std::endl;
                        std::cout << "FFmaBlas Potential " << potentialDiff << std::endl;
                        std::cout << "FFmaBlas Fx " << fx << std::endl;
                        std::cout << "FFmaBlas Fy " << fy << std::endl;
                        std::cout << "FFmaBlas Fz " << fz << std::endl;
                } // end FFmaBlas kernel
        //
        ////////////////////////////////////////////////////////////////////
        //
        {	// begin FFmaBlockBlas kernel
                        std::cout << "\nFFmaBlockBlas FMM ... P: " <<DevP << std::endl;

                        // typedefs
374 375 376 377 378
            typedef FSphericalCell<FReal>                 CellClass;
            typedef FP2PParticleContainerIndexed<FReal>         ContainerClass;
            typedef FSimpleLeaf<FReal,  ContainerClass >                     LeafClass;
            typedef FOctree<FReal,  CellClass, ContainerClass , LeafClass >  OctreeClass;
            typedef FSphericalBlockBlasKernel<FReal,  CellClass, ContainerClass > KernelClass;
379
                        typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
380

381 382 383
                        // init cell class and oct-tree
                        CellClass::Init(DevP, true); // only for blas
                        OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
384

385 386
                        { // -----------------------------------------------------
                                time.tic();
387

388
                for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406
                                        // put in tree
                                        tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
                                }

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

                        // -----------------------------------------------------
                        time.tic();
                        KernelClass kernels(DevP, TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
                        FmmClass algorithm(&tree, &kernels);
                        algorithm.execute();
                        time.tac();
                        std::cout << "(FFmaBlockBlas @Algorithm = " << time.elapsed() << " s)." << std::endl;
                        // -----------------------------------------------------

                        FReal energy = 0.0;
407 408
            FMath::FAccurater<FReal> potentialDiff;
            FMath::FAccurater<FReal> fx, fy, fz;
409 410 411 412 413 414 415 416
                        { // Check that each particle has been summed with all other

                                tree.forEachLeaf([&](LeafClass* leaf){
                                        const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
                                        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();
417 418
                    const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
                    const FVector<FSize>& indexes       = leaf->getTargets()->getIndexes();
419

420 421
                    for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                        const FSize indexPartOrig = indexes[idxPart];
422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437
                                                potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
                                                fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
                                                fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
                                                fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
                                                energy += potentials[idxPart]*physicalValues[idxPart] ;
                                        }
                                });
                        }

                        // Print for information
                        std::cout << "FFmaBlockBlas Energy "  << FMath::Abs(energy-energyD) /energyD << std::endl;
                        std::cout << "FFmaBlockBlas Potential " << potentialDiff << std::endl;
                        std::cout << "FFmaBlockBlas Fx " << fx << std::endl;
                        std::cout << "FFmaBlockBlas Fy " << fy << std::endl;
                        std::cout << "FFmaBlockBlas Fz " << fz << std::endl;
                } // end FFmaBlas kernel
438 439 440

#endif

441
#ifdef  SCALFMM_USE_FFT
442 443 444 445
        //
        ////////////////////////////////////////////////////////////////////
        //
        {	// begin Lagrange/Uniform Grid kernel
446

447
                // TODO
448

449 450 451
                // accuracy
                const unsigned int ORDER = 7;
                std::cout << "\nLagrange FMM ... ORDER " << ORDER <<std::endl;
452

453
                // typedefs
454

455 456 457 458 459 460
        typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
        typedef FSimpleLeaf<FReal,  ContainerClass >  LeafClass;
        typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
        typedef FUnifCell<FReal, ORDER> CellClass;
        typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
        typedef FUnifKernel<FReal, CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
461
                typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
462 463


464 465
                // init oct-tree
                OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
466 467 468 469

    // Create Matrix Kernel
    const MatrixKernelClass MatrixKernel;

470 471
                { // -----------------------------------------------------
                                time.tic();
472

473
            for(FSize idxPart = 0 ; idxPart <nbParticles ; ++idxPart){
474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491
                                // put in tree
                                tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
                        }

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

                { // -----------------------------------------------------
                        time.tic();
                        KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
                        FmmClass algorithm(&tree, &kernels);
                        algorithm.execute();
                        time.tac();
                        std::cout <<  "(Lagrange @Algorithm = " << time.elapsed() << " s)." << std::endl;
                } // -----------------------------------------------------

                FReal energy = 0.0;
492 493
        FMath::FAccurater<FReal> potentialDiff;
        FMath::FAccurater<FReal> fx, fy, fz;
494 495 496 497 498 499 500 501
                { // Check that each particle has been summed with all other

                        tree.forEachLeaf([&](LeafClass* leaf){
                                const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
                                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();
502 503
                const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
                const FVector<FSize>& indexes       = leaf->getTargets()->getIndexes();
504

505 506
                for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const FSize indexPartOrig = indexes[idxPart];
507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523
                                        potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
                                        fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
                                        fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
                                        fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
                                        energy += potentials[idxPart]*physicalValues[idxPart] ;
                                }
                        });
                }

                // Print for information
                std::cout << "Lagrange Energy "  << FMath::Abs(energy-energyD) /energyD << std::endl;
                std::cout << "Lagrange Potential " << potentialDiff << std::endl;
                std::cout << "Lagrange Fx " << fx << std::endl;
                std::cout << "Lagrange Fy " << fy << std::endl;
                std::cout << "Lagrange Fz " << fz << std::endl;

        } // end Lagrange/Uniform Grid kernel
524
#endif
525 526 527 528 529
        //
        //         Spherical approximation
        //
        {
                //const static int P = 10;
530 531 532 533 534
        typedef FSphericalCell<FReal>               CellClass;
        typedef FP2PParticleContainerIndexed<FReal>          ContainerClass;
        typedef FSimpleLeaf<FReal,  ContainerClass >                     LeafClass;
        typedef FOctree<FReal,  CellClass, ContainerClass , LeafClass >  OctreeClass;
        typedef FSphericalKernel<FReal,  CellClass, ContainerClass >     KernelClass;
535
                typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
536

537
#ifndef  SCALFMM_USE_BLAS
538
                CellClass::Init(DevP, true);
539
#endif
540 541
                OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
                std::cout << "\nFFmaSpherical FMM ... P: " << DevP<< std::endl;
542

543 544
                { // -----------------------------------------------------
                        time.tic();
545

546
            for(FSize idxPart = 0 ; idxPart < nbParticles; ++idxPart){
547 548 549
                                // put in tree
                                tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
                        }
550

551 552 553
                        time.tac();
                        std::cout << "(FFmaSpherical @Inserting Particles = "<< time.elapsed() << " s)." << std::endl;
                } // -----------------------------------------------------
554

555 556 557 558 559
                // -----------------------------------------------------
                time.tic();
            CellClass::Init(DevP);
                KernelClass *kernels = new KernelClass(DevP, TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
                FmmClass algorithm(&tree, kernels);
560 561


562 563 564 565
                algorithm.execute();
                time.tac();
                std::cout << "(FFmaSpherical @Algorithm = " << time.elapsed() << " s)." << std::endl;
                // -----------------------------------------------------
566

567
                FReal energy = 0.0;
568 569
        FMath::FAccurater<FReal> potentialDiff;
        FMath::FAccurater<FReal> fx, fy, fz;
570 571 572 573 574 575 576 577
                { // Check that each particle has been summed with all other

                        tree.forEachLeaf([&](LeafClass* leaf){
                                const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
                                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();
578 579
                const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
                const FVector<FSize>& indexes       = leaf->getTargets()->getIndexes();
580

581 582
                for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const FSize indexPartOrig = indexes[idxPart];
583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604
                                        potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
                                        fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
                                        fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
                                        fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
                                        energy += potentials[idxPart]*physicalValues[idxPart] ;
                                }
                        });
                }

                // Print for information
                std::cout << "FFmaSpherical Energy "  << FMath::Abs(energy-energyD) /energyD << std::endl;
                std::cout << "FFmaSpherical Potential " << potentialDiff << std::endl;
                std::cout << "FFmaSpherical Fx " << fx << std::endl;
                std::cout << "FFmaSpherical Fy " << fy << std::endl;
                std::cout << "FFmaSpherical Fz " << fz << std::endl;

        }
        //
        //         Spherical Rotation
        //
        {
                const static int P = 11;
605 606 607 608 609
        typedef FRotationCell<FReal, P>               CellClass;
        typedef FP2PParticleContainerIndexed<FReal>          ContainerClass;
        typedef FSimpleLeaf<FReal,  ContainerClass >                     LeafClass;
        typedef FOctree<FReal,  CellClass, ContainerClass , LeafClass >  OctreeClass;
        typedef FRotationKernel<FReal,  CellClass, ContainerClass , P>   KernelClass;
610
                typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
611

612 613
                OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
                std::cout << "\nFFmaRotation FMM ... P: " << P<< std::endl;
614

615 616
                { // -----------------------------------------------------
                        time.tic();
617

618
            for(FSize idxPart = 0 ; idxPart < nbParticles; ++idxPart){
619 620 621
                                // put in tree
                                tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
                        }
622

623 624 625
                        time.tac();
                        std::cout << "(FFmaRotation @Inserting Particles = "<< time.elapsed() << " s)." << std::endl;
                } // -----------------------------------------------------
626

627 628 629 630
                // -----------------------------------------------------
                time.tic();
                KernelClass *kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
           FmmClass algorithm(&tree, kernels);
631 632


633 634 635 636
                algorithm.execute();
                time.tac();
                std::cout << "(FFmaRotation @Algorithm = " << time.elapsed() << " s)." << std::endl;
                // -----------------------------------------------------
637

638
                FReal energy = 0.0;
639 640
        FMath::FAccurater<FReal> potentialDiff;
        FMath::FAccurater<FReal> fx, fy, fz;
641 642 643 644 645 646 647 648
                { // Check that each particle has been summed with all other

                        tree.forEachLeaf([&](LeafClass* leaf){
                                const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
                                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();
649 650
                const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
                const FVector<FSize>& indexes       = leaf->getTargets()->getIndexes();
651

652 653
                for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const FSize indexPartOrig = indexes[idxPart];
654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678
                                        potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
                                        fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
                                        fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
                                        fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
                                        energy += potentials[idxPart]*physicalValues[idxPart] ;
                                }
                        });
                }

                // Print for information
                std::cout << "FFmaRotation Energy "  << FMath::Abs(energy-energyD) /energyD << std::endl;
                std::cout << "FFmaRotation Potential " << potentialDiff << std::endl;
                std::cout << "FFmaRotation Fx " << fx << std::endl;
                std::cout << "FFmaRotation Fy " << fy << std::endl;
                std::cout << "FFmaRotation Fz " << fz << std::endl;

        }

        ////////////////////////////////////////////////////////////////////
        {	// begin Taylor kernel

                // accuracy
                const unsigned int ORDER = 10;

                // typedefs
679
        typedef FTaylorCell<FReal, ORDER,1>                                 CellClass;
680
                std::cout << "\nFFmaTaylor FMM ... ORDER: " << ORDER << std::endl;
681

682 683 684 685
        typedef FP2PParticleContainerIndexed<FReal>                          ContainerClass;
        typedef FSimpleLeaf<FReal,  ContainerClass >                         LeafClass;
        typedef FOctree<FReal,  CellClass, ContainerClass , LeafClass >      OctreeClass;
        typedef FTaylorKernel<FReal, CellClass,ContainerClass,ORDER,1>       KernelClass;
686
                typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
687

688 689
                // init cell class and oct-tree
                OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
690

691 692
                { // -----------------------------------------------------
                        time.tic();
693

694
            for(FSize idxPart = 0 ; idxPart <nbParticles ; ++idxPart){
695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712
                                // put in tree
                                tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
                        }

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

                // -----------------------------------------------------
                time.tic();
                KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
                FmmClass algorithm(&tree, &kernels);
                algorithm.execute();
                time.tac();
                std::cout << "(FFmaTaylor @Algorithm = " << time.elapsed() << " s)." << std::endl;
                // -----------------------------------------------------

                FReal energy = 0.0;
713 714
        FMath::FAccurater<FReal> potentialDiff;
        FMath::FAccurater<FReal> fx, fy, fz;
715 716 717 718 719 720 721 722
                { // Check that each particle has been summed with all other

                        tree.forEachLeaf([&](LeafClass* leaf){
                                const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
                                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();
723 724
                const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
                const FVector<FSize>& indexes       = leaf->getTargets()->getIndexes();
725

726 727
                for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const FSize indexPartOrig = indexes[idxPart];
728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746
                                        potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
                                        fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
                                        fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
                                        fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
                                        energy += potentials[idxPart]*physicalValues[idxPart] ;
                                }
                        });
                }

                // Print for information
                std::cout << "FFmaTaylor Energy "  << FMath::Abs(energy-energyD) /energyD << std::endl;
                std::cout << "FFmaTaylor Potential " << potentialDiff << std::endl;
                std::cout << "FFmaTaylor Fx " << fx << std::endl;
                std::cout << "FFmaTaylor Fy " << fy << std::endl;
                std::cout << "FFmaTaylor Fz " << fz << std::endl;
        } // end FFTaylor kernel
        delete[] particles;

        return 0;
747 748

}