FGroupSeqAlgorithm.hpp 21.8 KB
Newer Older
1 2 3

// Keep in private GIT
// @SCALFMM_PRIVATE
4 5 6
#ifndef FGROUPSEQALGORITHM_HPP
#define FGROUPSEQALGORITHM_HPP

7 8 9 10 11 12
#include "../../Utils/FGlobal.hpp"
#include "../../Core/FCoreCommon.hpp"
#include "../../Utils/FQuickSort.hpp"
#include "../../Containers/FTreeCoordinate.hpp"
#include "../../Utils/FLog.hpp"
#include "../../Utils/FTic.hpp"
13

14 15
#include "FOutOfBlockInteraction.hpp"

BRAMAS Berenger's avatar
BRAMAS Berenger committed
16
#include <vector>
17 18
#include <vector>

19
template <class OctreeClass, class CellContainerClass, class CellClass, class KernelClass, class ParticleGroupClass, class ParticleContainerClass>
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
class FGroupSeqAlgorithm {
protected:
    const int MaxThreads;         //< The number of threads
    OctreeClass*const tree;       //< The Tree
    KernelClass*const kernels;    //< The kernels

public:
    FGroupSeqAlgorithm(OctreeClass*const inTree, KernelClass* inKernels) : MaxThreads(1), tree(inTree), kernels(inKernels){
        FAssertLF(tree, "tree cannot be null");
        FAssertLF(kernels, "kernels cannot be null");

        FLOG(FLog::Controller << "FGroupSeqAlgorithm (Max Thread " << MaxThreads << ")\n");
    }

    ~FGroupSeqAlgorithm(){
    }

    void execute(const unsigned operationsToProceed = FFmmNearAndFarFields){
38 39
        FLOG( FLog::Controller << "\tStart FGroupSeqAlgorithm\n" );

40 41 42 43 44 45 46 47 48 49 50 51 52
        if(operationsToProceed & FFmmP2M) bottomPass();

        if(operationsToProceed & FFmmM2M) upwardPass();

        if(operationsToProceed & FFmmM2L) transferPass();

        if(operationsToProceed & FFmmL2L) downardPass();

        if( (operationsToProceed & FFmmP2P) || (operationsToProceed & FFmmL2P) ) directPass();
    }

protected:
    void bottomPass(){
53
        FLOG( FTic timer; );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
54 55
        typename OctreeClass::ParticleGroupIterator iterParticles = tree->leavesBegin();
        const typename OctreeClass::ParticleGroupIterator endParticles = tree->leavesEnd();
56

BRAMAS Berenger's avatar
BRAMAS Berenger committed
57 58
        typename OctreeClass::CellGroupIterator iterCells = tree->cellsBegin(tree->getHeight()-1);
        const typename OctreeClass::CellGroupIterator endCells = tree->cellsEnd(tree->getHeight()-1);
59 60 61 62 63 64 65 66 67

        while(iterParticles != endParticles && iterCells != endCells){
            { // Can be a task(in:iterParticles, out:iterCells)
                const MortonIndex blockStartIdx = (*iterCells)->getStartingIndex();
                const MortonIndex blockEndIdx = (*iterCells)->getEndingIndex();

                for(MortonIndex mindex = blockStartIdx ; mindex < blockEndIdx ; ++mindex){
                    CellClass* cell = (*iterCells)->getCell(mindex);
                    if(cell){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
68
                        FAssertLF(cell->getMortonIndex() == mindex);
69
                        ParticleContainerClass particles = (*iterParticles)->template getLeaf<ParticleContainerClass>(mindex);
70 71 72 73 74 75 76 77 78 79 80
                        FAssertLF(particles.isAttachedToSomething());
                        kernels->P2M(cell, &particles);
                    }
                }
            }

            ++iterParticles;
            ++iterCells;
        }

        FAssertLF(iterParticles == endParticles && iterCells == endCells);
81
        FLOG( FLog::Controller << "\t\t bottomPass in " << timer.tacAndElapsed() << "s\n" );
82 83 84
    }

