src/decompositionMethods/parMetisDecomp/parMetisDecomp.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 "parMetisDecomp.H"
  27 #include "metisDecomp.H"
  28 #include "scotchDecomp.H"
  29 #include "syncTools.H"
  30 #include "addToRunTimeSelectionTable.H"
  31 #include "floatScalar.H"
  32 #include "polyMesh.H"
  33 #include "Time.H"
  34 #include "labelIOField.H"
  35 #include "globalIndex.H"
  36
  37 #include <mpi.h>
  38
  39 extern "C"
  40 {
  41 #   include "parmetis.h"
  42 }
  43
  44 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  45
  46 namespace Foam
  47 {
  48     defineTypeNameAndDebug(parMetisDecomp, 0);
  49
  50     addToRunTimeSelectionTable
  51     (
  52         decompositionMethod,
  53         parMetisDecomp,
  54         dictionaryMesh
  55     );
  56 }
  57
  58
  59 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  60
  61 //- Does prevention of 0 cell domains and calls parmetis.
  62 Foam::label Foam::parMetisDecomp::decompose
  63 (
  64     Field<int>& xadj,
  65     Field<int>& adjncy,
  66     const pointField& cellCentres,
  67     Field<int>& cellWeights,
  68     Field<int>& faceWeights,
  69     const List<int>& options,
  70
  71     List<int>& finalDecomp
  72 )
  73 {
  74     // C style numbering
  75     int numFlag = 0;
  76
  77     // Number of dimensions
  78     int nDims = 3;
  79
  80
  81     if (cellCentres.size() != xadj.size()-1)
  82     {
  83         FatalErrorIn("parMetisDecomp::decompose(..)")
  84             << "cellCentres:" << cellCentres.size()
  85             << " xadj:" << xadj.size()
  86             << abort(FatalError);
  87     }
  88
  89
  90     // Get number of cells on all processors
  91     List<int> nLocalCells(Pstream::nProcs());
  92     nLocalCells[Pstream::myProcNo()] = xadj.size()-1;
  93     Pstream::gatherList(nLocalCells);
  94     Pstream::scatterList(nLocalCells);
  95
  96     // Get cell offsets.
  97     List<int> cellOffsets(Pstream::nProcs()+1);
  98     int nGlobalCells = 0;
  99     forAll(nLocalCells, procI)
 100     {
 101         cellOffsets[procI] = nGlobalCells;
 102         nGlobalCells += nLocalCells[procI];
 103     }
 104     cellOffsets[Pstream::nProcs()] = nGlobalCells;
 105
 106     // Convert pointField into float
 107     Field<floatScalar> xyz(3*cellCentres.size());
 108     int compI = 0;
 109     forAll(cellCentres, cellI)
 110     {
 111         const point& cc = cellCentres[cellI];
 112         xyz[compI++] = float(cc.x());
 113         xyz[compI++] = float(cc.y());
 114         xyz[compI++] = float(cc.z());
 115     }
 116
 117     // Make sure every domain has at least one cell
 118     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 119     // (Metis falls over with zero sized domains)
 120     // Trickle cells from processors that have them up to those that
 121     // don't.
 122
 123
 124     // Number of cells to send to the next processor
 125     // (is same as number of cells next processor has to receive)
 126     List<int> nSendCells(Pstream::nProcs(), 0);
 127
 128     for (label procI = nLocalCells.size()-1; procI >=1; procI--)
 129     {
 130         if (nLocalCells[procI]-nSendCells[procI] < 1)
 131         {
 132             nSendCells[procI-1] = nSendCells[procI]-nLocalCells[procI]+1;
 133         }
 134     }
 135
 136     // First receive (so increasing the sizes of all arrays)
 137
 138     if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
 139     {
 140         // Receive cells from previous processor
 141         IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
 142
 143         Field<int> prevXadj(fromPrevProc);
 144         Field<int> prevAdjncy(fromPrevProc);
 145         Field<floatScalar> prevXyz(fromPrevProc);
 146         Field<int> prevCellWeights(fromPrevProc);
 147         Field<int> prevFaceWeights(fromPrevProc);
 148
 149         if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1])
 150         {
 151             FatalErrorIn("parMetisDecomp::decompose(..)")
 152                 << "Expected from processor " << Pstream::myProcNo()-1
 153                 << " connectivity for " << nSendCells[Pstream::myProcNo()-1]
 154                 << " nCells but only received " << prevXadj.size()
 155                 << abort(FatalError);
 156         }
 157
 158         // Insert adjncy
 159         prepend(prevAdjncy, adjncy);
 160         // Adapt offsets and prepend xadj
 161         xadj += prevAdjncy.size();
 162         prepend(prevXadj, xadj);
 163         // Coords
 164         prepend(prevXyz, xyz);
 165         // Weights
 166         prepend(prevCellWeights, cellWeights);
 167         prepend(prevFaceWeights, faceWeights);
 168     }
 169
 170
 171     // Send to my next processor
 172
 173     if (nSendCells[Pstream::myProcNo()] > 0)
 174     {
 175         // Send cells to next processor
 176         OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1);
 177
 178         int nCells = nSendCells[Pstream::myProcNo()];
 179         int startCell = xadj.size()-1 - nCells;
 180         int startFace = xadj[startCell];
 181         int nFaces = adjncy.size()-startFace;
 182
 183         // Send for all cell data: last nCells elements
 184         // Send for all face data: last nFaces elements
 185         toNextProc
 186             << Field<int>::subField(xadj, nCells, startCell)-startFace
 187             << Field<int>::subField(adjncy, nFaces, startFace)
 188             << SubField<floatScalar>(xyz, nDims*nCells, nDims*startCell)
 189             <<
 190             (
 191                 cellWeights.size()
 192               ? static_cast<const Field<int>&>
 193                 (
 194                     Field<int>::subField(cellWeights, nCells, startCell)
 195                 )
 196               : Field<int>(0)
 197             )
 198             <<
 199             (
 200                 faceWeights.size()
 201               ? static_cast<const Field<int>&>
 202                 (
 203                     Field<int>::subField(faceWeights, nFaces, startFace)
 204                 )
 205               : Field<int>(0)
 206             );
 207
 208         // Remove data that has been sent
 209         if (faceWeights.size())
 210         {
 211             faceWeights.setSize(faceWeights.size()-nFaces);
 212         }
 213         if (cellWeights.size())
 214         {
 215             cellWeights.setSize(cellWeights.size()-nCells);
 216         }
 217         xyz.setSize(xyz.size()-nDims*nCells);
 218         adjncy.setSize(adjncy.size()-nFaces);
 219         xadj.setSize(xadj.size() - nCells);
 220     }
 221
 222
 223
 224     // Adapt number of cells
 225     forAll(nSendCells, procI)
 226     {
 227         // Sent cells
 228         nLocalCells[procI] -= nSendCells[procI];
 229
 230         if (procI >= 1)
 231         {
 232             // Received cells
 233             nLocalCells[procI] += nSendCells[procI-1];
 234         }
 235     }
 236     // Adapt cellOffsets
 237     nGlobalCells = 0;
 238     forAll(nLocalCells, procI)
 239     {
 240         cellOffsets[procI] = nGlobalCells;
 241         nGlobalCells += nLocalCells[procI];
 242     }
 243
 244
 245     if (nLocalCells[Pstream::myProcNo()] != (xadj.size()-1))
 246     {
 247         FatalErrorIn("parMetisDecomp::decompose(..)")
 248             << "Have connectivity for " << xadj.size()-1
 249             << " cells but nLocalCells:" << nLocalCells[Pstream::myProcNo()]
 250             << abort(FatalError);
 251     }
 252
 253     // Weight info
 254     int wgtFlag = 0;
 255     int* vwgtPtr = NULL;
 256     int* adjwgtPtr = NULL;
 257
 258     if (cellWeights.size())
 259     {
 260         vwgtPtr = cellWeights.begin();
 261         wgtFlag += 2;       // Weights on vertices
 262     }
 263     if (faceWeights.size())
 264     {
 265         adjwgtPtr = faceWeights.begin();
 266         wgtFlag += 1;       // Weights on edges
 267     }
 268
 269
 270     // Number of weights or balance constraints
 271     int nCon = 1;
 272     // Per processor, per constraint the weight
 273     Field<floatScalar> tpwgts(nCon*nProcessors_, 1./nProcessors_);
 274     // Imbalance tolerance
 275     Field<floatScalar> ubvec(nCon, 1.02);
 276     if (nProcessors_ == 1)
 277     {
 278         // If only one processor there is no imbalance.
 279         ubvec[0] = 1;
 280     }
 281
 282     MPI_Comm comm = MPI_COMM_WORLD;
 283
 284     // output: cell -> processor addressing
 285     finalDecomp.setSize(nLocalCells[Pstream::myProcNo()]);
 286
 287     // output: number of cut edges
 288     int edgeCut = 0;
 289
 290
 291     ParMETIS_V3_PartGeomKway
 292     (
 293         cellOffsets.begin(),    // vtxDist
 294         xadj.begin(),
 295         adjncy.begin(),
 296         vwgtPtr,                // vertexweights
 297         adjwgtPtr,              // edgeweights
 298         &wgtFlag,
 299         &numFlag,
 300         &nDims,
 301         xyz.begin(),
 302         &nCon,
 303         &nProcessors_,          // nParts
 304         tpwgts.begin(),
 305         ubvec.begin(),
 306         const_cast<List<int>&>(options).begin(),
 307         &edgeCut,
 308         finalDecomp.begin(),
 309         &comm
 310     );
 311
 312
 313     // If we sent cells across make sure we undo it
 314     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 315
 316     // Receive back from next processor if I sent something
 317     if (nSendCells[Pstream::myProcNo()] > 0)
 318     {
 319         IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1);
 320
 321         List<int> nextFinalDecomp(fromNextProc);
 322
 323         if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()])
 324         {
 325             FatalErrorIn("parMetisDecomp::decompose(..)")
 326                 << "Expected from processor " << Pstream::myProcNo()+1
 327                 << " decomposition for " << nSendCells[Pstream::myProcNo()]
 328                 << " nCells but only received " << nextFinalDecomp.size()
 329                 << abort(FatalError);
 330         }
 331
 332         append(nextFinalDecomp, finalDecomp);
 333     }
 334
 335     // Send back to previous processor.
 336     if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
 337     {
 338         OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
 339
 340         int nToPrevious = nSendCells[Pstream::myProcNo()-1];
 341
 342         toPrevProc <<
 343             SubList<int>
 344             (
 345                 finalDecomp,
 346                 nToPrevious,
 347                 finalDecomp.size()-nToPrevious
 348             );
 349
 350         // Remove locally what has been sent
 351         finalDecomp.setSize(finalDecomp.size()-nToPrevious);
 352     }
 353
 354     return edgeCut;
 355 }
 356
 357
 358 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 359
 360 Foam::parMetisDecomp::parMetisDecomp
 361 (
 362     const dictionary& decompositionDict,
 363     const polyMesh& mesh
 364 )
 365 :
 366     decompositionMethod(decompositionDict),
 367     mesh_(mesh)
 368 {}
 369
 370
 371 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 372
 373 Foam::labelList Foam::parMetisDecomp::decompose
 374 (
 375     const pointField& cc,
 376     const scalarField& cWeights
 377 )
 378 {
 379     if (cc.size() != mesh_.nCells())
 380     {
 381         FatalErrorIn
 382         (
 383             "parMetisDecomp::decompose"
 384             "(const pointField&, const scalarField&)"
 385         )   << "Can use this decomposition method only for the whole mesh"
 386             << endl
 387             << "and supply one coordinate (cellCentre) for every cell." << endl
 388             << "The number of coordinates " << cc.size() << endl
 389             << "The number of cells in the mesh " << mesh_.nCells()
 390             << exit(FatalError);
 391     }
 392
 393     // For running sequential ...
 394     if (Pstream::nProcs() <= 1)
 395     {
 396         return metisDecomp(decompositionDict_, mesh_).decompose
 397         (
 398             cc,
 399             cWeights
 400         );
 401     }
 402
 403
 404     // Connections
 405     Field<int> adjncy;
 406     // Offsets into adjncy
 407     Field<int> xadj;
 408     calcMetisDistributedCSR
 409     (
 410         mesh_,
 411         adjncy,
 412         xadj
 413     );
 414
 415
 416     // decomposition options. 0 = use defaults
 417     List<int> options(3, 0);
 418     //options[0] = 1;     // don't use defaults but use values below
 419     //options[1] = -1;    // full debug info
 420     //options[2] = 15;    // random number seed
 421
 422     // cell weights (so on the vertices of the dual)
 423     Field<int> cellWeights;
 424
 425     // face weights (so on the edges of the dual)
 426     Field<int> faceWeights;
 427
 428
 429     // Check for externally provided cellweights and if so initialise weights
 430     scalar minWeights = gMin(cWeights);
 431     if (cWeights.size() > 0)
 432     {
 433         if (minWeights <= 0)
 434         {
 435             WarningIn
 436             (
 437                 "metisDecomp::decompose"
 438                 "(const pointField&, const scalarField&)"
 439             )   << "Illegal minimum weight " << minWeights
 440                 << endl;
 441         }
 442
 443         if (cWeights.size() != mesh_.nCells())
 444         {
 445             FatalErrorIn
 446             (
 447                 "parMetisDecomp::decompose"
 448                 "(const pointField&, const scalarField&)"
 449             )   << "Number of cell weights " << cWeights.size()
 450                 << " does not equal number of cells " << mesh_.nCells()
 451                 << exit(FatalError);
 452         }
 453
 454         // Convert to integers.
 455         cellWeights.setSize(cWeights.size());
 456         forAll(cellWeights, i)
 457         {
 458             cellWeights[i] = int(cWeights[i]/minWeights);
 459         }
 460     }
 461
 462
 463     // Check for user supplied weights and decomp options
 464     if (decompositionDict_.found("metisCoeffs"))
 465     {
 466         const dictionary& metisCoeffs =
 467             decompositionDict_.subDict("metisCoeffs");
 468         word weightsFile;
 469
 470         if (metisCoeffs.readIfPresent("cellWeightsFile", weightsFile))
 471         {
 472             Info<< "parMetisDecomp : Using cell-based weights read from "
 473                 << weightsFile << endl;
 474
 475             labelIOField cellIOWeights
 476             (
 477                 IOobject
 478                 (
 479                     weightsFile,
 480                     mesh_.time().timeName(),
 481                     mesh_,
 482                     IOobject::MUST_READ,
 483                     IOobject::AUTO_WRITE
 484                 )
 485             );
 486             cellWeights.transfer(cellIOWeights);
 487
 488             if (cellWeights.size() != mesh_.nCells())
 489             {
 490                 FatalErrorIn
 491                 (
 492                     "parMetisDecomp::decompose"
 493                     "(const pointField&, const scalarField&)"
 494                 )   << "Number of cell weights " << cellWeights.size()
 495                     << " read from " << cellIOWeights.objectPath()
 496                     << " does not equal number of cells " << mesh_.nCells()
 497                     << exit(FatalError);
 498             }
 499         }
 500
 501         if (metisCoeffs.readIfPresent("faceWeightsFile", weightsFile))
 502         {
 503             Info<< "parMetisDecomp : Using face-based weights read from "
 504                 << weightsFile << endl;
 505
 506             labelIOField weights
 507             (
 508                 IOobject
 509                 (
 510                     weightsFile,
 511                     mesh_.time().timeName(),
 512                     mesh_,
 513                     IOobject::MUST_READ,
 514                     IOobject::AUTO_WRITE
 515                 )
 516             );
 517
 518             if (weights.size() != mesh_.nFaces())
 519             {
 520                 FatalErrorIn
 521                 (
 522                     "parMetisDecomp::decompose"
 523                     "(const pointField&, const scalarField&)"
 524                 )   << "Number of face weights " << weights.size()
 525                     << " does not equal number of internal and boundary faces "
 526                     << mesh_.nFaces()
 527                     << exit(FatalError);
 528             }
 529
 530             faceWeights.setSize(adjncy.size());
 531
 532             // Assume symmetric weights. Keep same ordering as adjncy.
 533             List<int> nFacesPerCell(mesh_.nCells(), 0);
 534
 535             // Handle internal faces
 536             for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 537             {
 538                 label w = weights[faceI];
 539
 540                 label own = mesh_.faceOwner()[faceI];
 541                 label nei = mesh_.faceNeighbour()[faceI];
 542
 543                 faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
 544                 faceWeights[xadj[nei] + nFacesPerCell[nei]++] = w;
 545             }
 546             // Coupled boundary faces
 547             const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 548
 549             forAll(patches, patchI)
 550             {
 551                 const polyPatch& pp = patches[patchI];
 552
 553                 if (pp.coupled())
 554                 {
 555                     label faceI = pp.start();
 556
 557                     forAll(pp, i)
 558                     {
 559                         label w = weights[faceI];
 560                         label own = mesh_.faceOwner()[faceI];
 561                         faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
 562                         faceI++;
 563                     }
 564                 }
 565             }
 566         }
 567
 568         if (metisCoeffs.readIfPresent("options", options))
 569         {
 570             Info<< "Using Metis options     " << options
 571                 << nl << endl;
 572
 573             if (options.size() != 3)
 574             {
 575                 FatalErrorIn
 576                 (
 577                     "parMetisDecomp::decompose"
 578                     "(const pointField&, const scalarField&)"
 579                 )   << "Number of options " << options.size()
 580                     << " should be three." << exit(FatalError);
 581             }
 582         }
 583     }
 584
 585
 586     // Do actual decomposition
 587     List<int> finalDecomp;
 588     decompose
 589     (
 590         xadj,
 591         adjncy,
 592         cc,
 593         cellWeights,
 594         faceWeights,
 595         options,
 596
 597         finalDecomp
 598     );
 599
 600     // Copy back to labelList
 601     labelList decomp(finalDecomp.size());
 602     forAll(decomp, i)
 603     {
 604         decomp[i] = finalDecomp[i];
 605     }
 606     return decomp;
 607 }
 608
 609
 610 Foam::labelList Foam::parMetisDecomp::decompose
 611 (
 612     const labelList& cellToRegion,
 613     const pointField& regionPoints,
 614     const scalarField& regionWeights
 615 )
 616 {
 617     const labelList& faceOwner = mesh_.faceOwner();
 618     const labelList& faceNeighbour = mesh_.faceNeighbour();
 619     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 620
 621     if (cellToRegion.size() != mesh_.nCells())
 622     {
 623         FatalErrorIn
 624         (
 625             "parMetisDecomp::decompose(const labelList&, const pointField&)"
 626         )   << "Size of cell-to-coarse map " << cellToRegion.size()
 627             << " differs from number of cells in mesh " << mesh_.nCells()
 628             << exit(FatalError);
 629     }
 630
 631
 632     // Global region numbering engine
 633     globalIndex globalRegions(regionPoints.size());
 634
 635
 636     // Get renumbered owner region on other side of coupled faces
 637     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 638
 639     List<int> globalNeighbour(mesh_.nFaces()-mesh_.nInternalFaces());
 640
 641     forAll(patches, patchI)
 642     {
 643         const polyPatch& pp = patches[patchI];
 644
 645         if (pp.coupled())
 646         {
 647             label faceI = pp.start();
 648             label bFaceI = pp.start() - mesh_.nInternalFaces();
 649
 650             forAll(pp, i)
 651             {
 652                 label ownRegion = cellToRegion[faceOwner[faceI]];
 653                 globalNeighbour[bFaceI++] = globalRegions.toGlobal(ownRegion);
 654                 faceI++;
 655             }
 656         }
 657     }
 658
 659     // Get the cell on the other side of coupled patches
 660     syncTools::swapBoundaryFaceList(mesh_, globalNeighbour, false);
 661
 662
 663     // Get globalCellCells on coarse mesh
 664     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 665
 666     labelListList globalRegionRegions;
 667     {
 668         List<DynamicList<label> > dynRegionRegions(regionPoints.size());
 669
 670         // Internal faces first
 671         forAll(faceNeighbour, faceI)
 672         {
 673             label ownRegion = cellToRegion[faceOwner[faceI]];
 674             label neiRegion = cellToRegion[faceNeighbour[faceI]];
 675
 676             if (ownRegion != neiRegion)
 677             {
 678                 label globalOwn = globalRegions.toGlobal(ownRegion);
 679                 label globalNei = globalRegions.toGlobal(neiRegion);
 680
 681                 if (findIndex(dynRegionRegions[ownRegion], globalNei) == -1)
 682                 {
 683                     dynRegionRegions[ownRegion].append(globalNei);
 684                 }
 685                 if (findIndex(dynRegionRegions[neiRegion], globalOwn) == -1)
 686                 {
 687                     dynRegionRegions[neiRegion].append(globalOwn);
 688                 }
 689             }
 690         }
 691
 692         // Coupled boundary faces
 693         forAll(patches, patchI)
 694         {
 695             const polyPatch& pp = patches[patchI];
 696
 697             if (pp.coupled())
 698             {
 699                 label faceI = pp.start();
 700                 label bFaceI = pp.start() - mesh_.nInternalFaces();
 701
 702                 forAll(pp, i)
 703                 {
 704                     label ownRegion = cellToRegion[faceOwner[faceI]];
 705                     label globalNei = globalNeighbour[bFaceI++];
 706                     faceI++;
 707
 708                     if (findIndex(dynRegionRegions[ownRegion], globalNei) == -1)
 709                     {
 710                         dynRegionRegions[ownRegion].append(globalNei);
 711                     }
 712                 }
 713             }
 714         }
 715
 716         globalRegionRegions.setSize(dynRegionRegions.size());
 717         forAll(dynRegionRegions, i)
 718         {
 719             globalRegionRegions[i].transfer(dynRegionRegions[i]);
 720         }
 721     }
 722
 723     labelList regionDecomp
 724     (
 725         decompose
 726         (
 727             globalRegionRegions,
 728             regionPoints,
 729             regionWeights
 730         )
 731     );
 732
 733     // Rework back into decomposition for original mesh_
 734     labelList cellDistribution(cellToRegion.size());
 735
 736     forAll(cellDistribution, cellI)
 737     {
 738         cellDistribution[cellI] = regionDecomp[cellToRegion[cellI]];
 739     }
 740     return cellDistribution;
 741 }
 742
 743
 744 Foam::labelList Foam::parMetisDecomp::decompose
 745 (
 746     const labelListList& globalCellCells,
 747     const pointField& cellCentres,
 748     const scalarField& cWeights
 749 )
 750 {
 751     if (cellCentres.size() != globalCellCells.size())
 752     {
 753         FatalErrorIn
 754         (
 755             "parMetisDecomp::decompose(const labelListList&"
 756             ", const pointField&, const scalarField&)"
 757         )   << "Inconsistent number of cells (" << globalCellCells.size()
 758             << ") and number of cell centres (" << cellCentres.size()
 759             << ") or weights (" << cWeights.size() << ")." << exit(FatalError);
 760     }
 761
 762     // For running sequential ...
 763     if (Pstream::nProcs() <= 1)
 764     {
 765         return metisDecomp(decompositionDict_, mesh_)
 766             .decompose(globalCellCells, cellCentres, cWeights);
 767     }
 768
 769
 770     // Make Metis Distributed CSR (Compressed Storage Format) storage
 771
 772     // Connections
 773     Field<int> adjncy;
 774     // Offsets into adjncy
 775     Field<int> xadj;
 776     scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);
 777
 778     // decomposition options. 0 = use defaults
 779     List<int> options(3, 0);
 780     //options[0] = 1;     // don't use defaults but use values below
 781     //options[1] = -1;    // full debug info
 782     //options[2] = 15;    // random number seed
 783
 784     // cell weights (so on the vertices of the dual)
 785     Field<int> cellWeights;
 786
 787     // face weights (so on the edges of the dual)
 788     Field<int> faceWeights;
 789
 790
 791     // Check for externally provided cellweights and if so initialise weights
 792     scalar minWeights = gMin(cWeights);
 793     if (cWeights.size() > 0)
 794     {
 795         if (minWeights <= 0)
 796         {
 797             WarningIn
 798             (
 799                 "parMetisDecomp::decompose(const labelListList&"
 800                 ", const pointField&, const scalarField&)"
 801             )   << "Illegal minimum weight " << minWeights
 802                 << endl;
 803         }
 804
 805         if (cWeights.size() != globalCellCells.size())
 806         {
 807             FatalErrorIn
 808             (
 809                 "parMetisDecomp::decompose(const labelListList&"
 810                 ", const pointField&, const scalarField&)"
 811             )   << "Number of cell weights " << cWeights.size()
 812                 << " does not equal number of cells " << globalCellCells.size()
 813                 << exit(FatalError);
 814         }
 815
 816         // Convert to integers.
 817         cellWeights.setSize(cWeights.size());
 818         forAll(cellWeights, i)
 819         {
 820             cellWeights[i] = int(cWeights[i]/minWeights);
 821         }
 822     }
 823
 824
 825     // Check for user supplied weights and decomp options
 826     if (decompositionDict_.found("metisCoeffs"))
 827     {
 828         const dictionary& metisCoeffs =
 829             decompositionDict_.subDict("metisCoeffs");
 830
 831         if (metisCoeffs.readIfPresent("options", options))
 832         {
 833             Info<< "Using Metis options     " << options
 834                 << nl << endl;
 835
 836             if (options.size() != 3)
 837             {
 838                 FatalErrorIn
 839                 (
 840                     "parMetisDecomp::decompose(const labelListList&"
 841                     ", const pointField&, const scalarField&)"
 842                 )   << "Number of options " << options.size()
 843                     << " should be three." << exit(FatalError);
 844             }
 845         }
 846     }
 847
 848
 849     // Do actual decomposition
 850     List<int> finalDecomp;
 851     decompose
 852     (
 853         xadj,
 854         adjncy,
 855         cellCentres,
 856         cellWeights,
 857         faceWeights,
 858         options,
 859
 860         finalDecomp
 861     );
 862
 863     // Copy back to labelList
 864     labelList decomp(finalDecomp.size());
 865     forAll(decomp, i)
 866     {
 867         decomp[i] = finalDecomp[i];
 868     }
 869     return decomp;
 870 }
 871
 872
 873 void Foam::parMetisDecomp::calcMetisDistributedCSR
 874 (
 875     const polyMesh& mesh,
 876     List<int>& adjncy,
 877     List<int>& xadj
 878 )
 879 {
 880     // Create global cell numbers
 881     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 882
 883     globalIndex globalCells(mesh.nCells());
 884
 885
 886     //
 887     // Make Metis Distributed CSR (Compressed Storage Format) storage
 888     //   adjncy      : contains cellCells (= edges in graph)
 889     //   xadj(celli) : start of information in adjncy for celli
 890     //
 891
 892
 893     const labelList& faceOwner = mesh.faceOwner();
 894     const labelList& faceNeighbour = mesh.faceNeighbour();
 895     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 896
 897
 898     // Get renumbered owner on other side of coupled faces
 899     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 900
 901     List<int> globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
 902
 903     forAll(patches, patchI)
 904     {
 905         const polyPatch& pp = patches[patchI];
 906
 907         if (pp.coupled())
 908         {
 909             label faceI = pp.start();
 910             label bFaceI = pp.start() - mesh.nInternalFaces();
 911
 912             forAll(pp, i)
 913             {
 914                 globalNeighbour[bFaceI++] = globalCells.toGlobal
 915                 (
 916                     faceOwner[faceI++]
 917                 );
 918             }
 919         }
 920     }
 921
 922     // Get the cell on the other side of coupled patches
 923     syncTools::swapBoundaryFaceList(mesh, globalNeighbour, false);
 924
 925
 926     // Count number of faces (internal + coupled)
 927     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 928
 929     // Number of faces per cell
 930     List<int> nFacesPerCell(mesh.nCells(), 0);
 931
 932     // Number of coupled faces
 933     label nCoupledFaces = 0;
 934
 935     for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
 936     {
 937         nFacesPerCell[faceOwner[faceI]]++;
 938         nFacesPerCell[faceNeighbour[faceI]]++;
 939     }
 940     // Handle coupled faces
 941     HashSet<edge, Hash<edge> > cellPair(mesh.nFaces()-mesh.nInternalFaces());
 942
 943     forAll(patches, patchI)
 944     {
 945         const polyPatch& pp = patches[patchI];
 946
 947         if (pp.coupled())
 948         {
 949             label faceI = pp.start();
 950
 951             forAll(pp, i)
 952             {
 953                 label own = faceOwner[faceI];
 954                 label globalNei = globalNeighbour[faceI-mesh.nInternalFaces()];
 955
 956                 if (cellPair.insert(edge(own, globalNei)))
 957                 {
 958                     nCoupledFaces++;
 959                     nFacesPerCell[faceOwner[faceI]]++;
 960                 }
 961                 faceI++;
 962             }
 963         }
 964     }
 965
 966
 967     // Fill in xadj
 968     // ~~~~~~~~~~~~
 969
 970     xadj.setSize(mesh.nCells()+1);
 971
 972     int freeAdj = 0;
 973
 974     for (label cellI = 0; cellI < mesh.nCells(); cellI++)
 975     {
 976         xadj[cellI] = freeAdj;
 977
 978         freeAdj += nFacesPerCell[cellI];
 979     }
 980     xadj[mesh.nCells()] = freeAdj;
 981
 982
 983
 984     // Fill in adjncy
 985     // ~~~~~~~~~~~~~~
 986
 987     adjncy.setSize(2*mesh.nInternalFaces() + nCoupledFaces);
 988
 989     nFacesPerCell = 0;
 990
 991     // For internal faces is just offsetted owner and neighbour
 992     for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
 993     {
 994         label own = faceOwner[faceI];
 995         label nei = faceNeighbour[faceI];
 996
 997         adjncy[xadj[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei);
 998         adjncy[xadj[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own);
 999     }
1000     // For boundary faces is offsetted coupled neighbour
1001     cellPair.clear();
1002     forAll(patches, patchI)
1003     {
1004         const polyPatch& pp = patches[patchI];
1005
1006         if (pp.coupled())
1007         {
1008             label faceI = pp.start();
1009             label bFaceI = pp.start()-mesh.nInternalFaces();
1010
1011             forAll(pp, i)
1012             {
1013                 label own = faceOwner[faceI];
1014                 label globalNei = globalNeighbour[bFaceI];
1015
1016                 if (cellPair.insert(edge(own, globalNei)))
1017                 {
1018                     adjncy[xadj[own] + nFacesPerCell[own]++] =
1019                         globalNeighbour[bFaceI];
1020                 }
1021
1022                 faceI++;
1023                 bFaceI++;
1024             }
1025         }
1026     }
1027 }
1028
1029
1030 // ************************************************************************* //