testAdaptiveChebSymFMMDebug.cpp 18.3 KB
Newer Older
COULAUD Olivier's avatar
COULAUD Olivier committed
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
COULAUD Olivier's avatar
COULAUD Olivier committed
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.
COULAUD Olivier's avatar
COULAUD Olivier committed
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".
COULAUD Olivier's avatar
COULAUD Olivier committed
19 20 21 22 23
// ===================================================================================

// ==== CMAKE =====
// @FUSE_BLAS
// ================
24 25 26
// Keep in private GIT
// @SCALFMM_PRIVATE

COULAUD Olivier's avatar
COULAUD Olivier committed
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54

#include <iostream>
#include <cstdio>


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

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

//#include "Components/FSimpleLeaf.hpp"

#include "Utils/FPoint.hpp"

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

#include "Components/FBasicKernels.hpp"
#include "Components/FSimpleIndexedLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"

#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebCell.hpp"

#include "Adaptive/FAdaptiveCell.hpp"
#include "Adaptive/FAdaptiveKernelWrapper.hpp"
#include "Adaptive/FAbstractAdaptiveKernel.hpp"
55
#include "Adaptive/FAdaptChebSymKernel.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
56 57 58 59 60 61
//
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Adaptive/FAdaptTools.hpp"
//#include "AdaptiveTree/FAdaptChebSymKernel.hpp"

62 63
#include "Adaptive/FAdaptivePrintKernel.hpp"

COULAUD Olivier's avatar
COULAUD Olivier committed
64 65 66 67 68
//
//
#include "Core/FFmmAlgorithm.hpp"
//#include "Core/FFmmAlgorithmThread.hpp"
//#include "Core/FFmmAlgorithmTask.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
69
#include "Utils/FParameterNames.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
70 71 72 73 74

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

BRAMAS Berenger's avatar
BRAMAS Berenger committed
75

COULAUD Olivier's avatar
COULAUD Olivier committed
76 77
// Simply create particles and try the kernels
int main(int argc, char ** argv){
78
    //
BRAMAS Berenger's avatar
BRAMAS Berenger committed
79 80
    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."};
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

    FHelpDescribeAndExit(argc, argv,
            "Test Adaptive kernel and compare it with the direct computation.",
            FParameterDefinitions::OctreeHeight,FParameterDefinitions::NbThreads,
            FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
            LocalOptionMinMultipoleThreshod,LocalOptionMinLocalThreshod);

    for (int i = 0 ; i< argc ; ++i){
        std::cout << argv[i] << "  " ;
    }
    std::cout << std::endl<< std::endl;
    // ---------------------------------------------
    const std::string fileName(FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options,   "../Data/noDistprolate50.out.fma"));
    const unsigned int TreeHeight      = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 3);
    const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
    //
    // accuracy
98
    typedef double FReal;
99
    const unsigned int P = 5 ;
100 101 102 103 104
    typedef FChebCell<FReal,P>                                        CellClass;
    typedef FP2PParticleContainerIndexed<FReal>            ContainerClass;
    typedef FSimpleIndexedLeaf<FReal,ContainerClass>    LeafClass;
    typedef FInterpMatrixKernelR<FReal>                               MatrixKernelClass;
    typedef FAdaptiveChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,P> KernelClass;
105 106
    typedef FAdaptiveCell< CellClass, ContainerClass >                                        CellWrapperClass;
    typedef FAdaptiveKernelWrapper< KernelClass, CellClass, ContainerClass >   KernelWrapperClass;
107
    typedef FOctree< FReal, CellWrapperClass, ContainerClass , LeafClass >                  OctreeClass;
108 109 110
    // FFmmAlgorithmTask FFmmAlgorithmThread
    typedef FFmmAlgorithm<OctreeClass, CellWrapperClass, ContainerClass, KernelWrapperClass, LeafClass >     FmmClass;

