testCompareKernels.cpp 12.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
// ===================================================================================
// 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
// ===================================================================================

17 18 19 20
// ==== CMAKE =====
// @FUSE_BLAS
// ================

21 22 23 24 25
#include <iostream>

#include <cstdio>
#include <cstdlib>

26 27
#include "../../Src/Utils/FTic.hpp"
#include "../../Src/Utils/FParameters.hpp"
28

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

31 32
#include "../../Src/Containers/FOctree.hpp"
#include "../../Src/Containers/FVector.hpp"
33

34 35
#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Core/FFmmAlgorithmThread.hpp"
36 37

// chebyshev kernel
38

39
#include "../../Src/Kernels/Chebyshev/FChebCell.hpp"
40
#include "../../Src/Kernels/Interpolation/FInterpMatrixKernel.hpp"
41 42
#include "../../Src/Kernels/Chebyshev/FChebKernel.hpp"
#include "../../Src/Kernels/Chebyshev/FChebSymKernel.hpp"
43 44

// spherical kernel
45 46 47 48 49
#include "../../Src/Kernels/Spherical/FSphericalKernel.hpp"
#include "../../Src/Kernels/Spherical/FSphericalBlasKernel.hpp"
#include "../../Src/Kernels/Spherical/FSphericalBlockBlasKernel.hpp"
#include "../../Src/Kernels/Spherical/FSphericalRotationKernel.hpp"
#include "../../Src/Kernels/Spherical/FSphericalCell.hpp"
50

51
#include "../../Src/Components/FSimpleLeaf.hpp"
52
#include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp"
53 54 55

/**
 * This program compares two different kernels, eg., the Chebyshev kernel with
berenger-bramas's avatar
berenger-bramas committed
56
 * the SphericalBlas kernel.
57 58 59 60 61 62
 */


// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
63
    // get info from commandline
64
    const char* const filename       = FParameters::getStr(argc,argv,"-f", "../Data/test20k.fma");
65 66
    const unsigned int TreeHeight    = FParameters::getValue(argc, argv, "-h", 5);
    const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, "-sh", 2);
67 68 69 70 71 72
    const unsigned int NbThreads     = FParameters::getValue(argc, argv, "-t", omp_get_max_threads());

    omp_set_num_threads(NbThreads);

    std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;

73 74 75
    // init timer
    FTic time;

76 77 78 79 80 81 82 83 84 85 86 87 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
    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
        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;
    }
    {
        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);
            }
        }
    }
113

114 115 116 117
    ////////////////////////////////////////////////////////////////////
    {	// begin Chebyshef kernel

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

        // typedefs
122
        typedef FP2PParticleContainerIndexed<> ContainerClass;
123
        typedef FSimpleLeaf<ContainerClass> LeafClass;
124
        typedef FInterpMatrixKernelR MatrixKernelClass;
125
        typedef FChebCell<ORDER> CellClass;
126
        typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
127

128 129 130 131
        //typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
        typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
        //typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
        typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
132

133 134 135 136 137

        // init oct-tree
        OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());

        { // -----------------------------------------------------
138
            std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
139
                      << " particles ..." << std::endl;
140 141
            std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
            time.tic();
142 143

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

148
            time.tac();
149
            std::cout << "Done  " << "(@Creating and Inserting Particles = "
150
                      << time.elapsed() << "s)." << std::endl;
151 152 153 154 155
        } // -----------------------------------------------------

        { // -----------------------------------------------------
            std::cout << "\nChebyshev FMM ... " << std::endl;
            time.tic();
156
            KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(), epsilon);
157 158 159 160 161 162
            FmmClass algorithm(&tree, &kernels);
            algorithm.execute();
            time.tac();
            std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
        } // -----------------------------------------------------

163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
        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]);
181
                }
182 183
            });
        }
184

185 186
        // Print for information
        std::cout << "Potential " << potentialDiff << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
187 188 189
        std::cout << "Fx " << fx << std::endl;
        std::cout << "Fy " << fy << std::endl;
        std::cout << "Fz " << fz << std::endl;
190 191 192 193 194 195 196 197

    } // end Chebyshev kernel


    ////////////////////////////////////////////////////////////////////
    {	// begin FFmaBlas kernel

        // accuracy
BRAMAS Berenger's avatar
BRAMAS Berenger committed
198
        const int DevP = FParameters::getValue(argc, argv, "-p", 6);
199 200 201

        // typedefs
        typedef FSphericalCell                 CellClass;
202
        typedef FP2PParticleContainerIndexed<>         ContainerClass;
203 204 205 206 207 208 209 210 211
        typedef FSimpleLeaf< ContainerClass >                     LeafClass;
        typedef FOctree< CellClass, ContainerClass , LeafClass >  OctreeClass;
        //typedef FSphericalBlasKernel< CellClass, ContainerClass > KernelClass;
        //typedef FSphericalBlockBlasKernel< CellClass, ContainerClass > KernelClass;
        //typedef FSphericalKernel< CellClass, ContainerClass > KernelClass;
        typedef FSphericalBlockBlasKernel< CellClass, ContainerClass > KernelClass;
        //typedef FSphericalRotationKernel< CellClass, ContainerClass > KernelClass;
        //typedef FFmmAlgorithm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
        typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
212 213

        // init cell class and oct-tree
214 215
        CellClass::Init(DevP, true); // only for blas
        CellClass::Init(DevP, false);
216 217 218
        OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());

        { // -----------------------------------------------------
219
            std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
220
                      << " particles ..." << std::endl;
221
            std::cout << "\tHeight : " << TreeHeight << " \t sub-height : "
222
                      << SubTreeHeight << std::endl;
223
            time.tic();
224 225

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

230
            time.tac();
231
            std::cout << "Done  " << "(@Creating and Inserting Particles = "
232
                      << time.elapsed() << "s)." << std::endl;
233 234 235 236 237 238 239 240 241 242 243 244
        } // -----------------------------------------------------

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

245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
        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]);
263
                }
264 265
            });
        }
266

267 268
        // Print for information
        std::cout << "Potential " << potentialDiff << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
269 270 271
        std::cout << "Fx " << fx << std::endl;
        std::cout << "Fy " << fy << std::endl;
        std::cout << "Fz " << fz << std::endl;
272 273
    } // end FFmaBlas kernel

274
    delete[] particles;
275 276

    return 0;
277
}