testNewCompareKernels.cpp 23.1 KB
Newer Older
COULAUD Olivier's avatar
COULAUD Olivier committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
// ===================================================================================
// Ce LOGICIEL "ScalFmm" est couvert par le copyright Inria 20xx-2012.
// Inria détient tous les droits de propriété sur le LOGICIEL, et souhaite que
// la communauté scientifique l'utilise afin de le tester et de l'évaluer.
// Inria donne gracieusement le droit d'utiliser ce LOGICIEL. Toute utilisation
// dans un but lucratif ou à des fins commerciales est interdite sauf autorisation
// expresse et préalable d'Inria.
// Toute utilisation hors des limites précisées ci-dessus et réalisée sans l'accord
// expresse préalable d'Inria constituerait donc le délit de contrefaçon.
// Le LOGICIEL étant un produit en cours de développement, Inria ne saurait assurer
// aucune responsabilité et notamment en aucune manière et en aucun cas, être tenu
// de répondre d'éventuels dommages directs ou indirects subits par l'utilisateur.
// Tout utilisateur du LOGICIEL s'engage à communiquer à Inria ses remarques
// relatives à l'usage du LOGICIEL
// ===================================================================================

// ==== CMAKE =====
//
// ================

#include <iostream>

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

#include "../../Src/Files/FFmaScanfLoader.hpp"

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

#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Core/FFmmAlgorithmThread.hpp"

#ifdef ScalFMM_USE_BLAS
// chebyshev kernel

#include "../../Src/Kernels/Chebyshev/FChebCell.hpp"
42
#include "../../Src/Kernels/Interpolation/FInterpMatrixKernel.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
43 44 45 46 47 48 49 50 51 52 53 54 55
#include "../../Src/Kernels/Chebyshev/FChebKernel.hpp"
#include "../../Src/Kernels/Chebyshev/FChebSymKernel.hpp"
#endif
//
// spherical kernel
#include "../../Src/Kernels/Spherical/FSphericalCell.hpp"
#ifdef ScalFMM_USE_BLAS
#include "../../Src/Kernels/Spherical/FSphericalBlasKernel.hpp"
#include "../../Src/Kernels/Spherical/FSphericalBlockBlasKernel.hpp"
#endif
//
// taylor kernel
#include "../../Src/Kernels/Taylor/FTaylorCell.hpp"
56
#include "../../Src/Kernels/Taylor/FTaylorKernel.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
57 58 59 60
//
#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp"

61 62 63 64
//Rotation kernel
#include "../../Src/Kernels/Rotation/FRotationKernel.hpp"
#include "../../Src/Kernels/Rotation/FRotationCell.hpp"

65
#ifdef ScalFMM_USE_FFT
66 67 68 69 70 71
// Uniform grid kernel
#include "../../Src/Kernels/Uniform/FUnifCell.hpp"
//#include "../../Src/Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "../../Src/Kernels/Uniform/FUnifKernel.hpp"
#endif

72

COULAUD Olivier's avatar
COULAUD Olivier committed
73 74 75 76 77 78 79 80 81 82 83 84 85 86
/**
 * This program compares two different kernels, eg., the Chebyshev kernel with
 * the SphericalBlas kernel.
 */


// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
    // get info from commandline
    const char* const filename       = FParameters::getStr(argc,argv,"-f", "../Data/test20k.fma");
    const unsigned int TreeHeight    = FParameters::getValue(argc, argv, "-h", 5);
    const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, "-sh", 2);
    const unsigned int NbThreads     = FParameters::getValue(argc, argv, "-t", omp_get_max_threads());
COULAUD Olivier's avatar
COULAUD Olivier committed
87
    const int DevP                   = FParameters::getValue(argc, argv, "-p", 7);
COULAUD Olivier's avatar
COULAUD Olivier committed
88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
#ifdef _OPENMP
    omp_set_num_threads(NbThreads);
    std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
