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