testSphericalGalaxyCsv.cpp 7.97 KB
Newer Older
1
// See LICENCE file at project root
2

3 4 5 6 7
#include <iostream>

#include <cstdio>
#include <cstdlib>

8 9
#include "Utils/FTic.hpp"
#include "Utils/FParameters.hpp"
10

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

14 15
#include "Core/FFmmAlgorithm.hpp"
#include "Core/FFmmAlgorithmThread.hpp"
16

17 18
#include "Kernels/Spherical/FSphericalKernel.hpp"
#include "Kernels/Spherical/FSphericalCell.hpp"
19

20 21 22
#include "Files/FTreeCsvSaver.hpp"
#include "Files/FFmaGenericLoader.hpp"
#include "Arranger/FOctreeArranger.hpp"
23

24 25
#include "Components/FSimpleLeaf.hpp"
#include "Components/FParticleType.hpp"
26

27
#include "Kernels/P2P/FP2PParticleContainer.hpp"
28

29 30
#include "Utils/FParameterNames.hpp"
#include "Arranger/FAbstractMover.hpp"
31 32


33 34 35 36
template <class FReal>
class VelocityContainer : public FP2PParticleContainer<FReal> {
    typedef FP2PParticleContainer<FReal> Parent;
    FVector<FPoint<FReal>> velocities;
37 38 39

public:
    template<typename... Args>
40
    void push(const FPoint<FReal>& inParticlePosition, const FPoint<FReal>& velocity, Args... args){
41 42 43 44
        Parent::push(inParticlePosition, args... );
        velocities.push(velocity);
    }

45
    const FVector<FPoint<FReal>>& getVelocities() const{
46 47 48
        return velocities;
    }

49
    FVector<FPoint<FReal>>& getVelocities() {
50 51 52
        return velocities;
    }

53
    void fillToCsv(const FSize partIdx, FReal values[4]) const {
54 55 56 57 58
        values[0] = Parent::getPositions()[0][partIdx];
        values[1] = Parent::getPositions()[1][partIdx];
        values[2] = Parent::getPositions()[2][partIdx];
        values[3] = Parent::getPotentials()[partIdx];
    }
59 60
};

61

62 63
template <class FReal>
class GalaxyLoader : public FFmaGenericLoader<FReal> {
64
public:
65
    GalaxyLoader(const char* const filename) : FFmaGenericLoader<FReal>(filename) {
66 67
    }

68
    void fillParticle(FPoint<FReal>* position, FReal* physivalValue, FPoint<FReal>* velocity){
69
        FReal x,y,z,data, vx, vy, vz;
70
        *(FFmaGenericLoader<FReal>::file) >> x >> y >> z >> data >> vx >> vy >> vz;
71 72 73
        position->setPosition(x,y,z);
        *physivalValue = (data);
        velocity->setPosition(vx,vy,vz);
74 75 76
    }
};

77 78
template<class FReal, class OctreeClass>
class GalaxyMover : public FAbstractMover<FReal, OctreeClass, VelocityContainer<FReal>>{
79
private:
80
    VelocityContainer<FReal> toStoreRemovedParts;
81 82

public:
83 84 85 86
    GalaxyMover() {
    }

    virtual ~GalaxyMover(){
87 88
    }

89
    /** To get the position of the particle at idx idxPart in leaf lf */
90
    void getParticlePosition(VelocityContainer<FReal>* lf, const FSize idxPart, FPoint<FReal>* particlePos){
91
        (*particlePos) = FPoint<FReal>(lf->getPositions()[0][idxPart],lf->getPositions()[1][idxPart],lf->getPositions()[2][idxPart]);
92 93 94
    }

    /** Remove a particle but keep it to reinsert it later*/
95
    void removeFromLeafAndKeep(VelocityContainer<FReal>* lf, const FPoint<FReal>& particlePos, const FSize idxPart,FParticleType /*type*/){
96 97
        std::array<typename VelocityContainer<FReal>::AttributesClass, VelocityContainer<FReal>::NbAttributes> particleValues;
        for(int idxAttr = 0 ; idxAttr < VelocityContainer<FReal>::NbAttributes ; ++idxAttr){
98 99 100 101 102 103 104 105 106 107 108
            particleValues[idxAttr] = lf->getAttribute(idxAttr)[idxPart];
        }

        toStoreRemovedParts.push(particlePos,lf->getVelocities()[idxPart],particleValues);

        lf->getVelocities().removeOne(idxPart);
        lf->removeParticles(&idxPart,1);
    }