111 112 113 114 115
    typedef FP2PParticleContainerIndexed<FReal> ContainerClassDebug;
    typedef FSimpleLeaf<FReal, ContainerClass >  LeafClassDebug;
    typedef FChebCell<FReal,P> CellClassDebug;
    typedef FOctree<FReal, CellClassDebug,ContainerClassDebug,LeafClassDebug> OctreeClassDebug;
    typedef FChebSymKernel<FReal,CellClassDebug,ContainerClassDebug,MatrixKernelClass,P> KernelClassDebug;
116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
    typedef FFmmAlgorithm<OctreeClassDebug,CellClassDebug,ContainerClassDebug,KernelClassDebug,LeafClassDebug> FmmClassDebug;

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

    //////////////////////////////////////////////////////////////////////////////////
    // Not Random Loader
    //////////////////////////////////////////////////////////////////////////////////
131
    FFmaGenericLoader<FReal> loader(fileName);
132
    const FSize NbPart = loader.getNumberOfParticles() ;
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
    //////////////////////////////////////////////////////////////////////////////////

    OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
    OctreeClassDebug treedebug(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());

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

    std::cout << "Creating & Inserting " << NbPart << " particles ..." << std::endl;
    std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
    std::cout 		<< "         criteria SM:  "<< sminM     <<std::endl
            << "         criteria SL:  "<< sminL     <<std::endl <<std::endl;
    //

    counter.tic();
    FReal L= loader.getBoxWidth();
    //
150
    FmaRWParticle<FReal, 8,8>* const particles = new FmaRWParticle<FReal, 8,8>[NbPart];
151

152
    FPoint<FReal> minPos(L,L,L), maxPos(-L,-L,-L);
153 154 155
    //
    loader.fillParticle(particles,NbPart);

156
    for(FSize idxPart = 0 ; idxPart < NbPart; ++idxPart){
157
        const FPoint<FReal> PP(particles[idxPart].getPosition() ) ;
158 159 160 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 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202
        //
        minPos.setX(FMath::Min(minPos.getX(),PP.getX())) ;
        minPos.setY(FMath::Min(minPos.getY(),PP.getY())) ;
        minPos.setZ(FMath::Min(minPos.getZ(),PP.getZ())) ;
        maxPos.setX(FMath::Max(maxPos.getX(),PP.getX())) ;
        maxPos.setY(FMath::Max(maxPos.getY(),PP.getY())) ;
        maxPos.setZ(FMath::Max(maxPos.getZ(),PP.getZ())) ;
        //
        tree.insert(PP, idxPart, particles[idxPart].getPhysicalValue());
        treedebug.insert(PP, idxPart, particles[idxPart].getPhysicalValue());

    }

    counter.tac();
    std::cout << "Data are inside the box delimited by "<<std::endl
            << "         Min corner:  "<< minPos<<std::endl
            << "         Max corner:  "<< maxPos<<std::endl <<std::endl;
    std::cout << "Done  " << "(@Creating and Inserting Particles = " << counter.elapsed() << " s)." << std::endl;

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

    std::cout << "Working on particles ..." << std::endl;
    // For debug purpose
    //  Set Global id
    //
    long int idCell  = setGlobalID(tree);
    counter.tic();
    const MatrixKernelClass MatrixKernel;
    KernelWrapperClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel,sminM,sminL);            // FTestKernels FBasicKernels
    FmmClass algo(&tree,&kernels);  //FFmmAlgorithm FFmmAlgorithmThread
    algo.execute();

    KernelClassDebug kernelsdebug(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);            // FTestKernels FBasicKernels
    FmmClassDebug algodebug(&treedebug,&kernelsdebug);  //FFmmAlgorithm FFmmAlgorithmThread
    algodebug.execute();

    counter.tac();
    std::cout << "Done  " << "(@Algorithm = " << counter.elapsed() << " s)." << std::endl;
    //
        FReal energyD = 0.0 ;
    /////////////////////////////////////////////////////////////////////////////////////////////////
    // Compute direct energy
    /////////////////////////////////////////////////////////////////////////////////////////////////

