testBlockedChebyshev.cpp 10.1 KB
Newer Older
1 2
// ==== CMAKE =====
// @FUSE_BLAS
3
// 
4
// ================
5 6
// Keep in private GIT

7

8
#include "../../Src/Utils/FGlobal.hpp"
9
#include "../../Src/GroupTree/Core/FGroupTree.hpp"
10 11

#include "../../Src/Components/FSimpleLeaf.hpp"
12
#include "../../Src/Components/FSymbolicData.hpp"
13 14 15 16 17
#include "../../Src/Containers/FVector.hpp"

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

#include "../../Src/Kernels/Chebyshev/FChebSymKernel.hpp"
18
#include "../../Src/Kernels/Chebyshev/FChebCell.hpp"
19 20 21 22 23 24 25 26 27
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"

#include "../../Src/Utils/FMath.hpp"
#include "../../Src/Utils/FMemUtils.hpp"
#include "../../Src/Utils/FParameters.hpp"

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

28 29
#include "../../Src/GroupTree/Core/FGroupSeqAlgorithm.hpp"
#include "../../Src/GroupTree/Core/FGroupTaskAlgorithm.hpp"
30
#ifdef SCALFMM_USE_OMP4
31
#include "../../Src/GroupTree/Core/FGroupTaskDepAlgorithm.hpp"
32
#endif
33
#ifdef SCALFMM_USE_STARPU
34 35
#include "../../Src/GroupTree/Core/FGroupTaskStarpuAlgorithm.hpp"
#include "../../Src/GroupTree/StarPUUtils/FStarPUKernelCapacities.hpp"
36
#endif
37

38
#include "../../Src/GroupTree/Core/FP2PGroupParticleContainer.hpp"
39 40 41 42 43 44

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

#include <memory>


45
//#define RANDOM_PARTICLES
46 47 48 49 50

int main(int argc, char* argv[]){
    const FParameterNames LocalOptionBlocSize { {"-bs"}, "The size of the block of the blocked tree"};
    const FParameterNames LocalOptionNoValidate { {"-no-validation"}, "To avoid comparing with direct computation"};
    FHelpDescribeAndExit(argc, argv, "Test the blocked tree by counting the particles.",
51 52 53 54 55 56
                         FParameterDefinitions::OctreeHeight,
#ifdef RANDOM_PARTICLES
                         FParameterDefinitions::NbParticles,
#else
                         FParameterDefinitions::InputFile,
#endif
57 58 59
                         LocalOptionBlocSize, LocalOptionNoValidate);

    // Initialize the types
60
    typedef double FReal;
61
    static const int ORDER = 6;
62
    typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
63

64 65 66 67
    using GroupCellClass     = FChebCell<FReal, ORDER>;
    using GroupCellUpClass   = typename GroupCellClass::multipole_t;
    using GroupCellDownClass = typename GroupCellClass::local_expansion_t;
    using GroupCellSymbClass = FSymbolicData;
68

69
    typedef FP2PGroupParticleContainer<FReal>          GroupContainerClass;
70
    typedef FGroupTree< FReal, GroupCellSymbClass, GroupCellUpClass, GroupCellDownClass, GroupContainerClass, 1, 4, FReal>  GroupOctreeClass;
71
#ifdef SCALFMM_USE_STARPU
72
    typedef FStarPUAllCpuCapacities<FChebSymKernel<FReal,GroupCellClass,GroupContainerClass,MatrixKernelClass,ORDER>> GroupKernelClass;
73
    typedef FStarPUCpuWrapper<typename GroupOctreeClass::CellGroupClass, GroupCellClass, GroupKernelClass, typename GroupOctreeClass::ParticleGroupClass, GroupContainerClass> GroupCpuWrapper;
Berenger Bramas's avatar
Berenger Bramas committed
74
    typedef FGroupTaskStarPUAlgorithm<GroupOctreeClass, typename GroupOctreeClass::CellGroupClass, GroupKernelClass, typename GroupOctreeClass::ParticleGroupClass, GroupCpuWrapper, GroupContainerClass > GroupAlgorithm;
75
#elif defined(SCALFMM_USE_OMP4)
76
    typedef FChebSymKernel<FReal,GroupCellClass,GroupContainerClass,MatrixKernelClass,ORDER> GroupKernelClass;
77 78
    typedef FGroupTaskDepAlgorithm<GroupOctreeClass, typename GroupOctreeClass::CellGroupClass, GroupCellClass,
            GroupCellSymbClass, GroupCellUpClass, GroupCellDownClass, GroupKernelClass, typename GroupOctreeClass::ParticleGroupClass, GroupContainerClass > GroupAlgorithm;
79
#else
80
    typedef FChebSymKernel<FReal,GroupCellClass,GroupContainerClass,MatrixKernelClass,ORDER> GroupKernelClass;
81 82 83
    //typedef FGroupSeqAlgorithm<GroupOctreeClass, typename GroupOctreeClass::CellGroupClass, GroupCellClass, GroupKernelClass, typename GroupOctreeClass::ParticleGroupClass, GroupContainerClass > GroupAlgorithm;
    typedef FGroupTaskAlgorithm<GroupOctreeClass, typename GroupOctreeClass::CellGroupClass, GroupCellClass, GroupKernelClass, typename GroupOctreeClass::ParticleGroupClass, GroupContainerClass > GroupAlgorithm;
#endif
84 85 86 87 88
    // Get params
    const int NbLevels      = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 5);
    const int groupSize     = FParameters::getValue(argc,argv,LocalOptionBlocSize.options, 250);

    // Load the particles