#else
    std::cout << "\n>> Sequential version.\n" << std::
#endif
    // init timer
    FTic time;

    struct TestParticle{
        FPoint position;
        FReal forces[3];
        FReal physicalValue;
        FReal potential;
    };
    // open particle file
    FFmaScanfLoader loader(filename);
    if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");

    TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
    for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
        FPoint position;
        FReal physicalValue = 0.0;
        loader.fillParticle(&position,&physicalValue);
        // get copy
COULAUD Olivier's avatar
COULAUD Olivier committed
113 114 115 116 117 118
        particles[idxPart].position       = position;
        particles[idxPart].physicalValue  = physicalValue;
        particles[idxPart].potential      = 0.0;
        particles[idxPart].forces[0]      = 0.0;
        particles[idxPart].forces[1]      = 0.0;
        particles[idxPart].forces[2]      = 0.0;
COULAUD Olivier's avatar
COULAUD Olivier committed
119
    }
120
    time.tic();
COULAUD Olivier's avatar
COULAUD Olivier committed
121 122 123 124 125 126 127 128 129 130 131 132 133 134
    {
        for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
            for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
                FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
                                      particles[idxTarget].position.getZ(),particles[idxTarget].physicalValue,
                                      &particles[idxTarget].forces[0],&particles[idxTarget].forces[1],
                                      &particles[idxTarget].forces[2],&particles[idxTarget].potential,
                                particles[idxOther].position.getX(), particles[idxOther].position.getY(),
                                particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
                                &particles[idxOther].forces[0],&particles[idxOther].forces[1],
                                &particles[idxOther].forces[2],&particles[idxOther].potential);
            }
        }
    }
135 136
    time.tac();
    printf("Elapsed Time for direct computation: %f\n",time.elapsed());
COULAUD Olivier's avatar
COULAUD Olivier committed
137 138 139 140
    //
    ////////////////////////////////////////////////////////////////////
    //
#ifdef  ScalFMM_USE_BLAS
141
    {	// begin Chebyshev kernel
COULAUD Olivier's avatar
COULAUD Olivier committed
142 143 144 145 146 147

        // accuracy
        const unsigned int ORDER = 7;
        const FReal epsilon = FReal(1e-7);

        // typedefs
148
        typedef FP2PParticleContainerIndexed<> ContainerClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
149
        typedef FSimpleLeaf<ContainerClass> LeafClass;
150
        typedef FInterpMatrixKernelR MatrixKernelClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
151 152 153 154
        typedef FChebCell<ORDER> CellClass;
        typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;

        typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
155
        //        typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
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
        typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;


        // 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();

            for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
                // put in tree
                tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
            }

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

        { // -----------------------------------------------------
            std::cout << "\nChebyshev FMM ... " << std::endl;
            time.tic();
181
            KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(), epsilon);
COULAUD Olivier's avatar
COULAUD Olivier committed
182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
            FmmClass algorithm(&tree, &kernels);
            algorithm.execute();
            time.tac();
            std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
        } // -----------------------------------------------------

        FMath::FAccurater potentialDiff;
        FMath::FAccurater fx, fy, fz;
        { // Check that each particle has been summed with all other

            tree.forEachLeaf([&](LeafClass* leaf){
                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();
                const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
                const FVector<int>& indexes = leaf->getTargets()->getIndexes();

                for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const int indexPartOrig = indexes[idxPart];
                    potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
                    fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
                    fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
                    fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
                }
            });
        }

        // Print for information
        std::cout << "Potential " << potentialDiff << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
212 213 214
        std::cout << "Fx " << fx << std::endl;
        std::cout << "Fy " << fy << std::endl;
        std::cout << "Fz " << fz << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
215 216 217 218

    } // end Chebyshev kernel
    //
    ////////////////////////////////////////////////////////////////////
COULAUD Olivier's avatar
COULAUD Olivier committed
219
    //
