testSphericalProcAlgorithm.cpp 17.7 KB
Newer Older
1
// See LICENCE file at project root
2

3 4 5 6
// ==== CMAKE =====
// @FUSE_MPI
// ================

7 8 9 10
#include "Utils/FTic.hpp"
#include "Utils/FMpi.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FMath.hpp"
11

12 13
#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"
14

15 16
#include "Kernels/Spherical/FSphericalKernel.hpp"
#include "Kernels/Spherical/FSphericalCell.hpp"
17

18 19
#include "Kernels/Rotation/FRotationKernel.hpp"
#include "Kernels/Rotation/FRotationCell.hpp"
20

21 22
#include "Core/FFmmAlgorithmThreadProc.hpp"
#include "Core/FFmmAlgorithmThread.hpp"
23

24 25
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainer.hpp"
26

27 28
#include "Files/FMpiFmaGenericLoader.hpp"
#include "Files/FMpiTreeBuilder.hpp"
29

30
#include "Utils/FLeafBalance.hpp"
31

32
#include "Utils/FParameterNames.hpp"
33

berenger-bramas's avatar
berenger-bramas committed
34
#include <iostream>
35

36 37
#include <cstdio>
#include <cstdlib>
berenger-bramas's avatar
berenger-bramas committed
38

39
// Uncoment to validate the FMM
40
#define VALIDATE_FMM
41

42
/** This program show an example of use of
43
  * the fmm basic algo it also check that eachh particles is little or longer
44 45 46
  * related that each other
  */

47
typedef double FReal;
48

49 50 51 52
#ifdef VALIDATE_FMM

static const FReal Epsilon = FReal(0.0005);

53 54 55 56
///////////////////////////////////////////////////////
// to test equality between good and potentialy bad solution
///////////////////////////////////////////////////////
/** To compare data */
57 58
template <class CellClass>
bool isEqualPole(const CellClass& me, const CellClass& other, FReal*const cumul){
59
    FMath::FAccurater<FReal> accurate;
60
    for(int idx = 0; idx < CellClass::multipole_t::Size; ++idx){
61 62
        accurate.add(me.getMultipoleData().get()[idx].imag(),other.getMultipoleData().get()[idx].imag());
        accurate.add(me.getMultipoleData().get()[idx].real(),other.getMultipoleData().get()[idx].real());
63
    }
64 65
    *cumul = accurate.getInfNorm()+ accurate.getL2Norm();
    return accurate.getInfNorm() < Epsilon && accurate.getL2Norm() < Epsilon;//FMath::LookEqual(cumul,FReal(0.0));
66
}
67

68
/** To compare data */
69
bool isEqualLocal(const FSphericalCell<FReal>& me, const FSphericalCell<FReal>& other,FReal*const cumul){
70
    FMath::FAccurater<FReal> accurate;
71
    for(int idx = 0; idx < FSphericalCell<FReal>::local_expansion_t::Size; ++idx){
72 73
        accurate.add(me.getLocalExpansionData().get()[idx].imag(),other.getLocalExpansionData().get()[idx].imag());
        accurate.add(me.getLocalExpansionData().get()[idx].real(),other.getLocalExpansionData().get()[idx].real());
74
    }
75 76
    *cumul = accurate.getInfNorm()+ accurate.getL2Norm();
    return accurate.getInfNorm() < Epsilon && accurate.getL2Norm() < Epsilon;//FMath::LookEqual(cumul,FReal(0.0));
77
}
78

79

