FFmmAlgorithmTsm.hpp 14.1 KB
Newer Older
1
// ===================================================================================
2
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
3 4 5 6 7 8 9 10 11 12 13 14
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.  
// 
// 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
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info". 
// "http://www.gnu.org/licenses".
15
// ===================================================================================
16 17
#ifndef FFMMALGORITHMTSM_HPP
#define FFMMALGORITHMTSM_HPP
18

19

BRAMAS Berenger's avatar
BRAMAS Berenger committed
20
#include "../Utils/FAssert.hpp"
BRAMAS Berenger's avatar
BRAMAS Berenger committed
21
#include "../Utils/FLog.hpp"
22

23 24 25
#include "../Utils/FTic.hpp"

#include "../Containers/FOctree.hpp"
26
#include "FCoreCommon.hpp"
27 28 29 30 31 32 33 34 35 36 37

/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FFmmAlgorithmTsm
* @brief
* Please read the license
*
* This class is a basic FMM algorithm
* It just iterates on a tree and call the kernels with good arguments.
*
* Of course this class does not deallocate pointer given in arguements.
berenger-bramas's avatar
berenger-bramas committed
38 39
*
* The differences with FmmAlgorithm is that it used target source model.
40
*/
41
template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
BRAMAS Berenger's avatar
BRAMAS Berenger committed
42
class FFmmAlgorithmTsm : public FAbstractAlgorithm{
43

berenger-bramas's avatar
berenger-bramas committed
44 45
    OctreeClass* const tree;                                                     //< The octree to work on
    KernelClass* const kernels;    //< The kernels
46 47

    const int OctreeHeight;
48

49 50
    const int leafLevelSeperationCriteria;

BRAMAS Berenger's avatar
BRAMAS Berenger committed
51 52
    FLOG(FTic counterTime);                                               //< In case of debug: to count the elapsed time
    FLOG(FTic computationCounter);                                        //< In case of debug: to  count computation time
53 54 55 56 57 58 59

public:
    /** The constructor need the octree and the kernels used for computation
      * @param inTree the octree to work on
      * @param inKernels the kernels to call
      * An assert is launched if one of the arguments is null
      */
60 61
    FFmmAlgorithmTsm(OctreeClass* const inTree, KernelClass* const inKernels, const int inLeafLevelSeperationCriteria = 1)
        : tree(inTree) , kernels(inKernels) , OctreeHeight(tree->getHeight()), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria){
62

BRAMAS Berenger's avatar
BRAMAS Berenger committed
63 64
        FAssertLF(tree, "tree cannot be null");
        FAssertLF(kernels, "kernels cannot be null");
65

66 67
        FAbstractAlgorithm::setNbLevelsInTree(tree->getHeight());

BRAMAS Berenger's avatar
BRAMAS Berenger committed
68
        FLOG(FLog::Controller << "FFmmAlgorithmTsm\n");
69 70 71 72 73 74
    }

    /** Default destructor */
    virtual ~FFmmAlgorithmTsm(){
    }

75
protected:
76 77 78 79
    /**
      * To execute the fmm algorithm
      * Call this function to run the complete algorithm
      */
