testSmallCase.cpp 9.25 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
// Keep in private GIT
// @SCALFMM_PRIVATE

24
// @FUSE_FFT
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
#include <iostream>
#include <cstdio>


#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Utils/FTic.hpp"

#include "../../Src/Containers/FOctree.hpp"
#include "../../Src/Containers/FVector.hpp"

#include "../../Src/Components/FSimpleIndexedLeaf.hpp"

#include "../../Src/Utils/FPoint.hpp"

#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Adaptive/FAdaptChebSymKernel.hpp"
43 44
#include "Kernels/Uniform/FUnifCell.hpp"
#include "Adaptive/FAdaptUnifKernel.hpp"
45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87

#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"

#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Core/FFmmAlgorithmThread.hpp"
#include "../../Src/Core/FFmmAlgorithmTask.hpp"

#include "../../Src/Components/FBasicKernels.hpp"

#include "../../Src/Files/FRandomLoader.hpp"

#include "../../Src/Adaptive/FAdaptiveCell.hpp"
#include "../../Src/Adaptive/FAdaptiveKernelWrapper.hpp"
#include "../../Src/Adaptive/FAbstractAdaptiveKernel.hpp"
#include "../../Src/Adaptive/FAdaptiveTestKernel.hpp"

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

#include "../../Src/Kernels/P2P/FP2PR.hpp"


/** This program show an example of use of the fmm basic algo
  * it also check that each particles is impacted each other particles
  */

// Simply create particles and try the kernels
int main(int argc, char ** argv){
    const FParameterNames LocalOptionMinMultipoleThreshod {
        {"-sM"},
        " s_min^M threshold for Multipole (l+1)^2 for Spherical harmonic."
    };
    const FParameterNames LocalOptionMinLocalThreshod {
        {"-SL"},
        " s_min^L threshold for Local  (l+1)^2 for Spherical harmonics."
    };

    FHelpDescribeAndExit(argc, argv,
                         "Test the adaptive FMM.",
                         FParameterDefinitions::NbParticles, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight,LocalOptionMinMultipoleThreshod,
                         LocalOptionMinLocalThreshod);


88
    typedef double FReal;
89
    const unsigned int P = 5 ;
90
    typedef FChebCell<FReal,P>                                        CellClass;
91
    //typedef FUnifCell<FReal,P>                                        CellClass;
92

93 94 95 96
    typedef FP2PParticleContainerIndexed<FReal>            ContainerClass;
    typedef FSimpleIndexedLeaf<FReal,ContainerClass>    LeafClass;
    typedef FInterpMatrixKernelR<FReal>                               MatrixKernelClass;
    typedef FAdaptiveChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,P> KernelClass;
97
    //typedef FAdaptiveUnifKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,P> KernelClass;
98 99
    typedef FAdaptiveCell< CellClass, ContainerClass >                                        CellWrapperClass;
    typedef FAdaptiveKernelWrapper< KernelClass, CellClass, ContainerClass >   KernelWrapperClass;
100
    typedef FOctree< FReal, CellWrapperClass, ContainerClass , LeafClass >                  OctreeClass;
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
    typedef FFmmAlgorithm<OctreeClass, CellWrapperClass, ContainerClass, KernelWrapperClass, LeafClass >     FmmClass;

    ///////////////////////What we do/////////////////////////////
    std::cout << ">> This executable has to be used to test the FMM algorithm.\n";
    //////////////////////////////////////////////////////////////

    const int NbLevels      = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 7);
    const int SizeSubLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 3);
    const int sminM    = FParameters::getValue(argc,argv,LocalOptionMinMultipoleThreshod.options, P*P*P);
    const int sminL     = FParameters::getValue(argc,argv,LocalOptionMinLocalThreshod.options, P*P*P);
    FTic counter;

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    const FReal boxWidth = 1.0;
117
    const FPoint<FReal> boxCenter(0.0, 0.0, 0.0);