80
template<class OctreeClass, class ContainerClass>
berenger-bramas's avatar
berenger-bramas committed
81
void ValidateFMMAlgoProc(OctreeClass* const badTree,
82
                         OctreeClass* const valideTree){
83 84
    std::cout << "Check Result\n";
    {
85
        const int OctreeHeight = valideTree->getHeight();
86
        typename OctreeClass::Iterator octreeIterator(badTree);
87 88
        octreeIterator.gotoBottomLeft();

89
        typename OctreeClass::Iterator octreeIteratorValide(valideTree);
90 91
        octreeIteratorValide.gotoBottomLeft();

92
        for(int level = OctreeHeight - 1 ; level > 1 ; --level){
93
            while(octreeIteratorValide.getCurrentGlobalIndex() != octreeIterator.getCurrentGlobalIndex()){
94 95 96
                octreeIteratorValide.moveRight();
            }

97
            do {
98 99 100 101 102
                if(octreeIterator.getCurrentGlobalIndex() != octreeIteratorValide.getCurrentGlobalIndex()){
                    std::cout << "Error index are not equal!" << std::endl;
                }
                else{
                    FReal cumul;
103
                    if( !isEqualPole(*octreeIterator.getCurrentCell(),*octreeIteratorValide.getCurrentCell(),&cumul) ){
104
                        std::cout << "Pole Data are different. Cumul " << cumul << " at level " << level << " index is " << octreeIterator.getCurrentGlobalIndex() << std::endl;
105
                    }
106
                    if( !isEqualLocal(*octreeIterator.getCurrentCell(),*octreeIteratorValide.getCurrentCell(),&cumul) ){
107
                        std::cout << "Local Data are different. Cumul " << cumul << " at level " << level << " index is " << octreeIterator.getCurrentGlobalIndex() << std::endl;
108 109 110
                    }
                }

111
            } while(octreeIterator.moveRight() && octreeIteratorValide.moveRight());
112 113 114 115 116 117 118 119 120

            octreeIterator.moveUp();
            octreeIterator.gotoLeft();

            octreeIteratorValide.moveUp();
            octreeIteratorValide.gotoLeft();
        }
    }
    {
121
        // Check that each particle has been summed with all other
122
        typename OctreeClass::Iterator octreeIterator(badTree);
123
        octreeIterator.gotoBottomLeft();
124

125
        typename OctreeClass::Iterator octreeIteratorValide(valideTree);
126
        octreeIteratorValide.gotoBottomLeft();
127

128 129 130
        while(octreeIteratorValide.getCurrentGlobalIndex() != octreeIterator.getCurrentGlobalIndex()){
            octreeIteratorValide.moveRight();
        }
131

132
        do {
133 134

            if( octreeIterator.getCurrentListSrc()->getNbParticles() != octreeIteratorValide.getCurrentListSrc()->getNbParticles()){
135 136 137 138 139
                std::cout << " Particules numbers is different " << std::endl;
            }
            if( octreeIterator.getCurrentGlobalIndex() != octreeIteratorValide.getCurrentGlobalIndex()){
                std::cout << " Index are differents " << std::endl;
            }
140

141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158
            ContainerClass* firstLeaf = octreeIterator.getCurrentListTargets();
            ContainerClass* valideLeaf = octreeIteratorValide.getCurrentListTargets();

            const FReal*const potentials = firstLeaf->getPotentials();
            const FReal*const forcesX = firstLeaf->getForcesX();
            const FReal*const forcesY = firstLeaf->getForcesY();
            const FReal*const forcesZ = firstLeaf->getForcesZ();
            const FReal*const positionX = firstLeaf->getPositions()[0];
            const FReal*const positionY = firstLeaf->getPositions()[1];
            const FReal*const positionZ = firstLeaf->getPositions()[2];
            const FReal*const validePositionX = valideLeaf->getPositions()[0];
            const FReal*const validePositionY = valideLeaf->getPositions()[1];
            const FReal*const validePositionZ = valideLeaf->getPositions()[2];
            const FReal*const validePotentials = valideLeaf->getPotentials();
            const FReal*const valideForcesX = valideLeaf->getForcesX();
            const FReal*const valideForcesY = valideLeaf->getForcesY();
            const FReal*const valideForcesZ = valideLeaf->getForcesZ();

159
            for(FSize idxLeaf = 0 ; idxLeaf < firstLeaf->getNbParticles() ; ++idxLeaf){
160 161 162 163

                int idxValideLeaf = 0;
                for(; idxValideLeaf < valideLeaf->getNbParticles() ; ++idxValideLeaf){
                    if( FMath::LookEqual(validePositionX[idxValideLeaf],positionX[idxLeaf]) &&
164 165
                            FMath::LookEqual(validePositionY[idxValideLeaf],positionY[idxLeaf]) &&
                            FMath::LookEqual(validePositionZ[idxValideLeaf],positionZ[idxLeaf]) ){
166 167
                        break;
                    }
168
                }
169

170
                if( idxValideLeaf < valideLeaf->getNbParticles() ){
171 172 173
                    // If a particles has been impacted by less than NbPart - 1 (the current particle)
                    // there is a problem
                    bool error = false;
174 175
                    if( FMath::RelatifDiff(validePotentials[idxValideLeaf] , potentials[idxLeaf])  > Epsilon ){
                        std::cout << " Potential error : " << validePotentials[idxValideLeaf] << " " << potentials[idxLeaf] << "\n";
176 177
                        error = true;
                    }
178 179 180 181
                    if( FMath::RelatifDiff(valideForcesX[idxValideLeaf],forcesX[idxLeaf]) > Epsilon
                            || FMath::RelatifDiff(valideForcesY[idxValideLeaf],forcesY[idxLeaf]) > Epsilon
                            || FMath::RelatifDiff(valideForcesZ[idxValideLeaf],forcesZ[idxLeaf]) > Epsilon){
                        std::cout << " Forces error : x " << valideForcesX[idxValideLeaf] << " " << forcesX[idxLeaf]
182 183
                                     << " y " << valideForcesY[idxValideLeaf]  << " " << forcesY[idxLeaf]
                                        << " z " << valideForcesZ[idxValideLeaf]  << " " << forcesZ[idxLeaf] << "\n";
184 185 186
                        error = true;
                    }
                    if( error ){
187 188
                        std::cout << "At position " << FPoint<FReal>(validePositionX[idxValideLeaf],validePositionY[idxValideLeaf],validePositionZ[idxValideLeaf])
                                  << " == " << FPoint<FReal>(positionX[idxLeaf],positionY[idxLeaf],positionZ[idxLeaf]) << std::endl;
189 190 191
                    }
                }
                else{
192
                    std::cout << "Particle not found " << FPoint<FReal>(positionX[idxLeaf],positionY[idxLeaf],positionZ[idxLeaf]) << std::endl;
193
                }
194
            }
195 196

        } while(octreeIterator.moveRight() && octreeIteratorValide.moveRight());
197 198 199 200
    }

    std::cout << "Done\n";
}
201
#endif
202

