testFmmAlgorithmProcRotation.cpp 18 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

34 35 36 37 38 39
#include <iostream>

#include <cstdio>
#include <cstdlib>

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

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

47
typedef double FReal;
48 49 50

#ifdef VALIDATE_FMM

51
static const FReal Epsilon = FReal(0.00005);
52 53 54 55 56 57 58

///////////////////////////////////////////////////////
// to test equality between good and potentialy bad solution
///////////////////////////////////////////////////////
/** To compare data */
template <class CellClass>
bool isEqualPole(const CellClass& me, const CellClass& other, FReal*const cumul){
59
    FMath::FAccurater<FReal> accurate;
60 61 62 63 64 65
    for(int idx = 0; idx < me.getArraySize(); ++idx){
        accurate.add(me.getMultipole()[idx].getImag(),other.getMultipole()[idx].getImag());
        accurate.add(me.getMultipole()[idx].getReal(),other.getMultipole()[idx].getReal());
    }
    *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 FRotationCell<7>& me, const FRotationCell<7>& other, FReal*const cumul){
70
    FMath::FAccurater<FReal> accurate;
71 72 73 74 75 76
    for(int idx = 0; idx < me.getArraySize(); ++idx){
        accurate.add(me.getLocal()[idx].getImag(),other.getLocal()[idx].getImag());
        accurate.add(me.getLocal()[idx].getReal(),other.getLocal()[idx].getReal());
    }
    *cumul = accurate.getInfNorm()+ accurate.getL2Norm();
    return accurate.getInfNorm() < Epsilon && accurate.getL2Norm() < Epsilon;//FMath::LookEqual(cumul,FReal(0.0));
77 78 79 80 81
}


template<class OctreeClass, class ContainerClass>
void ValidateFMMAlgoProc(OctreeClass* const badTree,
82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
                         OctreeClass* const valideTree,
                         FMpi * const app){
    std::cout << "Check Result\n";
    {
        const int OctreeHeight = valideTree->getHeight();
        typename OctreeClass::Iterator octreeIterator(badTree);
        octreeIterator.gotoBottomLeft();

        typename OctreeClass::Iterator octreeIteratorValide(valideTree);
        octreeIteratorValide.gotoBottomLeft();

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

            do {
                if(octreeIterator.getCurrentGlobalIndex() != octreeIteratorValide.getCurrentGlobalIndex()){
                    std::cout << "Error index are not equal!" << std::endl;
                }
                else{
                    FReal cumul;
                    if( !isEqualPole(*octreeIterator.getCurrentCell(),*octreeIteratorValide.getCurrentCell(),&cumul) ){
                        std::cout << "Pole Data are different for proc "<< app->global().processId() <<" Cumul " << cumul << " at level " << level << " index is " << octreeIterator.getCurrentGlobalIndex() << std::endl;
                    }
                    if( !isEqualLocal(*octreeIterator.getCurrentCell(),*octreeIteratorValide.getCurrentCell(),&cumul) ){
                        std::cout << "Local Data are different for proc "<< app->global().processId() <<"  Cumul " << cumul << " at level " << level << " index is " << octreeIterator.getCurrentGlobalIndex() << std::endl;
                    }
                }

            } while(octreeIterator.moveRight() && octreeIteratorValide.moveRight());

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

            octreeIteratorValide.moveUp();
            octreeIteratorValide.gotoLeft();
        }
120
    }
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
    {
        // Check that each particle has been summed with all other
        typename OctreeClass::Iterator octreeIterator(badTree);
        octreeIterator.gotoBottomLeft();

        typename OctreeClass::Iterator octreeIteratorValide(valideTree);
        octreeIteratorValide.gotoBottomLeft();

        while(octreeIteratorValide.getCurrentGlobalIndex() != octreeIterator.getCurrentGlobalIndex()){
            octreeIteratorValide.moveRight();
        }

        do {

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

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

160
            for(FSize idxLeaf = 0 ; idxLeaf < firstLeaf->getNbParticles() ; ++idxLeaf){
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187

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

                if( idxValideLeaf < valideLeaf->getNbParticles() ){
                    // If a particles has been impacted by less than NbPart - 1 (the current particle)
                    // there is a problem
                    bool error = false;
                    if( FMath::RelatifDiff(validePotentials[idxValideLeaf] , potentials[idxLeaf])  > Epsilon ){
                        std::cout << " Potential error : " << validePotentials[idxValideLeaf] << " " << potentials[idxLeaf] << "\n";
                        error = true;
                    }
                    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]
                                     << " y " << valideForcesY[idxValideLeaf]  << " " << forcesY[idxLeaf]
                                        << " z " << valideForcesZ[idxValideLeaf]  << " " << forcesZ[idxLeaf] << "\n";
                        error = true;
                    }
                    if( error ){
188 189
                        std::cout << "At position " << FPoint<FReal>(validePositionX[idxValideLeaf],validePositionY[idxValideLeaf],validePositionZ[idxValideLeaf])
                                  << " == " << FPoint<FReal>(positionX[idxLeaf],positionY[idxLeaf],positionZ[idxLeaf]) << std::endl;
190 191 192
                    }
                }
                else{
193
                    std::cout << "Particle not found " << FPoint<FReal>(positionX[idxLeaf],positionY[idxLeaf],positionZ[idxLeaf]) << std::endl;
194 195 196 197
                }
            }

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

