testBlockedUniform.cpp 10 KB
Newer Older
BRAMAS Berenger's avatar
BRAMAS Berenger committed
1 2 3 4

// ==== CMAKE =====
// @FUSE_BLAS
// @FUSE_FFT
5
// @FUSE_STARPU
BRAMAS Berenger's avatar
BRAMAS Berenger committed
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
// ================
// Keep in private GIT


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

#include "../../Src/GroupTree/Core/FGroupTree.hpp"

#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Containers/FVector.hpp"

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

#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "../../Src/Kernels/Uniform/FUnifKernel.hpp"

#include "../../Src/GroupTree/Uniform/FUnifCellPOD.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"

#include "../../Src/GroupTree/Core/FGroupSeqAlgorithm.hpp"
#include "../../Src/GroupTree/Core/FGroupTaskAlgorithm.hpp"
33
#ifdef SCALFMM_USE_OMP4
BRAMAS Berenger's avatar
BRAMAS Berenger committed
34 35
#include "../../Src/GroupTree/Core/FGroupTaskDepAlgorithm.hpp"
#endif
36
#ifdef SCALFMM_USE_STARPU
BRAMAS Berenger's avatar
BRAMAS Berenger committed
37 38 39 40 41 42 43 44 45 46
#include "../../Src/GroupTree/Core/FGroupTaskStarpuAlgorithm.hpp"
#include "../../Src/GroupTree/StarPUUtils/FStarPUKernelCapacities.hpp"
#endif
#include "../../Src/GroupTree/Core/FP2PGroupParticleContainer.hpp"

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

#include <memory>


47
#define RANDOM_PARTICLES
BRAMAS Berenger's avatar
BRAMAS Berenger committed
48

49 50 51 52
#ifdef STARPU_SIMGRID_MLR_MODELS
extern "C" {
#endif

BRAMAS Berenger's avatar
BRAMAS Berenger committed
53 54 55
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"};
56
    FHelpDescribeAndExit(argc, argv, "Perform Lagrange Kernel based simulation with StarPU",
BRAMAS Berenger's avatar
BRAMAS Berenger committed
57 58 59 60 61 62 63 64 65 66
                         FParameterDefinitions::OctreeHeight,
#ifdef RANDOM_PARTICLES
                         FParameterDefinitions::NbParticles,
#else
                         FParameterDefinitions::InputFile,
#endif
                         FParameterDefinitions::NbThreads,
                         LocalOptionBlocSize, LocalOptionNoValidate);

    // Initialize the types
67
    typedef double FReal;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
68
    static const int ORDER = 5;
69
    typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
70 71

    typedef FUnifCellPODCore         GroupCellSymbClass;
72 73 74
    typedef FUnifCellPODPole<FReal,ORDER>  GroupCellUpClass;
    typedef FUnifCellPODLocal<FReal,ORDER> GroupCellDownClass;
    typedef FUnifCellPOD<FReal,ORDER>      GroupCellClass;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
75

76 77
    typedef FP2PGroupParticleContainer<FReal>          GroupContainerClass;
    typedef FGroupTree< FReal, GroupCellClass, GroupCellSymbClass, GroupCellUpClass, GroupCellDownClass, GroupContainerClass, 1, 4, FReal>  GroupOctreeClass;
78
#ifdef SCALFMM_USE_STARPU
79
    typedef FStarPUAllCpuCapacities<FUnifKernel<FReal,GroupCellClass,GroupContainerClass,MatrixKernelClass,ORDER>> GroupKernelClass;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
80
    typedef FStarPUCpuWrapper<typename GroupOctreeClass::CellGroupClass, GroupCellClass, GroupKernelClass, typename GroupOctreeClass::ParticleGroupClass, GroupContainerClass> GroupCpuWrapper;
Berenger Bramas's avatar
Berenger Bramas committed
81
    typedef FGroupTaskStarPUAlgorithm<GroupOctreeClass, typename GroupOctreeClass::CellGroupClass, GroupKernelClass, typename GroupOctreeClass::ParticleGroupClass, GroupCpuWrapper, GroupContainerClass > GroupAlgorithm;
82
#elif defined(SCALFMM_USE_OMP4)
83
    typedef FUnifKernel<FReal,GroupCellClass,GroupContainerClass,MatrixKernelClass,ORDER> GroupKernelClass;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
84 85
    // Set the number of threads
    omp_set_num_threads(FParameters::getValue(argc,argv,FParameterDefinitions::NbThreads.options, omp_get_max_threads()));
86 87
    typedef FGroupTaskDepAlgorithm<GroupOctreeClass, typename GroupOctreeClass::CellGroupClass, GroupCellClass,
            GroupCellSymbClass, GroupCellUpClass, GroupCellDownClass, GroupKernelClass, typename GroupOctreeClass::ParticleGroupClass, GroupContainerClass > GroupAlgorithm;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
88
#else
89
    typedef FUnifKernel<FReal,GroupCellClass,GroupContainerClass,MatrixKernelClass,ORDER> GroupKernelClass;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
90 91 92 93 94 95 96 97 98
    //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
    // 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
#ifdef RANDOM_PARTICLES
99
    FRandomLoader<FReal> loader(FParameters::getValue(argc,argv,FParameterDefinitions::NbParticles.options, 2000), 1.0, FPoint<FReal>(0,0,0), 0);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
100 101
#else
    const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma");
102
    FFmaGenericLoader<FReal> loader(filename);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
103 104 105 106
#endif
    FAssertLF(loader.isOpen());
    FTic timer;

107
    FP2PParticleContainer<FReal> allParticles;
108
    for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
109
        FPoint<FReal> particlePosition;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
        FReal physicalValue;
#ifdef RANDOM_PARTICLES
        physicalValue = 0.10;
        loader.fillParticle(&particlePosition);
#else
        loader.fillParticle(&particlePosition, &physicalValue);
#endif
        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);