203

204 205
// Simply create particles and try the kernels
int main(int argc, char ** argv){
206 207 208 209 210 211
    FHelpDescribeAndExit(argc, argv,
                         "Run a Spherical Harmonic (Old Implementation) FMM kernel with mpir parallelization.\n"
                         "The input file should be in binary format to enable distributed access.",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::SHDevelopment);

212 213
    typedef FSphericalCell<FReal>         CellClass;
    typedef FP2PParticleContainer<FReal>         ContainerClass;
214

215 216 217
    typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
    typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
    typedef FSphericalKernel<FReal, CellClass, ContainerClass >     KernelClass;
218

219 220
    typedef FFmmAlgorithmThreadProc<OctreeClass,  CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
    typedef FFmmAlgorithmThread<OctreeClass,  CellClass, ContainerClass, KernelClass, LeafClass > FmmClassNoProc;
221 222


223 224 225
    ///////////////////////What we do/////////////////////////////
    std::cout << ">> This executable has to be used to test Spherical algorithm.\n";
    //////////////////////////////////////////////////////////////
226

227
    FMpi app( argc, argv);
228

229 230 231 232 233 234 235 236
    const int DevP = FParameters::getValue(argc,argv,FParameterDefinitions::SHDevelopment.options, 8);
    const int NbLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 5);
    const int SizeSubLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 3);
    FTic counter;
    const char* const defaultFilename = (sizeof(FReal) == sizeof(float))?
                "../Data/test20k.bin.fma.single":
                "../Data/test20k.bin.fma.double";
    const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFilename);
237

238 239
    std::cout << "Opening : " << filename << "\n";

240
    FMpiFmaGenericLoader<FReal> loader(filename, app.global());
241 242 243 244
    if(!loader.isOpen()){
        std::cout << "Loader Error, " << filename << " is missing\n";
        return 1;
    }
245

246
    CellClass::Init(DevP);
247

248