COULAUD Olivier's avatar
COULAUD Olivier committed
220 221 222 223
    {	// begin FFmaBlas kernel

        // typedefs
        typedef FSphericalCell                 CellClass;
224
        typedef FP2PParticleContainerIndexed<>         ContainerClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
225 226 227 228 229 230 231
        typedef FSimpleLeaf< ContainerClass >                     LeafClass;
        typedef FOctree< CellClass, ContainerClass , LeafClass >  OctreeClass;
        typedef FSphericalBlockBlasKernel< CellClass, ContainerClass > KernelClass;
        typedef FFmmAlgorithm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;

        // init cell class and oct-tree
        CellClass::Init(DevP, true); // only for blas
232
        //CellClass::Init(DevP, false);
COULAUD Olivier's avatar
COULAUD Olivier committed
233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
        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();

            for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
                // put in tree
                tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
            }

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

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

        FMath::FAccurater potentialDiff;
        FMath::FAccurater fx, fy, fz;
        { // Check that each particle has been summed with all other

            tree.forEachLeaf([&](LeafClass* leaf){
                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();
                const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
                const FVector<int>& indexes = leaf->getTargets()->getIndexes();

                for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const int indexPartOrig = indexes[idxPart];
                    potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
                    fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
                    fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
                    fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
                }
            });
        }

        // Print for information
        std::cout << "Potential " << potentialDiff << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
286 287 288
        std::cout << "Fx " << fx << std::endl;
        std::cout << "Fy " << fy << std::endl;
        std::cout << "Fz " << fz << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
289
    } // end FFmaBlas kernel
290 291 292
#endif

#ifdef  ScalFMM_USE_FFT
293 294 295 296 297 298 299 300 301 302 303
    //
    ////////////////////////////////////////////////////////////////////
    //
    {	// begin Lagrange/Uniform Grid kernel

      // TODO 

      // accuracy
      const unsigned int ORDER = 7;

      // typedefs
304
      typedef FP2PParticleContainerIndexed<> ContainerClass;
305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335
      typedef FSimpleLeaf<ContainerClass> LeafClass;
      typedef FInterpMatrixKernelR MatrixKernelClass;
      typedef FUnifCell<ORDER> CellClass;
      typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;

      typedef FUnifKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
      typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;


      // 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();

        for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
          // put in tree
          tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
        }

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

      { // -----------------------------------------------------
        std::cout << "\nLagrange FMM ... " << std::endl;
        time.tic();
336
        KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371
        FmmClass algorithm(&tree, &kernels);
        algorithm.execute();
        time.tac();
        std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
      } // -----------------------------------------------------

      FMath::FAccurater potentialDiff;
      FMath::FAccurater fx, fy, fz;
      { // Check that each particle has been summed with all other

        tree.forEachLeaf([&](LeafClass* leaf){
            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();
            const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
            const FVector<int>& indexes = leaf->getTargets()->getIndexes();

            for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
              const int indexPartOrig = indexes[idxPart];
              potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
              fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
              fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
              fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
            }
          });
      }

      // Print for information
      std::cout << "Potential " << potentialDiff << std::endl;
      std::cout << "Fx " << fx << std::endl;
      std::cout << "Fy " << fy << std::endl;
      std::cout << "Fz " << fz << std::endl;

    } // end Lagrange/Uniform Grid kernel
COULAUD Olivier's avatar
COULAUD Olivier committed
372 373
#endif