    void upwardPass(){
85
        FLOG( FTic timer; );
86
        for(int idxLevel = tree->getHeight()-2 ; idxLevel >= 2 ; --idxLevel){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
87 88
            typename OctreeClass::CellGroupIterator iterCells = tree->cellsBegin(idxLevel);
            const typename OctreeClass::CellGroupIterator endCells = tree->cellsEnd(idxLevel);
89

BRAMAS Berenger's avatar
BRAMAS Berenger committed
90 91
            typename OctreeClass::CellGroupIterator iterChildCells = tree->cellsBegin(idxLevel+1);
            const typename OctreeClass::CellGroupIterator endChildCells = tree->cellsEnd(idxLevel+1);
92 93 94 95 96 97 98 99 100

            while(iterCells != endCells && iterChildCells != endChildCells){
                { // Can be a task(in:iterParticles, out:iterChildCells ...)
                    const MortonIndex blockStartIdx = (*iterCells)->getStartingIndex();
                    const MortonIndex blockEndIdx = (*iterCells)->getEndingIndex();

                    for(MortonIndex mindex = blockStartIdx ; mindex < blockEndIdx && iterChildCells != endChildCells; ++mindex){
                        CellClass* cell = (*iterCells)->getCell(mindex);
                        if(cell){
101
                            FAssertLF(cell->getMortonIndex() == mindex);
102
                            CellClass* child[8] = {nullptr,nullptr,nullptr,nullptr,nullptr,nullptr,nullptr,nullptr};
103 104

                            for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
105
                                while(iterChildCells != endChildCells && (*iterChildCells)->getEndingIndex() <= ((mindex<<3)+idxChild) ){
106 107 108 109 110 111
                                    ++iterChildCells;
                                }
                                if( iterChildCells == endChildCells ){
                                    break;
                                }
                                child[idxChild] = (*iterChildCells)->getCell((mindex<<3)+idxChild);
112
                                FAssertLF(child[idxChild] == nullptr || child[idxChild]->getMortonIndex() == ((mindex<<3)+idxChild));
113 114 115 116 117 118 119 120 121 122
                            }

                            kernels->M2M(cell, child, idxLevel);
                        }
                    }
                }

                ++iterCells;
            }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
123 124
            FAssertLF(iterCells == endCells);
            FAssertLF((iterChildCells == endChildCells || (++iterChildCells) == endChildCells));
125 126
            FAssertLF(iterCells == endCells && (iterChildCells == endChildCells || (++iterChildCells) == endChildCells));
        }
127
        FLOG( FLog::Controller << "\t\t upwardPass in " << timer.tacAndElapsed() << "s\n" );
128 129 130
    }

    void transferPass(){
131
        FLOG( FTic timer; );
132
        for(int idxLevel = tree->getHeight()-1 ; idxLevel >= 2 ; --idxLevel){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
133 134
            typename OctreeClass::CellGroupIterator iterCells = tree->cellsBegin(idxLevel);
            const typename OctreeClass::CellGroupIterator endCells = tree->cellsEnd(idxLevel);
135 136 137 138 139 140 141 142 143 144 145

            while(iterCells != endCells){
                std::vector<OutOfBlockInteraction> outsideInteractions;

                { // Can be a task(inout:iterCells, out:outsideInteractions)
                    const MortonIndex blockStartIdx = (*iterCells)->getStartingIndex();
                    const MortonIndex blockEndIdx = (*iterCells)->getEndingIndex();

                    for(MortonIndex mindex = blockStartIdx ; mindex < blockEndIdx ; ++mindex){
                        CellClass* cell = (*iterCells)->getCell(mindex);
                        if(cell){
146
                            FAssertLF(cell->getMortonIndex() == mindex);
147 148
                            MortonIndex interactionsIndexes[189];
                            int interactionsPosition[189];
149
                            const FTreeCoordinate coord(cell->getCoordinate());
150 151
                            int counter = coord.getInteractionNeighbors(idxLevel,interactionsIndexes,interactionsPosition);

152
                            const CellClass* interactions[343];
153 154 155 156 157 158 159
                            memset(interactions, 0, 343*sizeof(CellClass*));
                            int counterExistingCell = 0;

                            for(int idxInter = 0 ; idxInter < counter ; ++idxInter){
                                if( blockStartIdx <= interactionsIndexes[idxInter] && interactionsIndexes[idxInter] < blockEndIdx ){
                                    CellClass* interCell = (*iterCells)->getCell(interactionsIndexes[idxInter]);
                                    if(interCell){
160
                                        FAssertLF(interCell->getMortonIndex() == interactionsIndexes[idxInter]);
161
                                        FAssertLF(interactions[interactionsPosition[idxInter]] == nullptr);
162 163 164 165
                                        interactions[interactionsPosition[idxInter]] = interCell;
                                        counterExistingCell += 1;
                                    }
                                }
166
                                else if(interactionsIndexes[idxInter] < mindex){
167 168 169 170 171 172 173 174 175 176 177 178 179 180 181
                                    OutOfBlockInteraction property;
                                    property.insideIndex = mindex;
                                    property.outIndex    = interactionsIndexes[idxInter];
                                    property.outPosition = interactionsPosition[idxInter];
                                    outsideInteractions.push_back(property);
                                }
                            }

                            kernels->M2L( cell , interactions, counterExistingCell, idxLevel);
                        }
                    }
                }


                // Manage outofblock interaction
182
                FQuickSort<OutOfBlockInteraction, int>::QsSequential(outsideInteractions.data(),int(outsideInteractions.size()));
183

BRAMAS Berenger's avatar
BRAMAS Berenger committed
184
                typename OctreeClass::CellGroupIterator iterLeftCells = tree->cellsBegin(idxLevel);
185
                int currentOutInteraction = 0;
186
                while(iterLeftCells != iterCells && currentOutInteraction < int(outsideInteractions.size())){
187 188 189
                    const MortonIndex blockStartIdx = (*iterLeftCells)->getStartingIndex();
                    const MortonIndex blockEndIdx = (*iterLeftCells)->getEndingIndex();

190
                    while(currentOutInteraction < int(outsideInteractions.size()) && outsideInteractions[currentOutInteraction].outIndex < blockStartIdx){
191 192 193
                        currentOutInteraction += 1;
                    }

194
                    int lastOutInteraction = currentOutInteraction;
195
                    while(lastOutInteraction < int(outsideInteractions.size()) && outsideInteractions[lastOutInteraction].outIndex < blockEndIdx){
196 197 198 199 200 201 202
                        lastOutInteraction += 1;
                    }

                    { // Can be a task(in:currentOutInteraction, in:outsideInteractions, in:lastOutInteraction, inout:iterLeftCells, inout:iterCells)
                        for(int outInterIdx = currentOutInteraction ; outInterIdx < lastOutInteraction ; ++outInterIdx){
                            CellClass* interCell = (*iterLeftCells)->getCell(outsideInteractions[outInterIdx].outIndex);
                            if(interCell){
203
                                FAssertLF(interCell->getMortonIndex() == outsideInteractions[outInterIdx].outIndex);
204 205
                                CellClass* cell = (*iterCells)->getCell(outsideInteractions[outInterIdx].insideIndex);
                                FAssertLF(cell);
206 207
                                FAssertLF(cell->getMortonIndex() == outsideInteractions[outInterIdx].insideIndex);

208
                                const CellClass* interactions[343];
209 210 211 212 213
                                memset(interactions, 0, 343*sizeof(CellClass*));
                                interactions[outsideInteractions[outInterIdx].outPosition] = interCell;
                                const int counter = 1;
                                kernels->M2L( cell , interactions, counter, idxLevel);

214
                                interactions[outsideInteractions[outInterIdx].outPosition] = nullptr;
215 216 217 218 219 220 221 222 223 224 225 226 227 228
                                interactions[getOppositeInterIndex(outsideInteractions[outInterIdx].outPosition)] = cell;
                                kernels->M2L( interCell , interactions, counter, idxLevel);
                            }
                        }
                    }

                    currentOutInteraction = lastOutInteraction;
                    ++iterLeftCells;
                }

                ++iterCells;
            }

        }
229
        FLOG( FLog::Controller << "\t\t transferPass in " << timer.tacAndElapsed() << "s\n" );
230 231 232
    }

    void downardPass(){
233
        FLOG( FTic timer; );
234
        for(int idxLevel = 2 ; idxLevel <= tree->getHeight()-2 ; ++idxLevel){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
235 236
            typename OctreeClass::CellGroupIterator iterCells = tree->cellsBegin(idxLevel);
            const typename OctreeClass::CellGroupIterator endCells = tree->cellsEnd(idxLevel);
237

BRAMAS Berenger's avatar
BRAMAS Berenger committed
238 239
            typename OctreeClass::CellGroupIterator iterChildCells = tree->cellsBegin(idxLevel+1);
            const typename OctreeClass::CellGroupIterator endChildCells = tree->cellsEnd(idxLevel+1);
240 241 242 243 244 245 246 247 248

            while(iterCells != endCells && iterChildCells != endChildCells){
                { // Can be a task(in:iterParticles, inout:iterChildCells ...)
                    const MortonIndex blockStartIdx = (*iterCells)->getStartingIndex();
                    const MortonIndex blockEndIdx = (*iterCells)->getEndingIndex();

                    for(MortonIndex mindex = blockStartIdx ; mindex < blockEndIdx && iterChildCells != endChildCells; ++mindex){
                        CellClass* cell = (*iterCells)->getCell(mindex);
                        if(cell){
249
                            FAssertLF(cell->getMortonIndex() == mindex);
250
                            CellClass* child[8] = {nullptr,nullptr,nullptr,nullptr,nullptr,nullptr,nullptr,nullptr};
251 252

                            for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
253
                                while(iterChildCells != endChildCells && (*iterChildCells)->getEndingIndex() <= ((mindex<<3)+idxChild) ){
254 255 256 257 258 259
                                    ++iterChildCells;
                                }
                                if( iterChildCells == endChildCells ){
                                    break;
                                }
                                child[idxChild] = (*iterChildCells)->getCell((mindex<<3)+idxChild);
260
                                FAssertLF(child[idxChild] == nullptr || child[idxChild]->getMortonIndex() == ((mindex<<3)+idxChild));
261 262 263 264 265 266 267 268 269 270 271 272
                            }

                            kernels->L2L(cell, child, idxLevel);
                        }
                    }
                }

                ++iterCells;
            }

            FAssertLF(iterCells == endCells && (iterChildCells == endChildCells || (++iterChildCells) == endChildCells));
        }
