testRotationPeriodicBench.cpp 8.12 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 25 26 27 28 29 30 31 32 33 34 35 36 37
// ===================================================================================
#include "../../Src/Containers/FOctree.hpp"
#include "../../Src/Containers/FVector.hpp"

#include "../../Src/Kernels/Rotation/FRotationCell.hpp"
#include "../../Src/Kernels/Rotation/FRotationKernel.hpp"
#include "../../Src/Kernels/Rotation/FRotationOriginalKernel.hpp"
#include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp"

#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Core/FFmmAlgorithmPeriodic.hpp"
#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Files/FRandomLoader.hpp"

#include "../../Src/Components/FTestCell.hpp"
#include "../../Src/Components/FTestKernels.hpp"

#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Utils/FMath.hpp"
38 39 40 41
#include "../../Src/Utils/FParameterNames.hpp"



42 43 44 45 46
int main(int argc, char** argv){
    const FParameterNames LocalOptionMinPer {
        {"-min"},
        "The starting periodicity"
    };
47

48 49 50 51
    const FParameterNames LocalOptionMaxPer {
        {"-max"},
        "The ending periodicity"
    };
52

53 54 55 56
    const FParameterNames LocalOptionNbPrint {
        {"-nbprint"},
        "Number of result to print"
    };
57 58 59 60 61 62 63
    FHelpDescribeAndExit(argc, argv,
                         "Run a Spherical Harmonic (Rotation) FMM kernel in periodic.",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight,LocalOptionMinPer,
                         LocalOptionMaxPer, LocalOptionNbPrint
                         );

64 65 66
    /////////////////////////////////////////////////
    // Types
    /////////////////////////////////////////////////
67
    typedef double FReal;
68
    static const int P = 6;
69 70 71 72 73
    typedef FRotationCell<FReal,P>            CellClass;
    typedef FP2PParticleContainerIndexed<FReal>  ContainerClass;
    typedef FRotationKernel<FReal, CellClass, ContainerClass, P >   KernelClass;
    typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
    typedef FOctree<FReal,CellClass,ContainerClass,LeafClass> OctreeClass;
74
    typedef FFmmAlgorithmPeriodic<FReal,OctreeClass , CellClass, ContainerClass, KernelClass, LeafClass > FmmClassPer;
75 76 77 78 79 80 81 82

    /////////////////////////////////////////////////
    // Parameters
    /////////////////////////////////////////////////

    // Parameters
    const int NbLevels       = 4;
    const int SizeSubLevels  = 2;
83 84
    const int MinLevelAbove  = FParameters::getValue(argc, argv, LocalOptionMinPer.options,-1);
    const int MaxLevelAbove  = FParameters::getValue(argc, argv, LocalOptionMaxPer.options,3);
85
    const int IncLevelAbove  = 1;
86 87
    const FSize NbParticles    = FParameters::getValue(argc, argv, FParameterDefinitions::NbParticles.options,FSize(6));
    const int NbParticlesPrint    = FParameters::getValue(argc, argv, LocalOptionNbPrint.options, 6);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
88
    FAssertLF(NbParticlesPrint <= NbParticles , "The number of printer particles cannot be higer than the number of particles.");
89

BRAMAS Berenger's avatar
BRAMAS Berenger committed
90
    std::cout << "The application will use " << NbParticles << " but studies only " << NbParticlesPrint << " of them." << std::endl;
91 92 93 94 95

    /////////////////////////////////////////////////
    // Insert particlePositions in tree and copy into array
    /////////////////////////////////////////////////

96 97
    FRandomLoader<FReal> loader(NbParticles);
    FPoint<FReal>* const particlePositions = new FPoint<FReal>[NbParticles];
98

99
    for(FSize idxPart = 0 ; idxPart < NbParticles ; ++idxPart){
100 101 102
        loader.fillParticle(&particlePositions[idxPart]);
    }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
103
    FReal* allPotential = new FReal[(MaxLevelAbove-MinLevelAbove)*NbParticlesPrint];
104 105 106 107 108 109 110 111

    /////////////////////////////////////////////////
    // Test for different periodicity for FMM!!
    /////////////////////////////////////////////////
    for(int idxLevelAbove = MinLevelAbove; idxLevelAbove < MaxLevelAbove ; idxLevelAbove += IncLevelAbove){
        OctreeClass treePer(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());

        // insert in tree
112
        for(FSize idxPart = 0 ; idxPart < NbParticles ; ++idxPart){
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
            // put in tree
            treePer.insert(particlePositions[idxPart], idxPart, (idxPart&1?0.00010:-0.00010));
        }

        // Run FMM
        FmmClassPer algoPer(&treePer,idxLevelAbove);
        KernelClass kernelsPer( algoPer.extendedTreeHeight(), algoPer.extendedBoxWidth(), algoPer.extendedBoxCenter());
        algoPer.setKernel(&kernelsPer);
        algoPer.execute();

        FTreeCoordinate min, max;
        algoPer.repetitionsIntervals(&min,&max);
        std::cout << "Nb Levels above root = " << idxLevelAbove << "\n";
        std::cout << "\t repetitions       = " << algoPer.theoricalRepetition() << std::endl;
        std::cout << "\t properties was height " << treePer.getHeight() << ", width " << loader.getBoxWidth()
                  << ", center " << loader.getCenterOfBox() << "\n";
        std::cout << "\t new properties height " << algoPer.extendedTreeHeight() << ", width " << algoPer.extendedBoxWidth()
                  << ", center " << algoPer.extendedBoxCenter() << "\n";
        std::cout << "min = " << min << "\n";
        std::cout << "max = " << max << "\n";
        std::cout << std::endl;

        treePer.forEachLeaf([&](LeafClass* leaf){
            const FReal*const potentials = leaf->getTargets()->getPotentials();
137 138
            const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
            for(FSize idxPart = 0 ; idxPart < leaf->getTargets()->getNbParticles() ; ++idxPart){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
139 140 141
                if( indexes[idxPart] < NbParticlesPrint){
                    allPotential[(idxLevelAbove - MinLevelAbove)*NbParticlesPrint + indexes[idxPart]] = potentials[idxPart];
                }
142 143 144 145 146 147 148 149 150 151 152 153 154 155
            }
        });
    }

    /////////////////////////////////////////////////
    // Print results
    /////////////////////////////////////////////////


    std::cout << "Part" << "\t";
    for(int idxLevelAbove = MinLevelAbove; idxLevelAbove < MaxLevelAbove ; idxLevelAbove += IncLevelAbove){
        std::cout << idxLevelAbove << "\t";
    }
    std::cout << "\n";
156
    for(FSize idxPart = 0 ; idxPart < NbParticlesPrint ; ++idxPart){
157 158
        std::cout << idxPart << "\t";
        for(int idxLevelAbove = MinLevelAbove; idxLevelAbove < MaxLevelAbove ; idxLevelAbove += IncLevelAbove){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
159
            std::cout << allPotential[(idxLevelAbove - MinLevelAbove)*NbParticlesPrint + idxPart] << "\t";
160 161 162 163 164 165 166 167 168 169
        }
        std::cout << "\n";
    }
    std::cout << "\n";

    std::cout << "Part" << "\t";
    for(int idxLevelAbove = MinLevelAbove+1; idxLevelAbove < MaxLevelAbove ; idxLevelAbove += IncLevelAbove){
        std::cout << idxLevelAbove << "/" << idxLevelAbove-1 << "\t";
    }
    std::cout << "\n";
170
    for(FSize idxPart = 0 ; idxPart < NbParticlesPrint ; ++idxPart){
171 172
        std::cout << idxPart << "\t";
        for(int idxLevelAbove = MinLevelAbove +1; idxLevelAbove < MaxLevelAbove ; idxLevelAbove += IncLevelAbove){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
173 174
            std::cout << FMath::Abs((allPotential[(idxLevelAbove - MinLevelAbove)*NbParticlesPrint + idxPart]
                         / allPotential[(idxLevelAbove - MinLevelAbove - 1)*NbParticlesPrint + idxPart])-1.0) << "\t";
175 176 177 178 179 180 181 182 183 184 185 186
        }
        std::cout << "\n";
    }

    delete[] allPotential;
    delete[] particlePositions;

    return 0;
}