src/parallel/decompose/decompositionMethods/multiLevelDecomp/multiLevelDecomp.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "multiLevelDecomp.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "IFstream.H"
  29 #include "globalIndex.H"
  30 #include "mapDistribute.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36     defineTypeNameAndDebug(multiLevelDecomp, 0);
  37
  38     addToRunTimeSelectionTable
  39     (
  40         decompositionMethod,
  41         multiLevelDecomp,
  42         dictionary
  43     );
  44 }
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 // Given a subset of cells determine the new global indices. The problem
  50 // is in the cells from neighbouring processors which need to be renumbered.
  51 void Foam::multiLevelDecomp::subsetGlobalCellCells
  52 (
  53     const label nDomains,
  54     const label domainI,
  55     const labelList& dist,
  56
  57     const labelListList& cellCells,
  58     const labelList& set,
  59     labelListList& subCellCells,
  60     labelList& cutConnections
  61 ) const
  62 {
  63     // Determine new index for cells by inverting subset
  64     labelList oldToNew(invert(cellCells.size(), set));
  65
  66     globalIndex globalCells(cellCells.size());
  67
  68     // Subset locally the elements for which I have data
  69     subCellCells = UIndirectList<labelList>(cellCells, set);
  70
  71     // Get new indices for neighbouring processors
  72     List<Map<label> > compactMap;
  73     mapDistribute map(globalCells, subCellCells, compactMap);
  74     map.distribute(oldToNew);
  75     labelList allDist(dist);
  76     map.distribute(allDist);
  77
  78     // Now we have:
  79     // oldToNew : the locally-compact numbering of all our cellCells. -1 if
  80     //            cellCell is not in set.
  81     // allDist  : destination domain for all our cellCells
  82     // subCellCells : indexes into oldToNew and allDist
  83
  84     // Globally compact numbering for cells in set.
  85     globalIndex globalSubCells(set.size());
  86
  87     // Now subCellCells contains indices into oldToNew which are the
  88     // new locations of the neighbouring cells.
  89
  90     cutConnections.setSize(nDomains);
  91     cutConnections = 0;
  92
  93     forAll(subCellCells, subCellI)
  94     {
  95         labelList& cCells = subCellCells[subCellI];
  96
  97         // Keep the connections to valid mapped cells
  98         label newI = 0;
  99         forAll(cCells, i)
 100         {
 101             // Get locally-compact cell index of neighbouring cell
 102             label nbrCellI = oldToNew[cCells[i]];
 103             if (nbrCellI == -1)
 104             {
 105                 cutConnections[allDist[cCells[i]]]++;
 106             }
 107             else
 108             {
 109                 // Reconvert local cell index into global one
 110
 111                 // Get original neighbour
 112                 label cellI = set[subCellI];
 113                 label oldNbrCellI = cellCells[cellI][i];
 114                 // Get processor from original neighbour
 115                 label procI = globalCells.whichProcID(oldNbrCellI);
 116                 // Convert into global compact numbering
 117                 cCells[newI++] = globalSubCells.toGlobal(procI, nbrCellI);
 118             }
 119         }
 120         cCells.setSize(newI);
 121     }
 122 }
 123
 124
 125 void Foam::multiLevelDecomp::decompose
 126 (
 127     const labelListList& pointPoints,
 128     const pointField& points,
 129     const scalarField& pointWeights,
 130     const labelList& pointMap,      // map back to original points
 131     const label levelI,
 132
 133     labelField& finalDecomp
 134 )
 135 {
 136     labelList dist
 137     (
 138         methods_[levelI].decompose
 139         (
 140             pointPoints,
 141             points,
 142             pointWeights
 143         )
 144     );
 145
 146     forAll(pointMap, i)
 147     {
 148         label orig = pointMap[i];
 149         finalDecomp[orig] += dist[i];
 150     }
 151
 152     if (levelI != methods_.size()-1)
 153     {
 154         // Recurse
 155
 156         // Determine points per domain
 157         label n = methods_[levelI].nDomains();
 158         labelListList domainToPoints(invertOneToMany(n, dist));
 159
 160         // 'Make space' for new levels of decomposition
 161         finalDecomp *= methods_[levelI+1].nDomains();
 162
 163         // Extract processor+local index from point-point addressing
 164         if (debug && Pstream::master())
 165         {
 166             Pout<< "Decomposition at level " << levelI << " :" << endl;
 167         }
 168
 169         for (label domainI = 0; domainI < n; domainI++)
 170         {
 171             // Extract elements for current domain
 172             const labelList domainPoints(findIndices(dist, domainI));
 173
 174             // Subset point-wise data.
 175             pointField subPoints(points, domainPoints);
 176             scalarField subWeights(pointWeights, domainPoints);
 177             labelList subPointMap(UIndirectList<label>(pointMap, domainPoints));
 178             // Subset point-point addressing (adapt global numbering)
 179             labelListList subPointPoints;
 180             labelList nOutsideConnections;
 181             subsetGlobalCellCells
 182             (
 183                 n,
 184                 domainI,
 185                 dist,
 186
 187                 pointPoints,
 188                 domainPoints,
 189
 190                 subPointPoints,
 191                 nOutsideConnections
 192             );
 193
 194             label nPoints = returnReduce(domainPoints.size(), plusOp<label>());
 195             Pstream::listCombineGather(nOutsideConnections, plusEqOp<label>());
 196             Pstream::listCombineScatter(nOutsideConnections);
 197             label nPatches = 0;
 198             label nFaces = 0;
 199             forAll(nOutsideConnections, i)
 200             {
 201                 if (nOutsideConnections[i] > 0)
 202                 {
 203                     nPatches++;
 204                     nFaces += nOutsideConnections[i];
 205                 }
 206             }
 207
 208             string oldPrefix;
 209             if (debug && Pstream::master())
 210             {
 211                 Pout<< "    Domain " << domainI << nl
 212                     << "        Number of cells = " << nPoints << nl
 213                     << "        Number of inter-domain patches = " << nPatches
 214                     << nl
 215                     << "        Number of inter-domain faces = " << nFaces << nl
 216                     << endl;
 217                 oldPrefix = Pout.prefix();
 218                 Pout.prefix() = "  " + oldPrefix;
 219             }
 220
 221             decompose
 222             (
 223                 subPointPoints,
 224                 subPoints,
 225                 subWeights,
 226                 subPointMap,
 227                 levelI+1,
 228
 229                 finalDecomp
 230             );
 231             if (debug && Pstream::master())
 232             {
 233                 Pout.prefix() = oldPrefix;
 234             }
 235         }
 236
 237
 238         if (debug)
 239         {
 240             // Do straight decompose of two levels
 241             label nNext = methods_[levelI+1].nDomains();
 242             label nTotal = n*nNext;
 243
 244             // Retrieve original level0 dictionary and modify number of domains
 245             dictionary::const_iterator iter =
 246                 decompositionDict_.subDict(typeName + "Coeffs").begin();
 247             dictionary myDict = iter().dict();
 248             myDict.set("numberOfSubdomains", nTotal);
 249
 250             if (debug && Pstream::master())
 251             {
 252                 Pout<< "Reference decomposition with " << myDict << " :"
 253                     << endl;
 254             }
 255
 256             autoPtr<decompositionMethod> method0 = decompositionMethod::New
 257             (
 258                 myDict
 259             );
 260             labelList dist
 261             (
 262                 method0().decompose
 263                 (
 264                     pointPoints,
 265                     points,
 266                     pointWeights
 267                 )
 268             );
 269
 270             for (label blockI = 0; blockI < n; blockI++)
 271             {
 272                 // Count the number inbetween blocks of nNext size
 273
 274                 label nPoints = 0;
 275                 labelList nOutsideConnections(n, 0);
 276                 forAll(pointPoints, pointI)
 277                 {
 278                     if ((dist[pointI] / nNext) == blockI)
 279                     {
 280                         nPoints++;
 281
 282                         const labelList& pPoints = pointPoints[pointI];
 283
 284                         forAll(pPoints, i)
 285                         {
 286                             label distBlockI = dist[pPoints[i]] / nNext;
 287                             if (distBlockI != blockI)
 288                             {
 289                                 nOutsideConnections[distBlockI]++;
 290                             }
 291                         }
 292                     }
 293                 }
 294
 295                 reduce(nPoints, plusOp<label>());
 296                 Pstream::listCombineGather
 297                 (
 298                     nOutsideConnections,
 299                     plusEqOp<label>()
 300                 );
 301                 Pstream::listCombineScatter(nOutsideConnections);
 302                 label nPatches = 0;
 303                 label nFaces = 0;
 304                 forAll(nOutsideConnections, i)
 305                 {
 306                     if (nOutsideConnections[i] > 0)
 307                     {
 308                         nPatches++;
 309                         nFaces += nOutsideConnections[i];
 310                     }
 311                 }
 312
 313                 if (debug && Pstream::master())
 314                 {
 315                     Pout<< "    Domain " << blockI << nl
 316                         << "        Number of cells = " << nPoints << nl
 317                         << "        Number of inter-domain patches = "
 318                         << nPatches << nl
 319                         << "        Number of inter-domain faces = " << nFaces
 320                         << nl << endl;
 321                 }
 322             }
 323         }
 324     }
 325 }
 326
 327
 328 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 329
 330 Foam::multiLevelDecomp::multiLevelDecomp(const dictionary& decompositionDict)
 331 :
 332     decompositionMethod(decompositionDict),
 333     methodsDict_(decompositionDict_.subDict(typeName + "Coeffs"))
 334 {
 335     methods_.setSize(methodsDict_.size());
 336     label i = 0;
 337     forAllConstIter(dictionary, methodsDict_, iter)
 338     {
 339         methods_.set(i++, decompositionMethod::New(iter().dict()));
 340     }
 341
 342     label n = 1;
 343     Info<< "decompositionMethod " << type() << " :" << endl;
 344     forAll(methods_, i)
 345     {
 346         Info<< "    level " << i << " decomposing with " << methods_[i].type()
 347             << " into " << methods_[i].nDomains() << " subdomains." << endl;
 348
 349         n *= methods_[i].nDomains();
 350     }
 351
 352     if (n != nDomains())
 353     {
 354         FatalErrorIn("multiLevelDecomp::multiLevelDecomp(const dictionary&)")
 355             << "Top level decomposition specifies " << nDomains()
 356             << " domains which is not equal to the product of"
 357             << " all sub domains " << n
 358             << exit(FatalError);
 359     }
 360 }
 361
 362
 363 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 364
 365 bool Foam::multiLevelDecomp::parallelAware() const
 366 {
 367     forAll(methods_, i)
 368     {
 369         if (!methods_[i].parallelAware())
 370         {
 371             return false;
 372         }
 373     }
 374     return true;
 375 }
 376
 377
 378 Foam::labelList Foam::multiLevelDecomp::decompose
 379 (
 380     const polyMesh& mesh,
 381     const pointField& cc,
 382     const scalarField& cWeights
 383 )
 384 {
 385     CompactListList<label> cellCells;
 386     calcCellCells(mesh, identity(cc.size()), cc.size(), cellCells);
 387
 388     labelField finalDecomp(cc.size(), 0);
 389     labelList cellMap(identity(cc.size()));
 390
 391     decompose
 392     (
 393         cellCells(),
 394         cc,
 395         cWeights,
 396         cellMap,      // map back to original cells
 397         0,
 398
 399         finalDecomp
 400     );
 401
 402     return finalDecomp;
 403 }
 404
 405
 406 Foam::labelList Foam::multiLevelDecomp::decompose
 407 (
 408     const labelListList& globalPointPoints,
 409     const pointField& points,
 410     const scalarField& pointWeights
 411 )
 412 {
 413     labelField finalDecomp(points.size(), 0);
 414     labelList pointMap(identity(points.size()));
 415
 416     decompose
 417     (
 418         globalPointPoints,
 419         points,
 420         pointWeights,
 421         pointMap,       // map back to original points
 422         0,
 423
 424         finalDecomp
 425     );
 426
 427     return finalDecomp;
 428 }
 429
 430
 431 // ************************************************************************* //