273
        FLOG( FLog::Controller << "\t\t downardPass in " << timer.tacAndElapsed() << "s\n" );
274 275 276
    }

    void directPass(){
277
        FLOG( FTic timer; );
278
        {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
279 280
            typename OctreeClass::ParticleGroupIterator iterParticles = tree->leavesBegin();
            const typename OctreeClass::ParticleGroupIterator endParticles = tree->leavesEnd();
281

BRAMAS Berenger's avatar
BRAMAS Berenger committed
282 283
            typename OctreeClass::CellGroupIterator iterCells = tree->cellsBegin(tree->getHeight()-1);
            const typename OctreeClass::CellGroupIterator endCells = tree->cellsEnd(tree->getHeight()-1);
284 285 286 287

            while(iterParticles != endParticles && iterCells != endCells){
                { // Can be a task(in:iterCells, inout:iterParticles)
                    const MortonIndex blockStartIdx = (*iterCells)->getStartingIndex();
BRAMAS Berenger's avatar
BRAMAS Berenger committed
288
                    const MortonIndex blockEndIdx = (*iterCells)->getEndingIndex();
289 290 291 292

                    for(MortonIndex mindex = blockStartIdx ; mindex < blockEndIdx ; ++mindex){
                        CellClass* cell = (*iterCells)->getCell(mindex);
                        if(cell){
293
                            ParticleContainerClass particles = (*iterParticles)->template getLeaf<ParticleContainerClass>(mindex);
294
                            FAssertLF(particles.isAttachedToSomething());
295
                            kernels->L2P(cell, &particles);
296 297 298 299 300 301 302 303 304 305 306
                        }
                    }
                }

                ++iterParticles;
                ++iterCells;
            }

            FAssertLF(iterParticles == endParticles && iterCells == endCells);
        }
        {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
307 308
            typename OctreeClass::ParticleGroupIterator iterParticles = tree->leavesBegin();
            const typename OctreeClass::ParticleGroupIterator endParticles = tree->leavesEnd();
309 310 311 312 313 314 315 316 317

            while(iterParticles != endParticles){
                typename std::vector<OutOfBlockInteraction> outsideInteractions;

                { // Can be a task(inout:iterCells, out:outsideInteractions)
                    const MortonIndex blockStartIdx = (*iterParticles)->getStartingIndex();
                    const MortonIndex blockEndIdx = (*iterParticles)->getEndingIndex();

                    for(MortonIndex mindex = blockStartIdx ; mindex < blockEndIdx ; ++mindex){
318
                        ParticleContainerClass particles = (*iterParticles)->template getLeaf<ParticleContainerClass>(mindex);
319 320 321 322 323 324 325 326
                        if(particles.isAttachedToSomething()){
                            MortonIndex interactionsIndexes[26];
                            int interactionsPosition[26];
                            FTreeCoordinate coord(mindex, tree->getHeight()-1);
                            int counter = coord.getNeighborsIndexes(tree->getHeight(),interactionsIndexes,interactionsPosition);

                            ParticleContainerClass interactionsObjects[27];
                            ParticleContainerClass* interactions[27];
327
                            memset(interactions, 0, 27*sizeof(ParticleContainerClass*));
328 329 330 331
                            int counterExistingCell = 0;

                            for(int idxInter = 0 ; idxInter < counter ; ++idxInter){
                                if( blockStartIdx <= interactionsIndexes[idxInter] && interactionsIndexes[idxInter] < blockEndIdx ){
332
                                    interactionsObjects[counterExistingCell] = (*iterParticles)->template getLeaf<ParticleContainerClass>(interactionsIndexes[idxInter]);
333
                                    if(interactionsObjects[counterExistingCell].isAttachedToSomething()){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
334
                                        FAssertLF(interactions[interactionsPosition[idxInter]] == nullptr);
335 336 337 338
                                        interactions[interactionsPosition[idxInter]] = &interactionsObjects[counterExistingCell];
                                        counterExistingCell += 1;
                                    }
                                }
339
                                else if(interactionsIndexes[idxInter] < mindex){
340 341 342 343 344 345 346 347 348 349 350 351 352 353 354
                                    OutOfBlockInteraction property;
                                    property.insideIndex = mindex;
                                    property.outIndex    = interactionsIndexes[idxInter];
                                    property.outPosition = interactionsPosition[idxInter];
                                    outsideInteractions.push_back(property);
                                }
                            }

                            kernels->P2P( coord, &particles, &particles , interactions, counterExistingCell);
                        }
                    }
                }


                // Manage outofblock interaction
355
                FQuickSort<OutOfBlockInteraction, int>::QsSequential(outsideInteractions.data(),int(outsideInteractions.size()));
356

BRAMAS Berenger's avatar
BRAMAS Berenger committed
357
                typename OctreeClass::ParticleGroupIterator iterLeftParticles = tree->leavesBegin();
358
                int currentOutInteraction = 0;
359
                while(iterLeftParticles != iterParticles && currentOutInteraction < int(outsideInteractions.size())){
360 361 362
                    const MortonIndex blockStartIdx = (*iterLeftParticles)->getStartingIndex();
                    const MortonIndex blockEndIdx = (*iterLeftParticles)->getEndingIndex();

363
                    while(currentOutInteraction < int(outsideInteractions.size()) && outsideInteractions[currentOutInteraction].outIndex < blockStartIdx){
364 365 366
                        currentOutInteraction += 1;
                    }

367
                    int lastOutInteraction = currentOutInteraction;
368
                    while(lastOutInteraction < int(outsideInteractions.size()) && outsideInteractions[lastOutInteraction].outIndex < blockEndIdx){
369 370 371 372 373
                        lastOutInteraction += 1;
                    }

                    { // Can be a task(in:currentOutInteraction, in:outsideInteractions, in:lastOutInteraction, inout:iterLeftParticles, inout:iterParticles)
                        for(int outInterIdx = currentOutInteraction ; outInterIdx < lastOutInteraction ; ++outInterIdx){
374
                            ParticleContainerClass interParticles = (*iterLeftParticles)->template getLeaf<ParticleContainerClass>(outsideInteractions[outInterIdx].outIndex);
375
                            if(interParticles.isAttachedToSomething()){
376
                                ParticleContainerClass particles = (*iterParticles)->template getLeaf<ParticleContainerClass>(outsideInteractions[outInterIdx].insideIndex);
377
                                FAssertLF(particles.isAttachedToSomething());
378
                                ParticleContainerClass* interactions[27];
379
                                memset(interactions, 0, 27*sizeof(ParticleContainerClass*));
380 381 382 383
                                interactions[outsideInteractions[outInterIdx].outPosition] = &interParticles;
                                const int counter = 1;
                                kernels->P2PRemote( FTreeCoordinate(outsideInteractions[outInterIdx].insideIndex, tree->getHeight()-1), &particles, &particles , interactions, counter);

384
                                interactions[outsideInteractions[outInterIdx].outPosition] = nullptr;
385 386 387 388 389 390 391 392 393 394 395 396 397
                                interactions[getOppositeNeighIndex(outsideInteractions[outInterIdx].outPosition)] = &particles;
                                kernels->P2PRemote( FTreeCoordinate(outsideInteractions[outInterIdx].outIndex, tree->getHeight()-1), &interParticles, &interParticles , interactions, counter);
                            }
                        }
                    }

                    currentOutInteraction = lastOutInteraction;
                    ++iterLeftParticles;
                }

                ++iterParticles;
            }
        }
398
        FLOG( FLog::Controller << "\t\t directPass in " << timer.tacAndElapsed() << "s\n" );
399 400 401 402
    }

    int getOppositeNeighIndex(const int index) const {
        // ((idxX+1)*3 + (idxY+1)) * 3 + (idxZ+1)
403
        return 27-index-1;
404 405 406 407
    }

    int getOppositeInterIndex(const int index) const {
        // ((( (xdiff+3) * 7) + (ydiff+3))) * 7 + zdiff + 3
408
        return 343-index-1;
409 410 411 412 413
    }
};


#endif // FGROUPSEQALGORITHM_HPP