src/decompositionMethods/metisDecomp/metisDecomp.C

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