200
    std::cout << "Done\n";
201 202 203 204 205 206
}
#endif


// Simply create particles and try the kernels
int main(int argc, char ** argv){
207 208 209 210
    FHelpDescribeAndExit(argc, argv,
                         "Test FMM distributed algorithm with the rotation kernel.",
                         FParameterDefinitions::OctreeHeight, FParameterDefinitions::OctreeSubHeight,
                         FParameterDefinitions::InputFile);
211

212
    const int P = 7;
213
    // For Rotation test ::
214
    typedef FRotationCell<FReal,P>         CellClass;
215
    typedef FP2PParticleContainer<FReal>         ContainerClass;
216

217 218
    typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
    typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
219
    typedef FRotationKernel< FReal, CellClass, ContainerClass,P >     KernelClass;
220

221
    typedef FFmmAlgorithmThreadProc<OctreeClass,  CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
222
#ifdef VALIDATE_FMM
223
    typedef FFmmAlgorithmThread<OctreeClass,  CellClass, ContainerClass, KernelClass, LeafClass > FmmClassNoProc;
224
#endif
225

226 227 228
    ///////////////////////What we do/////////////////////////////
    std::cout << ">> This executable has to be used to test Proc Rotation algorithm.\n";
    //////////////////////////////////////////////////////////////
229

230
    FMpi app( argc, argv);
231

232 233 234 235 236 237 238
    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);
239

240
    std::cout << "Opening : " << filename << "\n";
241

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

248
    // ----Modified For Rotation----------------------------
249

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

252
    // -----------------------------------------------------
253

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

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

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

273 274
        TestParticle* particles = new TestParticle[loader.getMyNumberOfParticles()];
        memset(particles, 0, (unsigned int) (sizeof(TestParticle) * loader.getMyNumberOfParticles()));
275

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

        FVector<TestParticle> finalParticles;
        FLeafBalance balancer;
282

283
        FMpiTreeBuilder< FReal,TestParticle >::DistributeArrayToContainer(app.global(),particles,
284 285 286 287 288 289 290 291 292 293 294 295 296 297 298
                                                                    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);
        }

        delete[] particles;

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

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

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

313 314 315 316
    // -----------------------------------------------------
    std::cout << "Create kernel..." << std::endl;
    {
        KernelClass kernels( NbLevels,loader.getBoxWidth(), loader.getCenterOfBox());
317

318
        std::cout << "Done  " << " in " << counter.elapsed() << "s)." << std::endl;
319

320
        // -----------------------------------------------------
321

322 323 324
        std::cout << "Working on particles ..." << std::endl;
        {
            FmmClass algo(app.global(),&tree,&kernels);
325

326 327 328
            counter.tic();
            algo.execute(FFmmM2M);
            counter.tac();
329

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

332
            { // get sum forces&potential
333

334 335
                FReal potential = 0;
                FReal fx = 0.0, fy = 0.0, fz = 0.0;
336

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

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

352
                std::cout << "My potential is " << potential << std::endl;
353

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


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

#ifdef VALIDATE_FMM
367 368 369
            {
                OctreeClass treeValide(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
                {
370 371
                    FFmaGenericLoader<FReal> loaderSeq(filename);
                    FPoint<FReal> position;
372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393
                    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();
394
                    const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
395

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

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

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

413
    // -----------------------------------------------------
414

415
    return 0;
416
}