249
    OctreeClass tree(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
250

251
    // -----------------------------------------------------
252

253 254
    std::cout << "Creating & Inserting " << loader.getNumberOfParticles() << " particles ..." << std::endl;
    std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
255
    counter.tic();
256

257 258 259 260 261 262 263 264
    if( app.global().processCount() != 1){
        //////////////////////////////////////////////////////////////////////////////////
        // Build tree from mpi loader
        //////////////////////////////////////////////////////////////////////////////////
        std::cout << "Build Tree ..." << std::endl;
        counter.tic();

        struct TestParticle{
265
            FPoint<FReal> position;
266
            FReal physicalValue;
267
            const FPoint<FReal>& getPosition(){
268 269 270
                return position;
            }
        };
271

272 273
        TestParticle* particles = new TestParticle[loader.getNumberOfParticles()];
        memset(particles, 0, sizeof(TestParticle) * loader.getNumberOfParticles());
274

275
        for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
276 277
            loader.fillParticle(&particles[idxPart].position,&particles[idxPart].physicalValue);
        }
278

279 280
        FVector<TestParticle> finalParticles;
        FLeafBalance balancer;
281
        // FMpiTreeBuilder< FReal,TestParticle >::ArrayToTree(app.global(), particles, loader.getNumberOfParticles(),
282 283 284
        // 						 tree.getBoxCenter(),
        // 						 tree.getBoxWidth(),
        // 						 tree.getHeight(), &finalParticles,&balancer);
285
        FMpiTreeBuilder< FReal,TestParticle >::DistributeArrayToContainer(app.global(),particles,
286 287 288 289 290 291 292
                                                                    loader.getMyNumberOfParticles(),
                                                                    tree.getBoxCenter(),
                                                                    tree.getBoxWidth(),tree.getHeight(),
                                                                    &finalParticles, &balancer);

        for(int idx = 0 ; idx < finalParticles.getSize(); ++idx){
            tree.insert(finalParticles[idx].position,finalParticles[idx].physicalValue);
293

294
        }
295

296
        delete[] particles;
297

298 299 300 301 302 303
        counter.tac();
        std::cout << "Done  " << "(" << counter.elapsed() << "s)." << std::endl;

        //////////////////////////////////////////////////////////////////////////////////
    }
    else{
304
        FPoint<FReal> position;
305 306 307 308 309
        FReal physicalValue;
        for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
            loader.fillParticle(&position,&physicalValue);
            tree.insert(position, physicalValue);
        }
310
    }
311

312 313
    counter.tac();
    std::cout << "Done  " << "(@Creating and Inserting Particles = " << counter.elapsed() << "s)." << std::endl;
314

315 316
    // -----------------------------------------------------
    std::cout << "Create kernel..." << std::endl;
317

318
    KernelClass kernels(DevP, NbLevels,loader.getBoxWidth(), loader.getCenterOfBox());
319

320
    std::cout << "Done  " << " in " << counter.elapsed() << "s)." << std::endl;
321

322
    // -----------------------------------------------------
323

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

326
    FmmClass algo(app.global(),&tree,&kernels);
327

328 329 330
    counter.tic();
    algo.execute();
    counter.tac();
331

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

334
    { // get sum forces&potential
335

336 337
        FReal potential = 0;
        FReal fx = 0.0, fy = 0.0, fz = 0.0;
338

339 340 341 342 343
        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();
344
            const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
345

346
            for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
347 348 349 350 351 352
                potential += potentials[idxPart];
                fx += forcesX[idxPart];
                fy += forcesY[idxPart];
                fz += forcesZ[idxPart];
            }
        });
353

354
        std::cout << "My potential is " << potential << std::endl;
355

356 357 358 359
        potential = app.global().reduceSum(potential);
        fx = app.global().reduceSum(fx);
        fy = app.global().reduceSum(fy);
        fz = app.global().reduceSum(fz);
360

361

362 363 364 365
        if(app.global().processId() == 0){
            std::cout << "Foces Sum  x = " << fx << " y = " << fy << " z = " << fz << std::endl;
            std::cout << "Potential Sum = " << potential << std::endl;
        }
366 367
    }

368
#ifdef VALIDATE_FMM
369
    {
370 371
        OctreeClass treeValide(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
        {
372 373
            FFmaGenericLoader<FReal> loaderSeq(filename);
            FPoint<FReal> position;
374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395
            FReal physicalValue;
            for(FSize idxPart = 0 ; idxPart < loaderSeq.getNumberOfParticles() ; ++idxPart){
                loaderSeq.fillParticle(&position,&physicalValue);
                treeValide.insert(position,physicalValue);
            }
        }

        std::cout << "Working on particles ..." << std::endl;
        FmmClassNoProc algoValide(&treeValide,&kernels);
        counter.tic();
        algoValide.execute();
        counter.tac();
        std::cout << "Done  " << "(@Algorithm = " << counter.elapsed() << "s)." << std::endl;

        FReal potential = 0;
        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();
396
            const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
397

398
            for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
399 400 401 402 403 404 405 406 407 408 409
                potential += potentials[idxPart];
                fx += forcesX[idxPart];
                fy += forcesY[idxPart];
                fz += forcesZ[idxPart];
            }
        });

        std::cout << "Foces Sum  x = " << fx << " y = " << fy << " z = " << fz << std::endl;
        std::cout << "Potential = " << potential << std::endl;

        ValidateFMMAlgoProc<OctreeClass,ContainerClass>(&tree,&treeValide);
410
    }
411
#endif
412 413


414
    // -----------------------------------------------------
415

416
    return 0;
417
}