118 119 120 121 122 123 124
    OctreeClass tree(NbLevels, SizeSubLevels, boxWidth, boxCenter);

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    struct Particle{
        int idx;
125
        FPoint<FReal> pos;
126 127 128 129 130 131 132 133 134
        FReal physicalValue;
        FReal forces[3];
        FReal pot;
    };

    Particle particles[3];
    memset(particles, 0, sizeof(Particle) * 3);
    {
        particles[0].idx = 0;
135
        particles[0].pos = FPoint<FReal>(0.1-0.5, 0.1-0.5, 0.1-0.5);
136 137
        particles[0].physicalValue = 1.0;
        particles[1].idx = 1;
138
        particles[1].pos = FPoint<FReal>(0.1-0.5+0.0625, 0.1-0.5+0.0625, 0.1-0.5+0.0625);
139 140
        particles[1].physicalValue = 1.0;
        particles[2].idx = 2;
141
        particles[2].pos = FPoint<FReal>(0.5-0.1, 0.5-0.1, 0.5-0.1);
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
        particles[2].physicalValue = 1.0;

        tree.insert(particles[0].pos, 0, particles[0].physicalValue);
        tree.insert(particles[1].pos, 1, particles[1].physicalValue);
        tree.insert(particles[2].pos, 2, particles[2].physicalValue);
    }

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    std::cout << "Working on particles ..." << std::endl;
    counter.tic();

    const MatrixKernelClass MatrixKernel;
    KernelWrapperClass kernels(NbLevels, boxWidth, boxCenter,&MatrixKernel,sminM,sminL);
    FmmClass algo(&tree,&kernels);  //FFmmAlgorithm FFmmAlgorithmThread
    algo.execute();

    counter.tac();
    std::cout << "Done  " << "(@Algorithm = " << counter.elapsed() << "s)." << std::endl;

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    for(int idxTarget = 0 ; idxTarget < 3 ; ++idxTarget){
        for(int idxSource = idxTarget+1 ; idxSource < 3 ; ++idxSource){
            FP2PR::MutualParticles(
                particles[idxTarget].pos.getX(),particles[idxTarget].pos.getY(), particles[idxTarget].pos.getZ(),
                particles[idxTarget].physicalValue, &particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
171
                &particles[idxTarget].forces[2], &particles[idxTarget].pot,
172 173
                particles[idxSource].pos.getX(),particles[idxSource].pos.getY(), particles[idxSource].pos.getZ(),
                particles[idxSource].physicalValue, &particles[idxSource].forces[0], &particles[idxSource].forces[1],
174
                &particles[idxSource].forces[2], &particles[idxSource].pot);
175 176 177 178 179 180 181
        }
    }

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    tree.forEachLeaf([&](LeafClass* leaf){
182
        const FSize nbParticles = leaf->getSrc()->getNbParticles();
183 184 185 186 187 188 189
        const FReal* posx = leaf->getSrc()->getPositions()[0];
        const FReal* posy = leaf->getSrc()->getPositions()[1];
        const FReal* posz = leaf->getSrc()->getPositions()[2];
        const FReal* fx = leaf->getSrc()->getForcesX();
        const FReal* fy = leaf->getSrc()->getForcesY();
        const FReal* fz = leaf->getSrc()->getForcesZ();
        const FReal* pot = leaf->getSrc()->getPotentials();
190 191
        const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
        for(FSize idxPart = 0 ; idxPart < nbParticles ; ++idxPart){
192 193 194 195 196 197 198 199 200 201 202 203 204 205
            std::cout << "[FMM] Particle pos " << posx[idxPart] << " " << posy[idxPart] << " " << posz[idxPart] << "\n";
            std::cout << "\t>> res pot " << pot[idxPart] << " forces " << fx[idxPart] << " " << fy[idxPart] << " " << fz[idxPart] << "\n";
            std::cout << "[Direct] Particle pos " << particles[indexes[idxPart]].pos.getX() << " " << particles[indexes[idxPart]].pos.getY() << " " << particles[indexes[idxPart]].pos.getZ() << "\n";
            std::cout << "\t>> res pot " << particles[indexes[idxPart]].pot << " forces " << particles[indexes[idxPart]].forces[0] << " " << particles[indexes[idxPart]].forces[1] << " " << particles[indexes[idxPart]].forces[2] << "\n";
            std::cout << "\n";
        }
    });

    return 0;
}