374 375 376
{
  const static int P = 7;
  typedef FRotationCell<P>               CellClass;
377
  typedef FP2PParticleContainerIndexed<>          ContainerClass;
378 379 380 381 382 383 384
  typedef FSimpleLeaf< ContainerClass >                     LeafClass;
  typedef FOctree< CellClass, ContainerClass , LeafClass >  OctreeClass;
  typedef FRotationKernel< CellClass, ContainerClass , P>   KernelClass;
  typedef FFmmAlgorithm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
  
  OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
  
COULAUD Olivier's avatar
COULAUD Olivier committed
385
        { // -----------------------------------------------------
386
	  printf("Rotation kernel\n");
COULAUD Olivier's avatar
COULAUD Olivier committed
387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
            std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
                      << " particles ..." << std::endl;
            std::cout << "\tHeight : " << TreeHeight << " \t sub-height : "
                      << SubTreeHeight << std::endl;
            time.tic();

            for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
                // put in tree
                tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
            }

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

        // -----------------------------------------------------
404
        std::cout << "\nFFmaRotation FMM ..." << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
405
        time.tic();
406
        KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
COULAUD Olivier's avatar
COULAUD Olivier committed
407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429
        FmmClass algorithm(&tree, &kernels);
        algorithm.execute();
        time.tac();
        std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
        // -----------------------------------------------------

        FMath::FAccurater potentialDiff;
        FMath::FAccurater fx, fy, fz;
        { // Check that each particle has been summed with all other

            tree.forEachLeaf([&](LeafClass* leaf){
                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();
                const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
                const FVector<int>& indexes = leaf->getTargets()->getIndexes();
                for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const int indexPartOrig = indexes[idxPart];
                    potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
                    fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
                    fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
                    fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
430
               }
COULAUD Olivier's avatar
COULAUD Olivier committed
431 432 433 434 435
            });
        }

        // Print for information
        std::cout << "Potential " << potentialDiff << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
436 437 438
        std::cout << "Fx " << fx << std::endl;
        std::cout << "Fy " << fy << std::endl;
        std::cout << "Fz " << fz << std::endl;
439
    }
COULAUD Olivier's avatar
COULAUD Olivier committed
440 441 442 443 444 445 446 447 448
    ////////////////////////////////////////////////////////////////////
    {	// begin Taylor kernel

        // accuracy
        const unsigned int ORDER = 7;

        // typedefs
	typedef FTaylorCell<ORDER,1>                                 CellClass;

449
        typedef FP2PParticleContainerIndexed<>                          ContainerClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
450 451 452 453 454 455 456
        typedef FSimpleLeaf< ContainerClass >                         LeafClass;
        typedef FOctree< CellClass, ContainerClass , LeafClass >      OctreeClass;
	typedef FTaylorKernel<CellClass,ContainerClass,ORDER,1>       KernelClass;
        typedef FFmmAlgorithm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;

        // init cell class and oct-tree
        OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
457
	
COULAUD Olivier's avatar
COULAUD Olivier committed
458
        { // -----------------------------------------------------
459 460
          printf("TAylor Kernel\n");
	  std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
COULAUD Olivier's avatar
COULAUD Olivier committed
461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502
                      << " particles ..." << std::endl;
            std::cout << "\tHeight : " << TreeHeight << " \t sub-height : "
                      << SubTreeHeight << std::endl;
            time.tic();

            for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
                // put in tree
                tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
            }

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

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

        FMath::FAccurater potentialDiff;
        FMath::FAccurater fx, fy, fz;
        { // Check that each particle has been summed with all other

            tree.forEachLeaf([&](LeafClass* leaf){
                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();
                const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
                const FVector<int>& indexes = leaf->getTargets()->getIndexes();
                for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const int indexPartOrig = indexes[idxPart];
                    potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
                    fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
                    fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
                    fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
COULAUD Olivier's avatar
COULAUD Olivier committed
503
               }
COULAUD Olivier's avatar
COULAUD Olivier committed
504 505 506 507 508
            });
        }

        // Print for information
        std::cout << "Potential " << potentialDiff << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
509 510 511
        std::cout << "Fx " << fx << std::endl;
        std::cout << "Fy " << fy << std::endl;
        std::cout << "Fz " << fz << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
512 513 514 515 516
    } // end FFTaylor kernel
    delete[] particles;

    return 0;
}