testTaylor.cpp 7.83 KB
Newer Older
1
// See LICENCE file at project root
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

// ==== 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"
23
#include "../../Src/Utils/FParameterNames.hpp"
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
24 25 26
//
// taylor kernel
#include "../../Src/Kernels/Taylor/FTaylorCell.hpp"
27
#include "../../Src/Kernels/Taylor/FTaylorKernel.hpp"
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
28 29 30
//
#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp"
31
#include "../../Src/Kernels/Interpolation/FInterpMatrixKernel.hpp"
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
32 33 34 35 36 37 38 39 40 41

/**
 * 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[])
{
42 43 44 45 46
    FHelpDescribeAndExit(argc, argv,
                         "Run a Taylor FMM kernel and compare the accuracy with a direct computation.",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::NbThreads);

47
    typedef double FReal;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
48
    // get info from commandline
49 50 51 52
    const char* const filename       = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma");
    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());
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
53 54 55 56 57 58 59 60 61 62
#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{
63
        FPoint<FReal> position;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
64 65 66 67 68
        FReal forces[3];
        FReal physicalValue;
        FReal potential;
    };
    // open particle file
69
    FFmaScanfLoader<FReal> loader(filename);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
70 71 72
    if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");

    TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
73
    for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
74
        FPoint<FReal> position;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
75 76 77 78 79 80 81 82 83 84 85
        FReal physicalValue = 0.0;
        loader.fillParticle(&position,&physicalValue);
        // get copy
        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;
    }
    {
86 87
        for(FSize idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
            for(FSize idxOther =  idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
88
                FP2PR::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
89 90 91 92 93 94
                                      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],
95
                                      &particles[idxOther].forces[2],&particles[idxOther].potential);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
96 97 98 99 100 101 102 103 104 105 106 107
            }
        }
    }
    //
    ////////////////////////////////////////////////////////////////////
    //
    {	// begin Taylor kernel

        // accuracy
        const unsigned int ORDER = 7;

        // typedefs
108
    typedef FTaylorCell<FReal,ORDER,1>                                 CellClass;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
109

110 111 112
  typedef FP2PParticleContainerIndexed<FReal>                          ContainerClass;
        typedef FSimpleLeaf<FReal, ContainerClass >                         LeafClass;
        typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >      OctreeClass;
113
	typedef FTaylorKernel<FReal,CellClass,ContainerClass,ORDER,1>       KernelClass;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
114 115 116 117 118 119 120 121 122 123 124 125 126
        typedef FFmmAlgorithm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
	//  typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;

        // init cell class and 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();

127
            for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
                // 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;
        // -----------------------------------------------------

147 148
        FMath::FAccurater<FReal> potentialDiff;
        FMath::FAccurater<FReal> fx, fy, fz;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
149 150 151 152 153 154 155
        { // 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();
156 157 158 159
                const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
                for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const FSize indexPartOrig = indexes[idxPart];
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
                    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 FFTaylor kernel
    delete[] particles;

    return 0;
}