132
    // Extended for Native vs SimGrid makespans comparison
BRAMAS Berenger's avatar
BRAMAS Berenger committed
133
    timer.tic();
134
    double start_time = starpu_timing_now();
BRAMAS Berenger's avatar
BRAMAS Berenger committed
135
    groupalgo.execute();
136
    double end_time = starpu_timing_now();
BRAMAS Berenger's avatar
BRAMAS Berenger committed
137
    std::cout << "Kernel executed in in " << timer.tacAndElapsed() << "s\n";
138
    std::cout << (end_time - start_time)/1000 << "\n";
BRAMAS Berenger's avatar
BRAMAS Berenger committed
139 140 141

    // Validate the result
    if(FParameters::existParameter(argc, argv, LocalOptionNoValidate.options) == false){
142
        FSize offsetParticles = 0;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
143 144 145 146 147
        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]);

148
        groupedTree.forEachCellLeaf<FP2PGroupParticleContainer<FReal> >([&](GroupCellClass cellTarget, FP2PGroupParticleContainer<FReal> * leafTarget){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
149 150 151 152
            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];
153
            const FSize nbPartsInLeafTarget = leafTarget->getNbParticles();
BRAMAS Berenger's avatar
BRAMAS Berenger committed
154

155
            for(FSize idxPart = 0 ; idxPart < nbPartsInLeafTarget ; ++idxPart){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
156 157 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
                allPhysicalValues[offsetParticles + idxPart] = physicalValues[idxPart];
                allPosX[offsetParticles + idxPart] = posX[idxPart];
                allPosY[offsetParticles + idxPart] = posY[idxPart];
                allPosZ[offsetParticles + idxPart] = posZ[idxPart];
            }

            offsetParticles += nbPartsInLeafTarget;
        });

        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]
                );
            }
        }

183 184
        FMath::FAccurater<FReal> potentialDiff;
        FMath::FAccurater<FReal> fx, fy, fz;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
185
        offsetParticles = 0;
186
        groupedTree.forEachCellLeaf<FP2PGroupParticleContainer<FReal> >([&](GroupCellClass cellTarget, FP2PGroupParticleContainer<FReal> * leafTarget){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
187 188 189 190
            const FReal*const potentials = leafTarget->getPotentials();
            const FReal*const forcesX = leafTarget->getForcesX();
            const FReal*const forcesY = leafTarget->getForcesY();
            const FReal*const forcesZ = leafTarget->getForcesZ();
191
            const FSize nbPartsInLeafTarget = leafTarget->getNbParticles();
BRAMAS Berenger's avatar
BRAMAS Berenger committed
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210

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

        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;
}
211

212 213 214
#ifdef STARPU_SIMGRID_MLR_MODELS
}
#endif