203
    for(FSize idx = 0 ; idx < loader.getNumberOfParticles()  ; ++idx){
204 205 206 207 208 209
        energyD +=  particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
    }
    /////////////////////////////////////////////////////////////////////////////////////////////////
    // Compare
    /////////////////////////////////////////////////////////////////////////////////////////////////
    {
210 211
        FMath::FAccurater<FReal> potentialDiff;
        FMath::FAccurater<FReal> fx, fy, fz;
212 213 214 215 216
        FReal energy= 0.0;
        { // Check that each particle has been summed with all other

            //    std::cout << "indexPartOrig || DIRECT V fx || FMM V fx" << std::endl;

217
            tree.forEachCellLeaf([&](CellWrapperClass* cell, LeafClass* leaf){
218 219
                FMath::FAccurater<FReal> potentialDiffLeaf;
                FMath::FAccurater<FReal> fxLeaf, fyLeaf, fzLeaf;
220

221 222 223 224 225
                const FReal*const potentials        = leaf->getTargets()->getPotentials();
                const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
                const FReal*const forcesX            = leaf->getTargets()->getForcesX();
                const FReal*const forcesY            = leaf->getTargets()->getForcesY();
                const FReal*const forcesZ            = leaf->getTargets()->getForcesZ();
226 227
                const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
228

229 230
                for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const FSize indexPartOrig = indexes[idxPart];
231 232 233 234 235 236 237 238 239 240 241
                    potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
                    fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
                    fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
                    fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
                    energy   += potentials[idxPart]*physicalValues[idxPart];

                    //          std::cout << indexPartOrig
                    //                    << " " << particles[indexPartOrig].getPotential() << " " << particles[indexPartOrig].getForces()[0]
                    //                    << " " << potentials[idxPart] << " " << forcesX[idxPart]
                    //                    << std::endl;

242 243 244 245
                    potentialDiffLeaf.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
                    fxLeaf.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
                    fyLeaf.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
                    fzLeaf.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
246
                }
247 248 249 250 251 252

                std::cout << "Cell id " << cell->getGlobalId( ) << " morton " << cell->getMortonIndex() << "\n";
                std::cout << "\t Potential " << potentialDiffLeaf << std::endl;
                std::cout << "\t Fx " << fxLeaf << std::endl;
                std::cout << "\t Fy " << fyLeaf << std::endl;
                std::cout << "\t Fz " << fzLeaf << std::endl;
253 254 255 256 257 258 259 260 261 262 263
            });
        }

        // Print for information
        std::cout << "Energy [relative L2 error] "  << FMath::Abs(energy-energyD) /energyD << std::endl;
        std::cout << "Potential " << potentialDiff << std::endl;
        std::cout << "Fx " << fx << std::endl;
        std::cout << "Fy " << fy << std::endl;
        std::cout << "Fz " << fz << std::endl;
    }
    {
264 265
        FMath::FAccurater<FReal> potentialDiff;
        FMath::FAccurater<FReal> fx, fy, fz;
266 267 268 269 270 271 272 273 274 275 276
        FReal energy= 0.0;
        { // Check that each particle has been summed with all other

            //    std::cout << "indexPartOrig || DIRECT V fx || FMM V fx" << std::endl;

            treedebug.forEachLeaf([&](LeafClassDebug* leaf){
                const FReal*const potentials        = leaf->getTargets()->getPotentials();
                const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
                const FReal*const forcesX            = leaf->getTargets()->getForcesX();
                const FReal*const forcesY            = leaf->getTargets()->getForcesY();
                const FReal*const forcesZ            = leaf->getTargets()->getForcesZ();
277 278
                const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
279

280 281
                for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                    const FSize indexPartOrig = indexes[idxPart];
282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314
                    potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
                    fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
                    fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
                    fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
                    energy   += potentials[idxPart]*physicalValues[idxPart];

                    //          std::cout << indexPartOrig
                    //                    << " " << particles[indexPartOrig].getPotential() << " " << particles[indexPartOrig].getForces()[0]
                    //                    << " " << potentials[idxPart] << " " << forcesX[idxPart]
                    //                    << std::endl;

                }
            });
        }

        // Print for information
        std::cout << "Energy [relative L2 error] "  << FMath::Abs(energy-energyD) /energyD << std::endl;
        std::cout << "Potential " << potentialDiff << std::endl;
        std::cout << "Fx " << fx << std::endl;
        std::cout << "Fy " << fy << std::endl;
        std::cout << "Fz " << fz << std::endl;
    }

    {

        tree.forEachCellWithLevel([&](CellWrapperClass* cell, int level){
            if(cell->hasDevelopment()){
                FAssertLF(cell->isAdaptive());
                CellClass* testcell = cell->getRealCell();
                FAssertLF(testcell);
                CellClassDebug* testcelltest = treedebug.getCell(testcell->getMortonIndex(), level);
                FAssertLF(testcelltest);

315
                FMath::FAccurater<FReal> pole, local;
316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347

                for(int idx = 0 ; idx < testcell->getVectorSize() ; ++idx){
                    pole.add(testcelltest->getMultipole(0)[idx], testcell->getMultipole(0)[idx]);
                    local.add(testcelltest->getLocal(0)[idx], testcell->getLocal(0)[idx]);
                }

                if(pole.getRelativeInfNorm() > 1e-6 || pole.getRelativeL2Norm() > 1e-6){
                    std::cout << " id " << cell->getGlobalId( ) << " pole has error at level " << level << " morton idx " << testcell->getMortonIndex() << " err " << pole << "\n";
                }
                if(local.getRelativeInfNorm() > 1e-6 || local.getRelativeL2Norm() > 1e-6){
                    std::cout << " id " << cell->getGlobalId( ) << " local has error at level " << level << " morton idx " << testcell->getMortonIndex() << " err " << local << "\n";
                }
            }
        });
    }

    delete[] particles;


    OctreeClass::Iterator octreeIterator(&tree);
    std::ofstream file("aa.tree", std::ofstream::out );

    //
    ////////////////////////////////////////////////////////////////////
    //              Export adaptive tree in our format
    ////////////////////////////////////////////////////////////////////
    //
    // -----------------------------------------------------
    //
    //
    //  Set Global id
    //
