FOctreeArrangerProc.hpp 14.5 KB
Newer Older
1
// ===================================================================================
2 3 4 5
// Copyright ScalFmm 2016 INRIA, Olivier Coulaud, Bérenger Bramas,
// Matthias Messner olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the
// FMM.
6
//
7
// This software is governed by the CeCILL-C and LGPL licenses and
8
// abiding by the rules of distribution of free software.
9 10 11
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
12
//
13 14 15
// 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
16 17 18
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
19
// ===================================================================================
20 21 22 23 24
#ifndef FOCTREEARRANGERPROC_HPP
#define FOCTREEARRANGERPROC_HPP

#include "../Utils/FGlobal.hpp"
#include "../Containers/FVector.hpp"
25
#include "../Utils/FAssert.hpp"
26 27
#include "../Utils/FMpi.hpp"

28 29
#include "../Utils/FGlobalPeriodic.hpp"

BRAMAS Berenger's avatar
BRAMAS Berenger committed
30 31 32 33
/**
* This example show how to use the FOctreeArrangerProc.
* @example testOctreeRearrangeProc.cpp
*/
34

BRAMAS Berenger's avatar
BRAMAS Berenger committed
35 36
/** @brief This class is an arranger, it move the particles that need to be hosted in a different leaf. This is the parallel version that use MPI.
  *
37 38 39 40
  * For example, if a simulation has been executed and the position
  * of the particles have been changed, then it may be better
  * to move the particles in the tree instead of building a new
  * tree.
41
  */
42
template <class FReal, class OctreeClass, class ContainerClass, class ParticleClass, class ConverterClass >
43
class FOctreeArrangerProc  {
44 45 46
    /** Interval is the min/max morton index
      * for a proc
      */
47 48 49 50 51 52
    struct Interval{
        MortonIndex min;
        MortonIndex max;
    };


53
    /** Find the interval that contains mindex */
54 55
    int getInterval(const MortonIndex mindex, const int size, const Interval intervals[]) const{
        for(int idxProc = 0 ; idxProc < size ; ++idxProc){
56
            // does it contains the index?
57 58 59 60
            if( intervals[idxProc].min <= mindex && mindex < intervals[idxProc].max){
                return idxProc;
            }
        }
61
        // if no interval found return the lastest one
62 63 64 65 66 67 68
        return size - 1;
    }

