src/decompositionMethods/metisDecomp/metisDecomp.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "metisDecomp.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "floatScalar.H"
  29 #include "foamTime.H"
  30 #include "scotchDecomp.H"
  31
  32 extern "C"
  33 {
  34 #define OMPI_SKIP_MPICXX
  35 #   include "metis.h"
  36 }
  37
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43     defineTypeNameAndDebug(metisDecomp, 0);
  44
  45     addToRunTimeSelectionTable
  46     (
  47         decompositionMethod,
  48         metisDecomp,
  49         dictionaryMesh
  50     );
  51 }
  52
  53
  54 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  55
  56 Foam::label Foam::metisDecomp::decompose
  57 (
  58     const List<int>& adjncy,
  59     const List<int>& xadj,
  60     const scalarField& cWeights,
  61     List<int>& finalDecomp
  62 )
  63 {
  64     // Method of decomposition
  65     // recursive: multi-level recursive bisection (default)
  66     // k-way: multi-level k-way
  67     word method("k-way");
  68
  69     int numCells = xadj.size()-1;
  70
  71     // decomposition options. 0 = use defaults
  72     idx_t options[METIS_NOPTIONS];
  73     METIS_SetDefaultOptions(options);
  74
  75     // processor weights initialised with no size, only used if specified in
  76     // a file - use the data type that metis expects here
  77     Field<real_t> processorWeights;
  78
  79     // cell weights (so on the vertices of the dual)
  80     List<int> cellWeights;
  81
  82     // face weights (so on the edges of the dual)
  83     List<int> faceWeights;
  84
  85
  86     // Check for externally provided cellweights and if so initialise weights
  87     scalar minWeights = gMin(cWeights);
  88     if (cWeights.size() > 0)
  89     {
  90         if (minWeights <= 0)
  91         {
  92             WarningIn
  93             (
  94                 "metisDecomp::decompose"
  95                 "(const pointField&, const scalarField&)"
  96             )   << "Illegal minimum weight " << minWeights
  97                 << endl;
  98         }
  99
 100         if (cWeights.size() != numCells)
 101         {
 102             FatalErrorIn
 103             (
 104                 "metisDecomp::decompose"
 105                 "(const pointField&, const scalarField&)"
 106             )   << "Number of cell weights " << cWeights.size()
 107                 << " does not equal number of cells " << numCells
 108                 << exit(FatalError);
 109         }
 110         // Convert to integers.
 111         cellWeights.setSize(cWeights.size());
 112         forAll(cellWeights, i)
 113         {
 114             cellWeights[i] = int(cWeights[i]/minWeights);
 115         }
 116     }
 117
 118
 119     // Check for user supplied weights and decomp options
 120     if (decompositionDict_.found("metisCoeffs"))
 121     {
 122         const dictionary& metisCoeffs =
 123             decompositionDict_.subDict("metisCoeffs");
 124         word weightsFile;
 125
 126         if (metisCoeffs.readIfPresent("method", method))
 127         {
 128             if (method != "recursive" && method != "k-way")
 129             {
 130                 FatalErrorIn("metisDecomp::decompose()")
 131                     << "Method " << method
 132                     << " in metisCoeffs in dictionary : "
 133                     << decompositionDict_.name()
 134                     << " should be 'recursive' or 'k-way'"
 135                     << exit(FatalError);
 136             }
 137
 138             Info<< "metisDecomp : Using Metis method     " << method
 139                 << nl << endl;
 140         }
 141
 142         List<int> mOptions;
 143         if (metisCoeffs.readIfPresent("options", mOptions))
 144         {
 145             if (mOptions.size() != METIS_NOPTIONS)
 146             {
 147                 FatalErrorIn("metisDecomp::decompose()")
 148                     << "Number of options in metisCoeffs in dictionary : "
 149                     << decompositionDict_.name()
 150                     << " should be " << METIS_NOPTIONS
 151                     << exit(FatalError);
 152             }
 153
 154             forAll(mOptions, i)
 155             {
 156                 options[i] = mOptions[i];
 157             }
 158
 159             Info<< "metisDecomp : Using Metis options     " << mOptions
 160                 << nl << endl;
 161         }
 162
 163         if (metisCoeffs.readIfPresent("processorWeights", processorWeights))
 164         {
 165             processorWeights /= sum(processorWeights);
 166
 167             if (processorWeights.size() != nProcessors_)
 168             {
 169                 FatalErrorIn("metisDecomp::decompose(const pointField&)")
 170                     << "Number of processor weights "
 171                     << processorWeights.size()
 172                     << " does not equal number of domains " << nProcessors_
 173                     << exit(FatalError);
 174             }
 175         }
 176
 177         if (metisCoeffs.readIfPresent("cellWeightsFile", weightsFile))
 178         {
 179             Info<< "metisDecomp : Using cell-based weights." << endl;
 180
 181             IOList<int> cellIOWeights
 182             (
 183                 IOobject
 184                 (
 185                     weightsFile,
 186                     mesh_.time().timeName(),
 187                     mesh_,
 188                     IOobject::MUST_READ,
 189                     IOobject::AUTO_WRITE
 190                 )
 191             );
 192             cellWeights.transfer(cellIOWeights);
 193
 194             if (cellWeights.size() != xadj.size()-1)
 195             {
 196                 FatalErrorIn("metisDecomp::decompose(const pointField&)")
 197                     << "Number of cell weights " << cellWeights.size()
 198                     << " does not equal number of cells " << xadj.size()-1
 199                     << exit(FatalError);
 200             }
 201         }
 202     }
 203
 204     int nProcs = nProcessors_;
 205
 206     // output: cell -> processor addressing
 207     finalDecomp.setSize(numCells);
 208
 209     // output: number of cut edges
 210     int edgeCut = 0;
 211
 212     // Vertex weight info
 213     int* vwgtPtr = NULL;
 214     int* adjwgtPtr = NULL;
 215
 216     if (cellWeights.size())
 217     {
 218         vwgtPtr = cellWeights.begin();
 219     }
 220     if (faceWeights.size())
 221     {
 222         adjwgtPtr = faceWeights.begin();
 223     }
 224
 225     int one = 1;
 226
 227     if (method == "recursive")
 228     {
 229         METIS_PartGraphRecursive
 230         (
 231             &numCells,         // num vertices in graph
 232             &one,
 233             const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
 234             const_cast<List<int>&>(adjncy).begin(), // neighbour info
 235             vwgtPtr,           // vertexweights
 236             NULL,
 237             adjwgtPtr,         // no edgeweights
 238             &nProcs,
 239             processorWeights.begin(),
 240             NULL,
 241             options,
 242             &edgeCut,
 243             finalDecomp.begin()
 244         );
 245     }
 246     else
 247     {
 248         METIS_PartGraphKway
 249         (
 250             &numCells,         // num vertices in graph
 251             &one,
 252             const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
 253             const_cast<List<int>&>(adjncy).begin(), // neighbour info
 254             vwgtPtr,           // vertexweights
 255             NULL,
 256             adjwgtPtr,         // no edgeweights
 257             &nProcs,
 258             processorWeights.begin(),
 259             NULL,
 260             options,
 261             &edgeCut,
 262             finalDecomp.begin()
 263         );
 264     }
 265
 266     return edgeCut;
 267 }
 268
 269
 270 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 271
 272 Foam::metisDecomp::metisDecomp
 273 (
 274     const dictionary& decompositionDict,
 275     const polyMesh& mesh
 276 )
 277 :
 278     decompositionMethod(decompositionDict),
 279     mesh_(mesh)
 280 {}
 281
 282
 283 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 284
 285 Foam::labelList Foam::metisDecomp::decompose
 286 (
 287     const pointField& points,
 288     const scalarField& pointWeights
 289 )
 290 {
 291     if (points.size() != mesh_.nCells())
 292     {
 293         FatalErrorIn
 294         (
 295             "metisDecomp::decompose(const pointField&,const scalarField&)"
 296         )   << "Can use this decomposition method only for the whole mesh"
 297             << endl
 298             << "and supply one coordinate (cellCentre) for every cell." << endl
 299             << "The number of coordinates " << points.size() << endl
 300             << "The number of cells in the mesh " << mesh_.nCells()
 301             << exit(FatalError);
 302     }
 303
 304     List<int> adjncy;
 305     List<int> xadj;
 306     calcCSR
 307     (
 308         mesh_,
 309         adjncy,
 310         xadj
 311     );
 312
 313     // Decompose using default weights
 314     List<int> finalDecomp;
 315     decompose(adjncy, xadj, pointWeights, finalDecomp);
 316
 317     // Copy back to labelList
 318     labelList decomp(finalDecomp.size());
 319     forAll(decomp, i)
 320     {
 321         decomp[i] = finalDecomp[i];
 322     }
 323     return decomp;
 324 }
 325
 326
 327 Foam::labelList Foam::metisDecomp::decompose
 328 (
 329     const labelList& fineToCoarse,
 330     const pointField& coarsePoints,
 331     const scalarField& coarseWeights
 332 )
 333 {
 334     if (fineToCoarse.size() != mesh_.nCells())
 335     {
 336         FatalErrorIn
 337         (
 338             "metisDecomp::decompose"
 339             "(const labelList&, const pointField&, const scalarField&)"
 340         )   << "Size of cell-to-coarse map " << fineToCoarse.size()
 341             << " differs from number of cells in mesh " << mesh_.nCells()
 342             << exit(FatalError);
 343     }
 344
 345     // Make Metis CSR (Compressed Storage Format) storage
 346     //   adjncy      : contains neighbours (= edges in graph)
 347     //   xadj(celli) : start of information in adjncy for celli
 348     List<int> adjncy;
 349     List<int> xadj;
 350     {
 351         // Get cellCells on coarse mesh.
 352         labelListList cellCells;
 353         calcCellCells
 354         (
 355             mesh_,
 356             fineToCoarse,
 357             coarsePoints.size(),
 358             cellCells
 359         );
 360
 361         calcCSR(cellCells, adjncy, xadj);
 362     }
 363
 364     // Decompose using default weights
 365     List<int> finalDecomp;
 366     decompose(adjncy, xadj, coarseWeights, finalDecomp);
 367
 368
 369     // Rework back into decomposition for original mesh_
 370     labelList fineDistribution(fineToCoarse.size());
 371
 372     forAll(fineDistribution, i)
 373     {
 374         fineDistribution[i] = finalDecomp[fineToCoarse[i]];
 375     }
 376
 377     return fineDistribution;
 378 }
 379
 380
 381 Foam::labelList Foam::metisDecomp::decompose
 382 (
 383     const labelListList& globalCellCells,
 384     const pointField& cc,
 385     const scalarField& cWeights
 386 )
 387 {
 388     if (cc.size() != globalCellCells.size())
 389     {
 390         FatalErrorIn
 391         (
 392             "metisDecomp::decompose"
 393             "(const pointField&, const labelListList&, const scalarField&)"
 394         )   << "Inconsistent number of cells (" << globalCellCells.size()
 395             << ") and number of cell centres (" << cc.size()
 396             << ")." << exit(FatalError);
 397     }
 398
 399
 400     // Make Metis CSR (Compressed Storage Format) storage
 401     //   adjncy      : contains neighbours (= edges in graph)
 402     //   xadj(celli) : start of information in adjncy for celli
 403
 404     List<int> adjncy;
 405     List<int> xadj;
 406     calcCSR(globalCellCells, adjncy, xadj);
 407
 408     // Decompose using default weights
 409     List<int> finalDecomp;
 410     decompose(adjncy, xadj, cWeights, finalDecomp);
 411
 412     // Copy back to labelList
 413     labelList decomp(finalDecomp.size());
 414
 415     forAll(decomp, i)
 416     {
 417         decomp[i] = finalDecomp[i];
 418     }
 419
 420     return decomp;
 421 }
 422
 423
 424 // ************************************************************************* //