    /** Reinsert the previously saved particles */
    void insertAllParticles(OctreeClass* tree){
109
        std::array<typename VelocityContainer<FReal>::AttributesClass, VelocityContainer<FReal>::NbAttributes> particleValues;
110

111
        for(FSize idxToInsert = 0; idxToInsert<toStoreRemovedParts.getNbParticles() ; ++idxToInsert){
112
            for(int idxAttr = 0 ; idxAttr < VelocityContainer<FReal>::NbAttributes ; ++idxAttr){
113 114
                particleValues[idxAttr] = toStoreRemovedParts.getAttribute(idxAttr)[idxToInsert];
            }
115
            const FPoint<FReal> particlePos(toStoreRemovedParts.getPositions()[0][idxToInsert],
116 117 118 119 120 121 122 123
                                     toStoreRemovedParts.getPositions()[1][idxToInsert],
                                     toStoreRemovedParts.getPositions()[2][idxToInsert]);

            tree->insert(particlePos, toStoreRemovedParts.getVelocities()[idxToInsert], particleValues);
        }

        toStoreRemovedParts.clear();
        toStoreRemovedParts.getVelocities().clear();
124 125 126
    }
};

127

128 129
// Simply create particles and try the kernels
int main(int argc, char ** argv){
130 131 132 133 134 135
    FHelpDescribeAndExit(argc, argv,
                         "Run a Spherical Harmonic (Old Implementation) FMM kernel with several time step.",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::SHDevelopment,
                         FParameterDefinitions::DeltaT, FParameterDefinitions::OutputFile);

136 137
    typedef double FReal;
    typedef FSphericalCell<FReal>          CellClass;
138
    typedef VelocityContainer<FReal>  ContainerClass;
139

140 141 142
    typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
    typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
    typedef FSphericalKernel<FReal, CellClass, ContainerClass >   KernelClass;
143

144
    typedef FFmmAlgorithmThread<OctreeClass,  CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
145

146 147
    typedef GalaxyMover<FReal, OctreeClass> MoverClass;
    typedef FOctreeArranger<FReal,OctreeClass, ContainerClass, MoverClass> ArrangerClass;
148
    ///////////////////////What we do/////////////////////////////
149
    std::cout << ">> This executable has to be used to test Spherical algorithm.\n";
150 151
    //////////////////////////////////////////////////////////////

152 153 154 155
    const int NbLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 6);
    const int SizeSubLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 3);
    const FReal DT          = FParameters::getValue(argc,argv,FParameterDefinitions::DeltaT.options, FReal(0.1));
    const int DevP          = FParameters::getValue(argc,argv,FParameterDefinitions::SHDevelopment.options, 5);
156

157
    FSphericalCell<FReal>::Init(DevP);
158

159
    GalaxyLoader<FReal> loader(FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/galaxy.fma.tmp"));
160 161 162 163 164 165 166 167 168 169 170

    // -----------------------------------------------------

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

    // -----------------------------------------------------

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

    {
171
        FPoint<FReal> position, velocity;
172
        FReal physicalValue;
173

174
        for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
175 176
            loader.fillParticle(&position, &physicalValue, &velocity);
            tree.insert(position, velocity, physicalValue);
177 178 179 180 181
        }
    }

    // -----------------------------------------------------

182
    KernelClass kernels( DevP, NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
183
    FmmClass algo( &tree, &kernels);
184
    ArrangerClass arranger(&tree);
185
    FTreeCsvSaver<FReal, OctreeClass, ContainerClass> saver(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "/tmp/test%d.csv"));
186 187 188 189

    for(int idx = 0; idx < 100 ; ++idx){
        algo.execute();
        { // update velocity and position
190
            OctreeClass::Iterator octreeIterator(&tree);
191 192
            octreeIterator.gotoBottomLeft();
            do{
193 194
                kernels.computeVelocity(octreeIterator.getCurrentListTargets(), DT);
                kernels.updatePosition(octreeIterator.getCurrentListTargets(), DT);
195 196 197
            } while(octreeIterator.moveRight());
        }
        // update tree and vtk
198
        arranger.rearrange();
199 200 201 202 203 204 205
        saver.exportTree(&tree);
    }

    // -----------------------------------------------------

    return 0;
}