COULAUD Olivier's avatar
COULAUD Olivier committed
348
//	long int idCell  = setGlobalID(tree);
349 350 351 352 353 354 355 356 357 358 359 360 361 362
    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    tree.forEachCellLeaf([&](CellWrapperClass* cell, LeafClass* leaf){
        file << "Cell Id " << cell->getGlobalId( ) << " Nb particles "<<  leaf->getSrc()->getNbParticles()<<std::endl;
    });

    octreeIterator.gotoTop() ;  // here we are at level 1 (first child)
    //	octreeIterator.moveDown() ;
    octreeIterator.gotoLeft();
    //	octreeIterator.moveDown() ; // We are at the levell 2
    std::cout << " Number of Cells: " << idCell <<std::endl;
    //
    std::cout << "Top of the octree " << octreeIterator.level() << std::endl ;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
363
    for(unsigned int idxLevel = 1; idxLevel < TreeHeight ;  ++idxLevel){
364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383
        file << std::endl << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<< std::endl;
        file << "  Level " << idxLevel <<"  Level  "<<  octreeIterator.level()<<  "  -- leave level " <<   std::boolalpha <<  octreeIterator.isAtLeafLevel() << std::endl;
        do{
            if(octreeIterator.getCurrentCell()->hasDevelopment()){
                file <<"Cell id  "<< octreeIterator.getCurrentCell()->getGlobalId( ) << "   "<<*(octreeIterator.getCurrentCell())<< std::endl ;
            }
        } while(octreeIterator.moveRight());
        octreeIterator.moveDown() ;
        octreeIterator.gotoLeft();
    }
    std::cout << "   END    " << std::endl;

    // Check
    octreeIterator.gotoBottomLeft();
    do {
        std::cout << " Level " <<octreeIterator.level() <<std::endl;
    }while(octreeIterator.moveUp() );
    std::cout << "   RETURN 0  " << std::endl;

    return 0;
COULAUD Olivier's avatar
COULAUD Olivier committed
384 385 386 387 388
}