src/decompositionMethods/decompositionMethods/metisDecomp/metisDecomp.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 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 "IFstream.H"
  31 #include "Time.H"
  32 #include "cyclicPolyPatch.H"
  33
  34 extern "C"
  35 {
  36 #define OMPI_SKIP_MPICXX
  37 #   include "metis.h"
  38 }
  39
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45     defineTypeNameAndDebug(metisDecomp, 0);
  46
  47     addToRunTimeSelectionTable
  48     (
  49         decompositionMethod,
  50         metisDecomp,
  51         dictionaryMesh
  52     );
  53 }
  54
  55
  56 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  57
  58 // Call Metis with options from dictionary.
  59 Foam::label Foam::metisDecomp::decompose
  60 (
  61     const List<int>& adjncy,
  62     const List<int>& xadj,
  63     const scalarField& cWeights,
  64
  65     List<int>& finalDecomp
  66 )
  67 {
  68     // C style numbering
  69     int numFlag = 0;
  70
  71     // Method of decomposition
  72     // recursive: multi-level recursive bisection (default)
  73     // k-way: multi-level k-way
  74     word method("k-way");
  75
  76     int numCells = xadj.size()-1;
  77
  78     // decomposition options. 0 = use defaults
  79     List<int> options(5, 0);
  80
  81     // processor weights initialised with no size, only used if specified in
  82     // a file
  83     Field<floatScalar> processorWeights;
  84
  85     // cell weights (so on the vertices of the dual)
  86     List<int> cellWeights;
  87
  88     // face weights (so on the edges of the dual)
  89     List<int> faceWeights;
  90
  91
  92     // Check for externally provided cellweights and if so initialise weights
  93     scalar minWeights = gMin(cWeights);
  94     if (cWeights.size() > 0)
  95     {
  96         if (minWeights <= 0)
  97         {
  98             WarningIn
  99             (
 100                 "metisDecomp::decompose"
 101                 "(const pointField&, const scalarField&)"
 102             )   << "Illegal minimum weight " << minWeights
 103                 << endl;
 104         }
 105
 106         if (cWeights.size() != numCells)
 107         {
 108             FatalErrorIn
 109             (
 110                 "metisDecomp::decompose"
 111                 "(const pointField&, const scalarField&)"
 112             )   << "Number of cell weights " << cWeights.size()
 113                 << " does not equal number of cells " << numCells
 114                 << exit(FatalError);
 115         }
 116         // Convert to integers.
 117         cellWeights.setSize(cWeights.size());
 118         forAll(cellWeights, i)
 119         {
 120             cellWeights[i] = int(cWeights[i]/minWeights);
 121         }
 122     }
 123
 124
 125     // Check for user supplied weights and decomp options
 126     if (decompositionDict_.found("metisCoeffs"))
 127     {
 128         const dictionary& metisCoeffs =
 129             decompositionDict_.subDict("metisCoeffs");
 130         word weightsFile;
 131
 132         if (metisCoeffs.readIfPresent("method", method))
 133         {
 134             if (method != "recursive" && method != "k-way")
 135             {
 136                 FatalErrorIn("metisDecomp::decompose()")
 137                     << "Method " << method << " in metisCoeffs in dictionary : "
 138                     << decompositionDict_.name()
 139                     << " should be 'recursive' or 'k-way'"
 140                     << exit(FatalError);
 141             }
 142
 143             Info<< "metisDecomp : Using Metis method     " << method
 144                 << nl << endl;
 145         }
 146
 147         if (metisCoeffs.readIfPresent("options", options))
 148         {
 149             if (options.size() != 5)
 150             {
 151                 FatalErrorIn("metisDecomp::decompose()")
 152                     << "Number of options in metisCoeffs in dictionary : "
 153                     << decompositionDict_.name()
 154                     << " should be 5"
 155                     << exit(FatalError);
 156             }
 157
 158             Info<< "metisDecomp : Using Metis options     " << options
 159                 << nl << endl;
 160         }
 161
 162         if (metisCoeffs.readIfPresent("processorWeights", processorWeights))
 163         {
 164             processorWeights /= sum(processorWeights);
 165
 166             if (processorWeights.size() != nProcessors_)
 167             {
 168                 FatalErrorIn("metisDecomp::decompose(const pointField&)")
 169                     << "Number of processor weights "
 170                     << processorWeights.size()
 171                     << " does not equal number of domains " << nProcessors_
 172                     << exit(FatalError);
 173             }
 174         }
 175
 176         if (metisCoeffs.readIfPresent("cellWeightsFile", weightsFile))
 177         {
 178             Info<< "metisDecomp : Using cell-based weights." << endl;
 179
 180             IOList<int> cellIOWeights
 181             (
 182                 IOobject
 183                 (
 184                     weightsFile,
 185                     mesh_.time().timeName(),
 186                     mesh_,
 187                     IOobject::MUST_READ,
 188                     IOobject::AUTO_WRITE
 189                 )
 190             );
 191             cellWeights.transfer(cellIOWeights);
 192
 193             if (cellWeights.size() != xadj.size()-1)
 194             {
 195                 FatalErrorIn("metisDecomp::decompose(const pointField&)")
 196                     << "Number of cell weights " << cellWeights.size()
 197                     << " does not equal number of cells " << xadj.size()-1
 198                     << exit(FatalError);
 199             }
 200         }
 201     }
 202
 203     int nProcs = nProcessors_;
 204
 205     // output: cell -> processor addressing
 206     finalDecomp.setSize(numCells);
 207
 208     // output: number of cut edges
 209     int edgeCut = 0;
 210
 211     // Vertex weight info
 212     int wgtFlag = 0;
 213     int* vwgtPtr = NULL;
 214     int* adjwgtPtr = NULL;
 215
 216     if (cellWeights.size())
 217     {
 218         vwgtPtr = cellWeights.begin();
 219         wgtFlag += 2;       // Weights on vertices
 220     }
 221     if (faceWeights.size())
 222     {
 223         adjwgtPtr = faceWeights.begin();
 224         wgtFlag += 1;       // Weights on edges
 225     }
 226
 227     if (method == "recursive")
 228     {
 229         if (processorWeights.size())
 230         {
 231             METIS_WPartGraphRecursive
 232             (
 233                 &numCells,         // num vertices in graph
 234                 const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
 235                 const_cast<List<int>&>(adjncy).begin(), // neighbour info
 236                 vwgtPtr,           // vertexweights
 237                 adjwgtPtr,         // no edgeweights
 238                 &wgtFlag,
 239                 &numFlag,
 240                 &nProcs,
 241                 processorWeights.begin(),
 242                 options.begin(),
 243                 &edgeCut,
 244                 finalDecomp.begin()
 245             );
 246         }
 247         else
 248         {
 249             METIS_PartGraphRecursive
 250             (
 251                 &numCells,         // num vertices in graph
 252                 const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
 253                 const_cast<List<int>&>(adjncy).begin(), // neighbour info
 254                 vwgtPtr,           // vertexweights
 255                 adjwgtPtr,         // no edgeweights
 256                 &wgtFlag,
 257                 &numFlag,
 258                 &nProcs,
 259                 options.begin(),
 260                 &edgeCut,
 261                 finalDecomp.begin()
 262             );
 263         }
 264     }
 265     else
 266     {
 267         if (processorWeights.size())
 268         {
 269             METIS_WPartGraphKway
 270             (
 271                 &numCells,         // num vertices in graph
 272                 const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
 273                 const_cast<List<int>&>(adjncy).begin(), // neighbour info
 274                 vwgtPtr,           // vertexweights
 275                 adjwgtPtr,         // no edgeweights
 276                 &wgtFlag,
 277                 &numFlag,
 278                 &nProcs,
 279                 processorWeights.begin(),
 280                 options.begin(),
 281                 &edgeCut,
 282                 finalDecomp.begin()
 283             );
 284         }
 285         else
 286         {
 287             METIS_PartGraphKway
 288             (
 289                 &numCells,         // num vertices in graph
 290                 const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
 291                 const_cast<List<int>&>(adjncy).begin(), // neighbour info
 292                 vwgtPtr,           // vertexweights
 293                 adjwgtPtr,         // no edgeweights
 294                 &wgtFlag,
 295                 &numFlag,
 296                 &nProcs,
 297                 options.begin(),
 298                 &edgeCut,
 299                 finalDecomp.begin()
 300             );
 301         }
 302     }
 303
 304     return edgeCut;
 305 }
 306
 307
 308 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 309
 310 Foam::metisDecomp::metisDecomp
 311 (
 312     const dictionary& decompositionDict,
 313     const polyMesh& mesh
 314 )
 315 :
 316     decompositionMethod(decompositionDict),
 317     mesh_(mesh)
 318 {}
 319
 320
 321 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 322
 323 Foam::labelList Foam::metisDecomp::decompose
 324 (
 325     const pointField& points,
 326     const scalarField& pointWeights
 327 )
 328 {
 329     if (points.size() != mesh_.nCells())
 330     {
 331         FatalErrorIn
 332         (
 333             "metisDecomp::decompose(const pointField&,const scalarField&)"
 334         )   << "Can use this decomposition method only for the whole mesh"
 335             << endl
 336             << "and supply one coordinate (cellCentre) for every cell." << endl
 337             << "The number of coordinates " << points.size() << endl
 338             << "The number of cells in the mesh " << mesh_.nCells()
 339             << exit(FatalError);
 340     }
 341
 342     List<int> adjncy;
 343     List<int> xadj;
 344     calcMetisCSR
 345     (
 346         mesh_,
 347         adjncy,
 348         xadj
 349     );
 350
 351     // Decompose using default weights
 352     List<int> finalDecomp;
 353     decompose(adjncy, xadj, pointWeights, finalDecomp);
 354
 355     // Copy back to labelList
 356     labelList decomp(finalDecomp.size());
 357     forAll(decomp, i)
 358     {
 359         decomp[i] = finalDecomp[i];
 360     }
 361     return decomp;
 362 }
 363
 364
 365 void Foam::metisDecomp::calcMetisCSR
 366 (
 367     const polyMesh& mesh,
 368     List<int>& adjncy,
 369     List<int>& xadj
 370 )
 371 {
 372     // Make Metis CSR (Compressed Storage Format) storage
 373     //   adjncy      : contains neighbours (= edges in graph)
 374     //   xadj(celli) : start of information in adjncy for celli
 375
 376     xadj.setSize(mesh.nCells()+1);
 377
 378     // Initialise the number of internal faces of the cells to twice the
 379     // number of internal faces
 380     label nInternalFaces = 2*mesh.nInternalFaces();
 381
 382     // Check the boundary for coupled patches and add to the number of
 383     // internal faces
 384     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 385
 386     forAll(pbm, patchi)
 387     {
 388         if (isA<cyclicPolyPatch>(pbm[patchi]))
 389         {
 390             nInternalFaces += pbm[patchi].size();
 391         }
 392     }
 393
 394     // Create the adjncy array the size of the total number of internal and
 395     // coupled faces
 396     adjncy.setSize(nInternalFaces);
 397
 398     // Fill in xadj
 399     // ~~~~~~~~~~~~
 400     label freeAdj = 0;
 401
 402     for (label cellI = 0; cellI < mesh.nCells(); cellI++)
 403     {
 404         xadj[cellI] = freeAdj;
 405
 406         const labelList& cFaces = mesh.cells()[cellI];
 407
 408         forAll(cFaces, i)
 409         {
 410             label faceI = cFaces[i];
 411
 412             if
 413             (
 414                 mesh.isInternalFace(faceI)
 415              || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
 416             )
 417             {
 418                 freeAdj++;
 419             }
 420         }
 421     }
 422     xadj[mesh.nCells()] = freeAdj;
 423
 424
 425     // Fill in adjncy
 426     // ~~~~~~~~~~~~~~
 427
 428     labelList nFacesPerCell(mesh.nCells(), 0);
 429
 430     // Internal faces
 431     for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
 432     {
 433         label own = mesh.faceOwner()[faceI];
 434         label nei = mesh.faceNeighbour()[faceI];
 435
 436         adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
 437         adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
 438     }
 439
 440     // Coupled faces. Only cyclics done.
 441     forAll(pbm, patchi)
 442     {
 443         if (isA<cyclicPolyPatch>(pbm[patchi]))
 444         {
 445             const unallocLabelList& faceCells = pbm[patchi].faceCells();
 446
 447             label sizeby2 = faceCells.size()/2;
 448
 449             for (label facei=0; facei<sizeby2; facei++)
 450             {
 451                 label own = faceCells[facei];
 452                 label nei = faceCells[facei + sizeby2];
 453
 454                 adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
 455                 adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
 456             }
 457         }
 458     }
 459 }
 460
 461
 462 // From cell-cell connections to Metis format (like CompactListList)
 463 void Foam::metisDecomp::calcMetisCSR
 464 (
 465     const labelListList& cellCells,
 466     List<int>& adjncy,
 467     List<int>& xadj
 468 )
 469 {
 470     // Count number of internal faces
 471     label nConnections = 0;
 472
 473     forAll(cellCells, coarseI)
 474     {
 475         nConnections += cellCells[coarseI].size();
 476     }
 477
 478     // Create the adjncy array as twice the size of the total number of
 479     // internal faces
 480     adjncy.setSize(nConnections);
 481
 482     xadj.setSize(cellCells.size()+1);
 483
 484
 485     // Fill in xadj
 486     // ~~~~~~~~~~~~
 487     label freeAdj = 0;
 488
 489     forAll(cellCells, coarseI)
 490     {
 491         xadj[coarseI] = freeAdj;
 492
 493         const labelList& cCells = cellCells[coarseI];
 494
 495         forAll(cCells, i)
 496         {
 497             adjncy[freeAdj++] = cCells[i];
 498         }
 499     }
 500     xadj[cellCells.size()] = freeAdj;
 501 }
 502
 503
 504 Foam::labelList Foam::metisDecomp::decompose
 505 (
 506     const labelList& agglom,
 507     const pointField& agglomPoints,
 508     const scalarField& agglomWeights
 509 )
 510 {
 511     if (agglom.size() != mesh_.nCells())
 512     {
 513         FatalErrorIn
 514         (
 515             "metisDecomp::decompose"
 516             "(const labelList&, const pointField&, const scalarField&)"
 517         )   << "Size of cell-to-coarse map " << agglom.size()
 518             << " differs from number of cells in mesh " << mesh_.nCells()
 519             << exit(FatalError);
 520     }
 521
 522     // Make Metis CSR (Compressed Storage Format) storage
 523     //   adjncy      : contains neighbours (= edges in graph)
 524     //   xadj(celli) : start of information in adjncy for celli
 525     List<int> adjncy;
 526     List<int> xadj;
 527     {
 528         // Get cellCells on coarse mesh.
 529         labelListList cellCells;
 530         calcCellCells
 531         (
 532             mesh_,
 533             agglom,
 534             agglomPoints.size(),
 535             cellCells
 536         );
 537
 538         calcMetisCSR(cellCells, adjncy, xadj);
 539     }
 540
 541     // Decompose using default weights
 542     List<int> finalDecomp;
 543     decompose(adjncy, xadj, agglomWeights, finalDecomp);
 544
 545
 546     // Rework back into decomposition for original mesh_
 547     labelList fineDistribution(agglom.size());
 548
 549     forAll(fineDistribution, i)
 550     {
 551         fineDistribution[i] = finalDecomp[agglom[i]];
 552     }
 553
 554     return fineDistribution;
 555 }
 556
 557
 558 Foam::labelList Foam::metisDecomp::decompose
 559 (
 560     const labelListList& globalCellCells,
 561     const pointField& cellCentres,
 562     const scalarField& cellWeights
 563 )
 564 {
 565     if (cellCentres.size() != globalCellCells.size())
 566     {
 567         FatalErrorIn
 568         (
 569             "metisDecomp::decompose"
 570             "(const pointField&, const labelListList&, const scalarField&)"
 571         )   << "Inconsistent number of cells (" << globalCellCells.size()
 572             << ") and number of cell centres (" << cellCentres.size()
 573             << ")." << exit(FatalError);
 574     }
 575
 576
 577     // Make Metis CSR (Compressed Storage Format) storage
 578     //   adjncy      : contains neighbours (= edges in graph)
 579     //   xadj(celli) : start of information in adjncy for celli
 580
 581     List<int> adjncy;
 582     List<int> xadj;
 583     calcMetisCSR(globalCellCells, adjncy, xadj);
 584
 585
 586     // Decompose using default weights
 587     List<int> finalDecomp;
 588     decompose(adjncy, xadj, cellWeights, finalDecomp);
 589
 590     // Copy back to labelList
 591     labelList decomp(finalDecomp.size());
 592     forAll(decomp, i)
 593     {
 594         decomp[i] = finalDecomp[i];
 595     }
 596     return decomp;
 597 }
 598
 599
 600 // ************************************************************************* //