src/decompositionMethods/metisDecomp/metisDecomp.C

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