    OctreeClass* const tree;


public:
69
    /** Basic constructor */
70
    FOctreeArrangerProc(OctreeClass* const inTree) : tree(inTree) {
71
        FAssertLF(tree, "Tree cannot be null");
72 73
    }

74
    /** return false if the tree is empty after processing */
75
    bool rearrange(const FMpi::FComm& comm, const int isPeriodic = DirNone){
76
        // interval of each procs
77 78
        Interval*const intervals = new Interval[comm.processCount()];
        memset(intervals, 0, sizeof(Interval) * comm.processCount());
79 80 81

        {   // We need to exchange interval of each process, this interval
            // will be based on the current morton min max
82 83
            Interval myLastInterval;

84
            // take fist index
85 86 87
            typename OctreeClass::Iterator octreeIterator(tree);
            octreeIterator.gotoBottomLeft();
            myLastInterval.min = octreeIterator.getCurrentGlobalIndex();
88
            // take last index
89 90 91 92
            octreeIterator.gotoRight();
            myLastInterval.max = octreeIterator.getCurrentGlobalIndex();

            // We get the min/max indexes from each procs
93
            FMpi::MpiAssert( MPI_Allgather( &myLastInterval, sizeof(Interval), MPI_BYTE, intervals, sizeof(Interval), MPI_BYTE, comm.getComm()),  __LINE__ );
94

95
            // increase interval in the empty morton index
96 97 98 99 100 101 102 103 104
            intervals[0].min = 0;
            for(int idxProc = 1 ; idxProc < comm.processCount() ; ++idxProc){
                intervals[idxProc].min = ((intervals[idxProc].min - intervals[idxProc-1].max)/2) + intervals[idxProc-1].max;
                intervals[idxProc-1].max = intervals[idxProc].min;
            }

            intervals[comm.processCount() - 1].max = ((1 << (3*(tree->getHeight()-1))) - 1);
        }

105
        // Particles that move
106 107 108
        FVector<ParticleClass>*const toMove = new FVector<ParticleClass>[comm.processCount()];

        { // iterate on the leafs and found particle to remove or to send
109 110
            // For periodic
            const FReal boxWidth = tree->getBoxWidth();
111 112
            const FPoint<FReal> min(tree->getBoxCenter(),-boxWidth/2);
            const FPoint<FReal> max(tree->getBoxCenter(),boxWidth/2);
113

114
            FVector<FSize> indexesToExtract;
115

116 117 118 119
            typename OctreeClass::Iterator octreeIterator(tree);
            octreeIterator.gotoBottomLeft();
            do{
                const MortonIndex currentIndex = octreeIterator.getCurrentGlobalIndex();
120
                ContainerClass* particles = octreeIterator.getCurrentLeaf()->getSrc();
121
                //IdxPart is incremented at the end of the loop
122
                for(FSize idxPart = 0 ; idxPart < particles->getNbParticles(); /*++idxPart*/){
123
                    FPoint<FReal> partPos( particles->getPositions()[0][idxPart],
124 125
                            particles->getPositions()[1][idxPart],
                            particles->getPositions()[2][idxPart] );
126 127
                    // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
                    if( TestPeriodicCondition(isPeriodic, DirPlusX) ){
128
                        while(partPos.getX() >= max.getX()){
129 130
                            partPos.incX(-boxWidth);
                        }
131 132 133 134 135 136 137 138
                    }
                    else if(partPos.getX() >= max.getX()){
                        printf("Error, particle out of Box in +X, index %lld\n", currentIndex);
                        printf("Application is exiting...\n");
                    }
                    if( TestPeriodicCondition(isPeriodic, DirMinusX) ){
                        while(partPos.getX() < min.getX()){
                            partPos.incX(boxWidth);
139
                        }
140 141 142 143 144 145 146
                    }
                    else if(partPos.getX() < min.getX()){
                        printf("Error, particle out of Box in -X, index %lld\n", currentIndex);
                        printf("Application is exiting...\n");
                    }
                    // YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY
                    if( TestPeriodicCondition(isPeriodic, DirPlusY) ){
147
                        while(partPos.getY() >= max.getY()){
148 149
                            partPos.incY(-boxWidth);
                        }
150 151 152 153 154 155 156 157
                    }
                    else if(partPos.getY() >= max.getY()){
                        printf("Error, particle out of Box in +Y, index %lld\n", currentIndex);
                        printf("Application is exiting...\n");
                    }
                    if( TestPeriodicCondition(isPeriodic, DirMinusY) ){
                        while(partPos.getY() < min.getY()){
                            partPos.incY(boxWidth);
158
                        }
159 160 161 162 163 164 165
                    }
                    else if(partPos.getY() < min.getY()){
                        printf("Error, particle out of Box in -Y, index %lld\n", currentIndex);
                        printf("Application is exiting...\n");
                    }
                    // ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
                    if( TestPeriodicCondition(isPeriodic, DirPlusX) ){
166
                        while(partPos.getZ() >= max.getZ()){
167 168 169
                            partPos.incZ(-boxWidth);
                        }
                    }
170 171 172 173 174 175 176 177 178 179 180 181 182 183
                    else if(partPos.getZ() >= max.getZ()){
                        printf("Error, particle out of Box in +Z, index %lld\n", currentIndex);
                        printf("Application is exiting...\n");
                    }
                    if( TestPeriodicCondition(isPeriodic, DirMinusX) ){
                        while(partPos.getZ() < min.getZ()){
                            partPos.incZ(boxWidth);
                        }
                    }
                    else if(partPos.getZ() < min.getZ()){
                        printf("Error, particle out of Box in -Z, index %lld\n", currentIndex);
                        printf("Application is exiting...\n");
                    }
                    // set pos
184 185 186
                    particles->getWPositions()[0][idxPart] = partPos.getX();
                    particles->getWPositions()[1][idxPart] = partPos.getY();
                    particles->getWPositions()[2][idxPart] = partPos.getZ();
187

188
                    const MortonIndex particuleIndex = tree->getMortonFromPosition(partPos);
189 190 191 192
                    // is this particle need to be changed from its leaf
                    if(particuleIndex != currentIndex){
                        // find the right interval
                        const int procConcerned = getInterval( particuleIndex, comm.processCount(), intervals);
193
                        toMove[procConcerned].push(ConverterClass::GetParticleAndRemove(particles,idxPart));
194
                        indexesToExtract.push(idxPart);
195 196 197 198
                        //No need to increment idxPart, since the array has been staggered
                    }
                    else{
                        idxPart++;
199 200
                    }
                }
201 202 203 204

                particles->removeParticles(indexesToExtract.data(), indexesToExtract.getSize());
                indexesToExtract.clear();

205 206 207 208
            } while(octreeIterator.moveRight());
        }

        // To send and recv
209
        ParticleClass* toReceive = nullptr;
210 211
        MPI_Request*const requests = new MPI_Request[comm.processCount()*2];
        memset(requests, 0, sizeof(MPI_Request) * comm.processCount() * 2);
212 213
        FSize*const indexToReceive = new FSize[comm.processCount() + 1];
        memset(indexToReceive, 0, sizeof(FSize) * comm.processCount() + 1);
214 215 216 217 218 219

        int iterRequests = 0;
        int limitRecvSend = 0;
        int hasToRecvFrom = 0;

        { // gather what to send to who + isend data
220 221
            FSize*const counter = new FSize[comm.processCount()];
            memset(counter, 0, sizeof(FSize) * comm.processCount());
222 223 224 225 226

            for(int idxProc = 0 ; idxProc < comm.processCount() ; ++idxProc){
                counter[idxProc] = toMove[idxProc].getSize();
            }
            // say who send to who
227 228 229
            FSize*const allcounter = new FSize[comm.processCount()*comm.processCount()];
            FMpi::MpiAssert( MPI_Allgather( counter, comm.processCount(), FMpi::GetType(*counter), allcounter, comm.processCount(),
                                            FMpi::GetType(*counter), comm.getComm()),  __LINE__ );
230 231

            // prepare buffer to receive
232
            FSize sumToRecv = 0;
233 234 235 236 237 238 239 240 241 242 243 244
            indexToReceive[0] = 0;
            for(int idxProc = 0 ; idxProc < comm.processCount() ; ++idxProc){
                if( idxProc != comm.processId()){
                    sumToRecv += allcounter[idxProc * comm.processCount() + comm.processId()];
                }
                indexToReceive[idxProc + 1] = sumToRecv;
            }
            toReceive = new ParticleClass[sumToRecv];

            // send
            for(int idxProc = 0 ; idxProc < comm.processCount() ; ++idxProc){
                if(idxProc != comm.processId() && allcounter[idxProc * comm.processCount() + comm.processId()]){
245 246
                    FAssertLF( allcounter[idxProc * comm.processCount() + comm.processId()] * sizeof(ParticleClass) < std::numeric_limits<int>::max());
                    FMpi::MpiAssert( MPI_Irecv(&toReceive[indexToReceive[idxProc]], int(allcounter[idxProc * comm.processCount() + comm.processId()] * sizeof(ParticleClass)), MPI_BYTE,
247 248 249 250 251 252 253 254 255 256
                              idxProc, 0, comm.getComm(), &requests[iterRequests++]),  __LINE__ );
                    hasToRecvFrom += 1;
                }
            }

            limitRecvSend = iterRequests;

            // recv
            for(int idxProc = 0 ; idxProc < comm.processCount() ; ++idxProc){
                if(idxProc != comm.processId() && toMove[idxProc].getSize()){
257 258
                    FAssertLF( toMove[idxProc].getSize() * sizeof(ParticleClass) < std::numeric_limits<int>::max());
                    FMpi::MpiAssert( MPI_Isend(toMove[idxProc].data(), int(toMove[idxProc].getSize() * sizeof(ParticleClass)), MPI_BYTE,
259 260 261 262 263 264 265
                              idxProc, 0, comm.getComm(), &requests[iterRequests++]),  __LINE__ );
                }
            }

            delete[] allcounter;
            delete[] counter;
        }
266

267
        { // insert particles that moved
268
            for(FSize idxPart = 0 ; idxPart < toMove[comm.processId()].getSize() ; ++idxPart){
269
                ConverterClass::Insert( tree , toMove[comm.processId()][idxPart]);
270 271
            }
        }
272

273 274 275 276 277
        {   // wait any recv or send
            // if it is a recv then insert particles
            MPI_Status status;
            while( hasToRecvFrom ){
                int done = 0;
278
                FMpi::MpiAssert( MPI_Waitany( iterRequests, requests, &done, &status ),  __LINE__ );
279 280
                if( done < limitRecvSend ){
                    const int source = status.MPI_SOURCE;
281
                    for(FSize idxPart = indexToReceive[source] ; idxPart < indexToReceive[source+1] ; ++idxPart){
282
                        ConverterClass::Insert( tree , toReceive[idxPart]);
283 284 285
                    }
                    hasToRecvFrom -= 1;
                }
286
            }
287
        }
288

289 290 291 292 293 294 295 296
        int counterLeavesAlive = 0;
        { // Remove empty leaves
            typename OctreeClass::Iterator octreeIterator(tree);
            octreeIterator.gotoBottomLeft();
            bool workOnNext = true;

            do{
                // Empty leaf
297
                if( octreeIterator.getCurrentListTargets()->getNbParticles() == 0 ){
298 299 300 301 302 303 304
                    const MortonIndex currentIndex = octreeIterator.getCurrentGlobalIndex();
                    workOnNext = octreeIterator.moveRight();
                    tree->removeLeaf( currentIndex );
                }
                // Not empty, just continue
                else {
                    workOnNext = octreeIterator.moveRight();
305
                    counterLeavesAlive += 1;
306 307 308 309 310
                }
            } while( workOnNext );
        }

        // wait all send
311
        FMpi::MpiAssert( MPI_Waitall( iterRequests, requests, MPI_STATUSES_IGNORE),  __LINE__ );
312 313 314 315 316 317 318 319 320 321 322 323 324 325

        delete[] intervals;
        delete[] toMove;
        delete[] requests;
        delete[] toReceive;
        delete[] indexToReceive;

        // return false if tree is empty
        return counterLeavesAlive != 0;
    }

};

#endif // FOCTREEARRANGERPROC_HPP