89
#ifdef RANDOM_PARTICLES
90
    FRandomLoader<FReal> loader(FParameters::getValue(argc,argv,FParameterDefinitions::NbParticles.options, 2000), 1.0, FPoint<FReal>(0,0,0), 0);
91 92
#else
    const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma");
93
    FFmaGenericLoader<FReal> loader(filename);
94
#endif
95 96 97
    FAssertLF(loader.isOpen());
    FTic timer;

98
    FP2PParticleContainer<FReal> allParticles;
99
    for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
100
        FPoint<FReal> particlePosition;
101 102 103 104 105
        FReal physicalValue;
#ifdef RANDOM_PARTICLES
        physicalValue = 0.10;
        loader.fillParticle(&particlePosition);
#else
106
        loader.fillParticle(&particlePosition, &physicalValue);
107
#endif
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
        allParticles.push(particlePosition, physicalValue);
    }
    std::cout << "Particles loaded in " << timer.tacAndElapsed() << "s\n";

    // Put the data into the tree
    timer.tic();
    GroupOctreeClass groupedTree(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(), groupSize, &allParticles);
    groupedTree.printInfoBlocks();
    std::cout << "Tree created in " << timer.tacAndElapsed() << "s\n";

    // Run the algorithm
    const MatrixKernelClass MatrixKernel;
    GroupKernelClass groupkernel(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(), &MatrixKernel);
    GroupAlgorithm groupalgo(&groupedTree,&groupkernel);

    timer.tic();
124 125 126
#ifdef SCALFMM_USE_STARPU
    starpu_fxt_start_profiling();
#endif
127
    groupalgo.execute();
128 129 130
#ifdef SCALFMM_USE_STARPU
    starpu_fxt_stop_profiling();