80
    void executeCore(const unsigned operationsToProceed) override {
81

82
        if(operationsToProceed & FFmmP2M) bottomPass();
83

84
        if(operationsToProceed & FFmmM2M) upwardPass();
85

86
        if(operationsToProceed & FFmmM2L) transferPass();
87

88
        if(operationsToProceed & FFmmL2L) downardPass();
89

90
        if((operationsToProceed & FFmmP2P) || (operationsToProceed & FFmmL2P)) directPass((operationsToProceed & FFmmP2P),(operationsToProceed & FFmmL2P));
91 92 93 94
    }

    /** P2M */
    void bottomPass(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
95
        FLOG( FLog::Controller.write("\tStart Bottom Pass\n").write(FLog::Flush) );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
96 97
        FLOG( counterTime.tic() );
        FLOG( double totalComputation = 0 );
98

berenger-bramas's avatar
berenger-bramas committed
99
        typename OctreeClass::Iterator octreeIterator(tree);
100 101 102 103 104 105

        // Iterate on leafs
        octreeIterator.gotoBottomLeft();
        do{
            // We need the current cell that represent the leaf
            // and the list of particles
BRAMAS Berenger's avatar
BRAMAS Berenger committed
106
            FLOG(computationCounter.tic());
berenger-bramas's avatar
berenger-bramas committed
107
            ContainerClass* const sources = octreeIterator.getCurrentListSrc();
108
            if(sources->getNbParticles()){
109 110 111
                octreeIterator.getCurrentCell()->setSrcChildTrue();
                kernels->P2M( octreeIterator.getCurrentCell() , sources);
            }
112
            if(octreeIterator.getCurrentListTargets()->getNbParticles()){
113 114
                octreeIterator.getCurrentCell()->setTargetsChildTrue();
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
115 116
            FLOG(computationCounter.tac());
            FLOG(totalComputation += computationCounter.elapsed());
117 118
        } while(octreeIterator.moveRight());

BRAMAS Berenger's avatar
BRAMAS Berenger committed
119
        FLOG( counterTime.tac() );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
120 121
        FLOG( FLog::Controller << "\tFinished (@Bottom Pass (P2M) = "  << counterTime.elapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << totalComputation << " s\n" );
122

123 124 125 126
    }

    /** M2M */
    void upwardPass(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
127
        FLOG( FLog::Controller.write("\tStart Upward Pass\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
128 129
        FLOG( counterTime.tic() );
        FLOG( double totalComputation = 0 );
130 131

        // Start from leal level - 1
berenger-bramas's avatar
berenger-bramas committed
132
        typename OctreeClass::Iterator octreeIterator(tree);
133 134 135
        octreeIterator.gotoBottomLeft();
        octreeIterator.moveUp();

136
        for(int idxLevel = OctreeHeight - 2 ; idxLevel > FAbstractAlgorithm::lowerWorkingLevel-1 ; --idxLevel){
137 138 139
            octreeIterator.moveUp();
        }

berenger-bramas's avatar
berenger-bramas committed
140
        typename OctreeClass::Iterator avoidGotoLeftIterator(octreeIterator);
141 142

        // for each levels
143
        for(int idxLevel = FMath::Min(OctreeHeight - 2, FAbstractAlgorithm::lowerWorkingLevel - 1) ; idxLevel >= FAbstractAlgorithm::upperWorkingLevel ; --idxLevel ){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
144
            FLOG(FTic counterTimeLevel);
145 146 147 148
            // for each cells
            do{
                // We need the current cell and the child
                // child is an array (of 8 child) that may be null
BRAMAS Berenger's avatar
BRAMAS Berenger committed
149
                FLOG(computationCounter.tic());
150 151 152 153 154

                CellClass* potentialChild[8];
                CellClass** const realChild = octreeIterator.getCurrentChild();
                CellClass* const currentCell = octreeIterator.getCurrentCell();
                for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
155
                    potentialChild[idxChild] = nullptr;
156 157 158 159 160 161 162 163 164 165 166 167
                    if(realChild[idxChild]){
                        if(realChild[idxChild]->hasSrcChild()){
                            currentCell->setSrcChildTrue();
                            potentialChild[idxChild] = realChild[idxChild];
                        }
                        if(realChild[idxChild]->hasTargetsChild()){
                            currentCell->setTargetsChildTrue();
                        }
                    }
                }
                kernels->M2M( currentCell , potentialChild, idxLevel);

BRAMAS Berenger's avatar
BRAMAS Berenger committed
168 169
                FLOG(computationCounter.tac());
                FLOG(totalComputation += computationCounter.elapsed());
170 171 172 173
            } while(octreeIterator.moveRight());

            avoidGotoLeftIterator.moveUp();
            octreeIterator = avoidGotoLeftIterator;// equal octreeIterator.moveUp(); octreeIterator.gotoLeft();
BRAMAS Berenger's avatar
BRAMAS Berenger committed
174
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << " = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
175 176
        }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
177
        FLOG( counterTime.tac() );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
178 179
        FLOG( FLog::Controller << "\tFinished (@Upward Pass (M2M) = "  << counterTime.elapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << totalComputation << " s\n" );
180

181 182
    }

183 184
    /** M2L */
    void transferPass(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
185
        FLOG( FLog::Controller.write("\tStart Downward Pass (M2L)\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
186 187
        FLOG( counterTime.tic() );
        FLOG( double totalComputation = 0 );
188

189 190 191
        typename OctreeClass::Iterator octreeIterator(tree);
        octreeIterator.moveDown();

192
        for(int idxLevel = 2 ; idxLevel < FAbstractAlgorithm::upperWorkingLevel ; ++idxLevel){
193 194 195
            octreeIterator.moveDown();
        }

196 197
        typename OctreeClass::Iterator avoidGotoLeftIterator(octreeIterator);

198
        const CellClass* neighbors[343];
199 200

        // for each levels
201
        for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
202
            FLOG(FTic counterTimeLevel);
203
            const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeperationCriteria);
204 205
            // for each cells
            do{
BRAMAS Berenger's avatar
BRAMAS Berenger committed
206
                FLOG(computationCounter.tic());
207
                CellClass* const currentCell = octreeIterator.getCurrentCell();
BRAMAS Berenger's avatar
BRAMAS Berenger committed
208

209
                if(currentCell->hasTargetsChild()){
210
                    const int counter = tree->getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(),idxLevel, separationCriteria);
211 212 213 214 215 216 217
                    if( counter ){
                        int counterWithSrc = 0;
                        for(int idxRealNeighbors = 0 ; idxRealNeighbors < 343 ; ++idxRealNeighbors ){
                            if(neighbors[idxRealNeighbors] && neighbors[idxRealNeighbors]->hasSrcChild()){
                                ++counterWithSrc;
                            }
                            else{
218
                                neighbors[idxRealNeighbors] = nullptr;
219 220
                            }
                        }
221 222
                        if(counterWithSrc){
                            kernels->M2L( currentCell , neighbors, counterWithSrc, idxLevel);
223 224
                        }
                    }
225
                }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
226 227
                FLOG(computationCounter.tac());
                FLOG(totalComputation += computationCounter.elapsed());
228
            } while(octreeIterator.moveRight());
229

BRAMAS Berenger's avatar
BRAMAS Berenger committed
230
            FLOG(computationCounter.tic());
231
            kernels->finishedLevelM2L(idxLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
232
            FLOG(computationCounter.tac());
233

234 235
            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
236
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << " = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
237
        }
238

BRAMAS Berenger's avatar
BRAMAS Berenger committed
239
        FLOG( counterTime.tac() );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
240 241
        FLOG( FLog::Controller << "\tFinished (@Downward Pass (M2L) = "  << counterTime.elapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << totalComputation << " s\n" );
242
    }
243

244 245
    /** L2L */
    void downardPass(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
246
        FLOG( FLog::Controller.write("\tStart Downward Pass (L2L)\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
247 248
        FLOG( counterTime.tic() );
        FLOG( double totalComputation = 0 );
249 250 251 252

        typename OctreeClass::Iterator octreeIterator(tree);
        octreeIterator.moveDown();

253
        for(int idxLevel = 2 ; idxLevel < FAbstractAlgorithm::upperWorkingLevel ; ++idxLevel){
254 255 256
            octreeIterator.moveDown();
        }

257 258
        typename OctreeClass::Iterator avoidGotoLeftIterator(octreeIterator);

259
        const int heightMinusOne = FAbstractAlgorithm::lowerWorkingLevel - 1;
260
        // for each levels exepted leaf level
261
        for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < heightMinusOne ; ++idxLevel ){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
262
            FLOG(FTic counterTimeLevel);
263 264
            // for each cells
            do{
BRAMAS Berenger's avatar
BRAMAS Berenger committed
265 266 267 268 269 270 271 272 273 274
                if( octreeIterator.getCurrentCell()->hasTargetsChild() ){
                    FLOG(computationCounter.tic());
                    CellClass* potentialChild[8];
                    CellClass** const realChild = octreeIterator.getCurrentChild();
                    CellClass* const currentCell = octreeIterator.getCurrentCell();
                    for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
                        if(realChild[idxChild] && realChild[idxChild]->hasTargetsChild()){
                            potentialChild[idxChild] = realChild[idxChild];
                        }
                        else{
275
                            potentialChild[idxChild] = nullptr;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
276
                        }
277
                    }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
278 279 280
                    kernels->L2L( currentCell , potentialChild, idxLevel);
                    FLOG(computationCounter.tac());
                    FLOG(totalComputation += computationCounter.elapsed());
281 282
                }
            } while(octreeIterator.moveRight());
283

284 285
            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
286
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << " = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
287 288
        }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
289
        FLOG( counterTime.tac() );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
290 291
        FLOG( FLog::Controller << "\tFinished (@Downward Pass (L2L) = "  << counterTime.elapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << totalComputation << " s\n" );
292 293
    }

294 295


296
    /** P2P */
297
    void directPass(const bool p2pEnabled, const bool l2pEnabled){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
298
        FLOG( FLog::Controller.write("\tStart Direct Pass\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
299 300
        FLOG( counterTime.tic() );
        FLOG( double totalComputation = 0 );
301 302 303

        const int heightMinusOne = OctreeHeight - 1;

berenger-bramas's avatar
berenger-bramas committed
304
        typename OctreeClass::Iterator octreeIterator(tree);
305 306
        octreeIterator.gotoBottomLeft();
        // There is a maximum of 26 neighbors
berenger-bramas's avatar
berenger-bramas committed
307
        ContainerClass* neighbors[27];
308 309
        // for each leafs
        do{
BRAMAS Berenger's avatar
BRAMAS Berenger committed
310 311
            if( octreeIterator.getCurrentCell()->hasTargetsChild() ){
                FLOG(computationCounter.tic());
312 313 314 315 316 317 318 319
                if(l2pEnabled){
                    kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
                }
                if(p2pEnabled){
                    // need the current particles and neighbors particles
                    const int counter = tree->getLeafsNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), heightMinusOne);
                    neighbors[13] = octreeIterator.getCurrentListSrc();
                    kernels->P2PRemote( octreeIterator.getCurrentGlobalCoordinate(), octreeIterator.getCurrentListTargets(),
BRAMAS Berenger's avatar
BRAMAS Berenger committed
320
                              octreeIterator.getCurrentListSrc() , neighbors, counter + 1);
321
                }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
322 323 324
                FLOG(computationCounter.tac());
                FLOG(totalComputation += computationCounter.elapsed());
            }
325 326
        } while(octreeIterator.moveRight());

BRAMAS Berenger's avatar
BRAMAS Berenger committed
327
        FLOG( counterTime.tac() );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
328 329
        FLOG( FLog::Controller << "\tFinished (@Direct Pass (L2P + P2P) = "  << counterTime.elapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation L2P + P2P : " << totalComputation << " s\n" );
330

331 332 333 334 335 336 337
    }

};


#endif //FFMMALGORITHMTSM_HPP

338