testTsmAlgorithm.cpp 8.05 KB
Newer Older
1
// ===================================================================================
2 3 4 5
// Copyright ScalFmm 2016 INRIA, Olivier Coulaud, Bérenger Bramas,
// Matthias Messner olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the
// FMM.
6
//
7
// This software is governed by the CeCILL-C and LGPL licenses and
8
// abiding by the rules of distribution of free software.
9 10 11
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
12
//
13 14 15
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 17 18
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
19
// ===================================================================================
20 21 22

#include <iostream>

23 24
#include <cstdio>
#include <cstdlib>
25

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

29
#include "../../Src/Components/FTypedLeaf.hpp"
30

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

34
#include "../../Src/Core/FFmmAlgorithmTsm.hpp"
35

36 37
#include "../../Src/Kernels/Spherical/FSphericalKernel.hpp"
#include "../../Src/Kernels/Spherical/FSphericalCell.hpp"
38

39 40 41
#include "../../Src/Kernels/Rotation/FRotationKernel.hpp"
#include "../../Src/Kernels/Rotation/FRotationCell.hpp"

42
#include "../../Src/Files/FFmaTsmLoader.hpp"
43

44 45
#include "../../Src/Kernels/P2P/FP2PParticleContainer.hpp"

46 47
#include "../../Src/Utils/FParameterNames.hpp"

48 49
/** This program show an example of use of
  * the fmm basic algo
50
  * it also check that eachh particles is little or longer
51 52 53 54
  * related that each other
  */


55
// Simply create particles and try the kernels
56
template <class FReal, class CellClass, class ContainerClass, class LeafClass, class OctreeClass,
57
          class KernelClass, class FmmClass, typename... Args>
58
int testFunction(int argc, char ** argv, Args ... kernelPreArgs){
59 60
    FTic counter;
    // Retrieve parameters
61 62
    const int NbLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 5);
    const int SizeSubLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 3);
63
    // Get working file
64
    const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.tsm.fma");
65
    std::cout << "Opening : " << filename << "\n";
66
    // Create particles loader
67
    FFmaTsmLoader<FReal> loader(filename);
68
    if(!loader.isOpen()){
69 70 71 72 73
        std::cout << "Loader Error, " << filename << " is missing\n";
        return 1;
    }

    // -----------------------------------------------------
74
    // Build the tree
berenger-bramas's avatar
berenger-bramas committed
75
    OctreeClass tree(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
76 77 78

    // -----------------------------------------------------

79
    std::cout << "Creating & Inserting " << loader.getNumberOfParticles() << " particles ..." << std::endl;
80
    std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
81 82
    counter.tic();

83
    for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
84
        FPoint<FReal> particlePosition;
85
        FReal physicalValue = 0.0;
86 87 88
        FParticleType particleType;
        loader.fillParticle(&particlePosition,&physicalValue,&particleType);
        tree.insert(particlePosition, particleType, physicalValue );
89
    }
90 91

    counter.tac();
92
    std::cout << "Done  " << "(@Creating and Inserting Particles = " << counter.elapsed() << "s)." << std::endl;
93 94 95

    // -----------------------------------------------------

96
    std::cout << "Create kernel ..." << std::endl;
97
    counter.tic();
98 99
    //    KernelClass kernels( kernelPreArgs... , NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
    KernelClass kernels( kernelPreArgs... , loader.getBoxWidth(), loader.getCenterOfBox());
berenger-bramas's avatar
berenger-bramas committed
100

101 102 103 104 105 106 107
    counter.tac();
    std::cout << "Done  " << " in " << counter.elapsed() << "s)." << std::endl;

    // -----------------------------------------------------

    std::cout << "Working on particles ..." << std::endl;

berenger-bramas's avatar
berenger-bramas committed
108
    FmmClass algo(&tree,&kernels);
109

110 111
    counter.tic();
    algo.execute();
112
    counter.tac();
113

114
    std::cout << "Done  " << "(@Algorithm = " << counter.elapsed() << "s)." << std::endl;
115 116 117

    { // get sum forces&potential
        FReal potential = 0;
118 119 120 121 122 123 124
        FReal fx = 0.0, fy = 0.0, fz = 0.0;

        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();
125
            const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
126

127
            for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
128 129 130 131
                potential += potentials[idxPart];
                fx += forcesX[idxPart];
                fy += forcesY[idxPart];
                fz += forcesZ[idxPart];
132
            }
133
        });
134

135
        std::cout << "Foces Sum  x = " << fx << " y = " << fy << " z = " << fz << std::endl;
136 137 138 139 140 141 142 143 144
        std::cout << "Potential = " << potential << std::endl;
    }


    // -----------------------------------------------------

    return 0;
}

145 146
// This is the real main!
int main(int argc, char ** argv){
147 148 149 150 151 152
    FHelpDescribeAndExit(argc, argv,
                         "Test the TSM (target source model) using the Rotation or the Spherical Harmonic Old implementations.",
                         FParameterDefinitions::OctreeHeight,FParameterDefinitions::SHDevelopment,
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
                         FParameterDefinitions::SphericalKernel, FParameterDefinitions::RotationKernel);

153
    std::cout << "[PARAM] Use Parameters -spherical -rotation -chebyshev\n";
154
    const int NbLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 5);
155

156
    if( FParameters::existParameter(argc,argv,FParameterDefinitions::SphericalKernel.options) ){
157 158
        std::cout << "[INFO] -spherical is used\n";
        // Create template
159
        typedef double FReal;
160
        typedef FTypedSphericalCell< FReal>            CellClass;
161
        typedef FP2PParticleContainer<FReal>         ContainerClass;
162

163
        typedef FTypedLeaf< FReal, ContainerClass >                      LeafClass;
164
        typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
165
        typedef FSphericalKernel< FReal, CellClass, ContainerClass >          KernelClass;
166 167

        typedef FFmmAlgorithmTsm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
168

169
        const int DevP = FParameters::getValue(argc,argv,FParameterDefinitions::SHDevelopment.options, 8);
170
        CellClass::Init(DevP);
171

172
        // Call Main function
173
        testFunction< FReal, CellClass, ContainerClass, LeafClass, OctreeClass, KernelClass, FmmClass>(argc, argv, DevP,NbLevels);
174 175
    }

176
    if( FParameters::existParameter(argc,argv,FParameterDefinitions::RotationKernel.options) ){
177 178
        std::cout << "[INFO] -rotation is used\n";
        // Create template
179
        typedef double FReal;
180
        static const int P = 9;
181
        typedef FTypedRotationCell<FReal,P>            CellClass;
182
        typedef FP2PParticleContainer<FReal>         ContainerClass;
183

184
        typedef FTypedLeaf< FReal, ContainerClass >                      LeafClass;
185
        typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
186
        typedef FRotationKernel< FReal, CellClass, ContainerClass, P >          KernelClass;
187 188 189 190

        typedef FFmmAlgorithmTsm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;

        // Call Main function
191
        testFunction< FReal, CellClass, ContainerClass, LeafClass, OctreeClass, KernelClass, FmmClass>(argc, argv,NbLevels);
192 193 194 195
    }

    return 0;
}