#endif
131 132 133 134
    std::cout << "Kernel executed in in " << timer.tacAndElapsed() << "s\n";

    // Validate the result
    if(FParameters::existParameter(argc, argv, LocalOptionNoValidate.options) == false){
135
        FSize offsetParticles = 0;
136 137 138 139 140
        FReal*const allPhysicalValues = allParticles.getPhysicalValues();
        FReal*const allPosX = const_cast<FReal*>( allParticles.getPositions()[0]);
        FReal*const allPosY = const_cast<FReal*>( allParticles.getPositions()[1]);
        FReal*const allPosZ = const_cast<FReal*>( allParticles.getPositions()[2]);

141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
        groupedTree.forEachCellLeaf<FP2PGroupParticleContainer<FReal> >(
            [&](GroupCellSymbClass*  /* gsymb */,
                GroupCellUpClass*    /* gmul  */,
                GroupCellDownClass*  /* gloc  */,
                FP2PGroupParticleContainer<FReal> * leafTarget)
            {
                const FReal*const physicalValues = leafTarget->getPhysicalValues();
                const FReal*const posX = leafTarget->getPositions()[0];
                const FReal*const posY = leafTarget->getPositions()[1];
                const FReal*const posZ = leafTarget->getPositions()[2];
                const FSize nbPartsInLeafTarget = leafTarget->getNbParticles();

                for(FSize idxPart = 0 ; idxPart < nbPartsInLeafTarget ; ++idxPart){
                    allPhysicalValues[offsetParticles + idxPart] = physicalValues[idxPart];
                    allPosX[offsetParticles + idxPart] = posX[idxPart];
                    allPosY[offsetParticles + idxPart] = posY[idxPart];
                    allPosZ[offsetParticles + idxPart] = posZ[idxPart];
                }

                offsetParticles += nbPartsInLeafTarget;
            });
162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180

        FAssertLF(offsetParticles == loader.getNumberOfParticles());

        FReal*const allDirectPotentials = allParticles.getPotentials();
        FReal*const allDirectforcesX = allParticles.getForcesX();
        FReal*const allDirectforcesY = allParticles.getForcesY();
        FReal*const allDirectforcesZ = allParticles.getForcesZ();

        for(int idxTgt = 0 ; idxTgt < offsetParticles ; ++idxTgt){
            for(int idxMutual = idxTgt + 1 ; idxMutual < offsetParticles ; ++idxMutual){
                FP2PR::MutualParticles(
                    allPosX[idxTgt],allPosY[idxTgt],allPosZ[idxTgt], allPhysicalValues[idxTgt],
                    &allDirectforcesX[idxTgt], &allDirectforcesY[idxTgt], &allDirectforcesZ[idxTgt], &allDirectPotentials[idxTgt],
                    allPosX[idxMutual],allPosY[idxMutual],allPosZ[idxMutual], allPhysicalValues[idxMutual],
                    &allDirectforcesX[idxMutual], &allDirectforcesY[idxMutual], &allDirectforcesZ[idxMutual], &allDirectPotentials[idxMutual]
                );
            }
        }

181 182
        FMath::FAccurater<FReal> potentialDiff;
        FMath::FAccurater<FReal> fx, fy, fz;
183
        offsetParticles = 0;
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203
        groupedTree.forEachCellLeaf<FP2PGroupParticleContainer<FReal> >(
            [&](GroupCellSymbClass*  /* gsymb */,
                GroupCellUpClass*    /* gmul  */,
                GroupCellDownClass*  /* gloc  */,
                FP2PGroupParticleContainer<FReal> * leafTarget){
                const FReal*const potentials = leafTarget->getPotentials();
                const FReal*const forcesX = leafTarget->getForcesX();
                const FReal*const forcesY = leafTarget->getForcesY();
                const FReal*const forcesZ = leafTarget->getForcesZ();
                const FSize nbPartsInLeafTarget = leafTarget->getNbParticles();

                for(int idxTgt = 0 ; idxTgt < nbPartsInLeafTarget ; ++idxTgt){
                    potentialDiff.add(allDirectPotentials[idxTgt + offsetParticles], potentials[idxTgt]);
                    fx.add(allDirectforcesX[idxTgt + offsetParticles], forcesX[idxTgt]);
                    fy.add(allDirectforcesY[idxTgt + offsetParticles], forcesY[idxTgt]);
                    fz.add(allDirectforcesZ[idxTgt + offsetParticles], forcesZ[idxTgt]);
                }

                offsetParticles += nbPartsInLeafTarget;
            });
204 205 206 207 208 209 210 211 212

        std::cout << "Error : Potential " << potentialDiff << "\n";
        std::cout << "Error : fx " << fx << "\n";
        std::cout << "Error : fy " << fy << "\n";
        std::cout << "Error : fz " << fz << "\n";
    }

    return 0;
}