src/meshLibrary/utilities/octrees/meshOctree/meshOctreeModifier/meshOctreeModifierLoadDistribution.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "meshOctreeModifier.H"
  29
  30 #include <map>
  31
  32 # ifdef USE_OMP
  33 #include <omp.h>
  34 # endif
  35
  36 //#define OCTREETiming
  37 //#define DEBUGBalancing
  38
  39 # ifdef DEBUGBalancing
  40 #include <sstream>
  41 # endif
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 namespace Foam
  46 {
  47
  48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49 // Private member functions
  50
  51 void meshOctreeModifier::loadDistribution(const direction usedType)
  52 {
  53     if( octree_.neiProcs().size() == 0 )
  54         return;
  55
  56     # ifdef OCTREETiming
  57     returnReduce(1, sumOp<label>());
  58     const scalar startTime = omp_get_wtime();
  59
  60     returnReduce(1, sumOp<label>());
  61     const scalar t1 = omp_get_wtime();
  62     Info << "Creation of list of leaves lasted " << t1-startTime << endl;
  63     # endif
  64
  65     const LongList<meshOctreeCube*>& leaves = octree_.leaves_;
  66
  67     label localNumWeighs(0);
  68     labelList globalLeafWeight(leaves.size());
  69
  70     if( usedType )
  71     {
  72         forAll(leaves, leafI)
  73         {
  74             if( leaves[leafI]->cubeType() & usedType )
  75             {
  76                 globalLeafWeight[leafI] = localNumWeighs;
  77                 ++localNumWeighs;
  78             }
  79             else
  80             {
  81                 globalLeafWeight[leafI] = -1;
  82             }
  83         }
  84     }
  85     else
  86     {
  87         forAll(leaves, leafI)
  88         {
  89             globalLeafWeight[leafI] = localNumWeighs;
  90             ++localNumWeighs;
  91         }
  92     }
  93
  94     # ifdef OCTREETiming
  95     returnReduce(1, sumOp<label>());
  96     const scalar t2 = omp_get_wtime();
  97     Info << "Creation of global leaf weights lasted " << t2-t1 << endl;
  98     # endif
  99
 100     const label totalNumWeights = returnReduce(localNumWeighs, sumOp<label>());
 101     const label nWeightsPerProcessor = totalNumWeights / Pstream::nProcs();
 102
 103     //- check if balancing should be performed
 104     //- the tolerance is set to 5% difference in the number of boxes
 105     //- from the ideal one
 106     label doBalancing(0);
 107     if(
 108         mag
 109         (
 110             scalar(localNumWeighs - nWeightsPerProcessor) /
 111             nWeightsPerProcessor
 112         ) > 0.05
 113     )
 114         doBalancing = 1;
 115
 116     reduce(doBalancing, maxOp<label>());
 117
 118     if( doBalancing == 0 )
 119         return;
 120
 121     Info << "Distributing load between processors" << endl;
 122
 123     //- start calculating new partitions
 124     //- find global labels of the leaf boxes
 125     doBalancing = 0;
 126
 127     labelList procWeights(Pstream::nProcs());
 128     procWeights[Pstream::myProcNo()] = localNumWeighs;
 129     Pstream::gatherList(procWeights);
 130     Pstream::scatterList(procWeights);
 131
 132     for(label procI=0;procI<Pstream::myProcNo();++procI)
 133         doBalancing += procWeights[procI];
 134
 135     forAll(globalLeafWeight, lI)
 136     {
 137         if( globalLeafWeight[lI] != -1 )
 138             globalLeafWeight[lI] += doBalancing;
 139     }
 140
 141     //- leaf boxes which are not in the range for the current processor
 142     //- shall be migrated to other processors
 143     std::map<label, labelLongList> leavesToSend;
 144
 145     bool oneRemainingBox(false);
 146     forAll(globalLeafWeight, leafI)
 147     {
 148         if( globalLeafWeight[leafI] == -1 )
 149             continue;
 150         if( !oneRemainingBox && (leafI == leaves.size() -1) )
 151             continue;
 152
 153         const label newProc =
 154             Foam::min
 155             (
 156                 globalLeafWeight[leafI] / nWeightsPerProcessor,
 157                 Pstream::nProcs()-1
 158             );
 159
 160         if( newProc != Pstream::myProcNo() )
 161         {
 162             leavesToSend[newProc].append(leafI);
 163             leaves[leafI]->setProcNo(newProc);
 164
 165             # ifdef DEBUGBalancing
 166             if( leaves[leafI]->hasContainedElements() )
 167                 Serr << Pstream::myProcNo() << "Deleting a DATA cube "
 168                 << leaves[leafI]->coordinates() << " data is "
 169                 << leaves[leafI]->containedElements() << endl;
 170             # endif
 171         }
 172         else
 173         {
 174             oneRemainingBox = true;
 175         }
 176     }
 177
 178     # ifdef OCTREETiming
 179     returnReduce(1, sumOp<label>());
 180     const scalar t3 = omp_get_wtime();
 181     Info << "Completed assignment of leaves to processors in " << t3-t2 << endl;
 182     # endif
 183
 184     //- send the information to other processors
 185     //- all processors shall received a list containing the same information
 186     //- each processor informs which other processors shall receive data from
 187     //- that processor
 188     labelListList sendToProcesssors(Pstream::nProcs());
 189     sendToProcesssors[Pstream::myProcNo()].setSize(leavesToSend.size());
 190     label counter(0);
 191     for
 192     (
 193         std::map<label, labelLongList>::const_iterator it=leavesToSend.begin();
 194         it!=leavesToSend.end();
 195         ++it
 196     )
 197         sendToProcesssors[Pstream::myProcNo()][counter++] = it->first;
 198
 199     Pstream::gatherList(sendToProcesssors);
 200     Pstream::scatterList(sendToProcesssors);
 201
 202     labelHashSet receiveFrom;
 203     forAll(sendToProcesssors, procI)
 204         forAll(sendToProcesssors[procI], neiI)
 205             if( sendToProcesssors[procI][neiI] == Pstream::myProcNo() )
 206                 receiveFrom.insert(procI);
 207
 208     //- send the coordinates of the boxes to other processors
 209     const labelList& sendToProcs = sendToProcesssors[Pstream::myProcNo()];
 210     forAll(sendToProcs, i)
 211     {
 212         const label procI = sendToProcs[i];
 213
 214         List<meshOctreeCubeBasic> sendCoordinates
 215         (
 216             leavesToSend[procI].size()
 217         );
 218
 219         forAll(leavesToSend[procI], lI)
 220         {
 221             const meshOctreeCube& oc = *leaves[leavesToSend[procI][lI]];
 222             sendCoordinates[lI] =
 223                 meshOctreeCubeBasic
 224                 (
 225                     oc.coordinates(),
 226                     oc.cubeType()
 227                 );
 228         }
 229
 230         OPstream toOtherProc
 231         (
 232             Pstream::blocking,
 233             procI,
 234             sendCoordinates.byteSize()
 235         );
 236
 237         toOtherProc << sendCoordinates;
 238     }
 239
 240     //- receive data sent from other processors
 241     LongList<meshOctreeCubeBasic> migratedCubes;
 242     forAllConstIter(labelHashSet, receiveFrom, iter)
 243     {
 244         List<meshOctreeCubeBasic> mc;
 245
 246         IPstream fromOtherProc(Pstream::blocking, iter.key());
 247
 248         fromOtherProc >> mc;
 249
 250         label currSize = migratedCubes.size();
 251         migratedCubes.setSize(currSize+mc.size());
 252         forAll(mc, mcI)
 253         {
 254             migratedCubes[currSize] = mc[mcI];
 255             ++currSize;
 256         }
 257     }
 258
 259     # ifdef OCTREETiming
 260     returnReduce(1, sumOp<label>());
 261     const scalar t4 = omp_get_wtime();
 262     Info << "Data exchange lasted " << t4-t3 << endl;
 263     # endif
 264
 265     //- delete cubes which have been moved to other processors
 266     octree_.initialCubePtr_->purgeProcessorCubes(Pstream::myProcNo());
 267
 268     # ifdef OCTREETiming
 269     returnReduce(1, sumOp<label>());
 270     const scalar t5 = omp_get_wtime();
 271     Info << "Purging lasted " << t5-t4 << endl;
 272     # endif
 273
 274     //- create boxes from the received coordinates
 275     forAll(migratedCubes, mcI)
 276     {
 277         refineTreeForCoordinates
 278         (
 279             migratedCubes[mcI].coordinates(),
 280             Pstream::myProcNo(),
 281             migratedCubes[mcI].cubeType()
 282         );
 283     }
 284
 285     createListOfLeaves();
 286
 287     # ifdef OCTREETiming
 288     returnReduce(1, sumOp<label>());
 289     const scalar t6 = omp_get_wtime();
 290     Info << "Tree refinement lasted " << t6-t5 << endl;
 291     # endif
 292
 293     //- update the communication pattern
 294     updateCommunicationPattern();
 295
 296     # ifdef OCTREETiming
 297     returnReduce(1, sumOp<label>());
 298     const scalar endTime = omp_get_wtime();
 299     Info << "Updating of communication pattern lasted " << endTime-t6 << endl;
 300     Info << "Time for load balancing is " << endTime-startTime << endl;
 301     # endif
 302
 303     Info << "Finished distributing load between processors" << endl;
 304 }
 305
 306 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 307
 308 } // End namespace Foam
 309
 310 // ************************************************************************* //