testChebAlgorithm.cpp 10.2 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 23 24
// ==== CMAKE =====
// @FUSE_BLAS
// ================

25 26 27 28
#include <iostream>

#include <cstdio>
#include <cstdlib>
29

30
#include "Files/FFmaGenericLoader.hpp"
31

32 33
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
34
#include "Kernels/Interpolation/FInterpMatrixKernel_Covariance.hpp"
35 36
#include "Kernels/Chebyshev/FChebKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
37

38 39
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
40

41 42
#include "Utils/FParameters.hpp"
#include "Utils/FMemUtils.hpp"
43

44 45
#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"
46

47
#include "Core/FFmmAlgorithm.hpp"
48
#ifdef _OPENMP
49
#include "Core/FFmmAlgorithmThread.hpp"
50
#endif
51

52 53
#include "../../Src/Utils/FParameterNames.hpp"

54 55 56
/**
 * This program runs the FMM Algorithm with the Chebyshev kernel and compares the results with a direct computation.
 */
57 58 59 60

// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
61 62 63 64 65
    FHelpDescribeAndExit(argc, argv,
                         "Test the chebyshev FMM and compare it to the direct computation for debugging purpose.",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::NbThreads);

66
    typedef double FReal;
67 68 69 70
    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, 4);
    const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
    const unsigned int NbThreads     = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
71
    FTic time;
72 73

#ifdef _OPENMP
74 75
    omp_set_num_threads(NbThreads);
    std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
76
#else
77
    std::cout << "\n>> Sequential version.\n" << std::
78 79
#endif

80
    // interaction kernel evaluator
81
    typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
82
    const MatrixKernelClass MatrixKernel;
83 84
    //typedef FInterpMatrixKernelGauss<FReal> MatrixKernelClass;
    //const MatrixKernelClass MatrixKernel(.5);
85 86 87

    // init particles position and physical value
    struct TestParticle{
88
        FPoint<FReal> position;
89 90 91 92 93 94
        FReal forces[3];
        FReal physicalValue;
        FReal potential;
    };

    // open particle file
95
    FFmaGenericLoader<FReal> loader(filename);
96 97 98
    if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");

    TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
99
    for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
100
        FPoint<FReal> position;
101 102 103 104 105 106 107 108 109
        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;
110
    }
111

112 113 114 115 116 117
    ////////////////////////////////////////////////////////////////////

    { // begin direct computation

        time.tic();
        {
118 119
            for(FSize idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
                for(FSize idxOther =  idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
120 121 122
                    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],
123 124 125 126 127
                                          &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,&MatrixKernel);
128 129 130 131 132 133 134
                }
            }
        }
        time.tac();
        printf("Elapsed Time for direct computation: %f\n",time.elapsed());

    } // end direct computation
135

136 137 138
    ////////////////////////////////////////////////////////////////////

    {	// begin Chebyshev kernel
139 140 141 142 143

        // accuracy
        const unsigned int ORDER = 7 ;

        // typedefs
144 145 146 147
        typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
        typedef FSimpleLeaf<FReal, ContainerClass >  LeafClass;
        typedef FChebCell<FReal,ORDER> CellClass;
        typedef FOctree<FReal,CellClass,ContainerClass,LeafClass> OctreeClass;
148 149
        typedef FChebKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
        //typedef FChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
#ifdef _OPENMP
        typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#else
        typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#endif

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

166
            for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
167 168 169 170 171 172 173 174 175 176 177 178
                // 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 (ORDER="<< ORDER << "... " << std::endl;
            time.tic();
179
            KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel,FReal(1.e-15));
180 181 182 183 184 185 186 187 188
            FmmClass algorithm(&tree, &kernels);
            algorithm.execute();
            time.tac();
            std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
        } // -----------------------------------------------------


        { // -----------------------------------------------------
            std::cout << "\nError computation ... " << std::endl;
189 190
            FMath::FAccurater<FReal> potentialDiff;
            FMath::FAccurater<FReal> fx, fy, fz;
191 192 193 194 195 196 197
            { // 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();
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
                        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 Chebyshev kernel
219

220 221 222
    std::cout << "Memory used at the end " << FMemStats::controler.getCurrentAllocated() << " Bytes (" << FMemStats::controler.getCurrentAllocatedMB() << "MB)\n";
    std::cout << "Max memory used " << FMemStats::controler.getMaxAllocated() << " Bytes (" << FMemStats::controler.getMaxAllocatedMB() << "MB)\n";
    std::cout << "Total memory used " << FMemStats::controler.getTotalAllocated() << " Bytes (" << FMemStats::controler.getTotalAllocatedMB() << "MB)\n";
223

224
    return 0;
225
}