src/decompositionMethods/scotchDecomp/scotchDecomp.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     From scotch forum:
  25
  26     By: Francois PELLEGRINI RE: Graph mapping 'strategy' string [ reply ]
  27     2008-08-22 10:09 Strategy handling in Scotch is a bit tricky. In order
  28     not to be confused, you must have a clear view of how they are built.
  29     Here are some rules:
  30
  31     1- Strategies are made up of "methods" which are combined by means of
  32     "operators".
  33
  34     2- A method is of the form "m{param=value,param=value,...}", where "m"
  35     is a single character (this is your first error: "f" is a method name,
  36     not a parameter name).
  37
  38     3- There exist different sort of strategies : bipartitioning strategies,
  39     mapping strategies, ordering strategies, which cannot be mixed. For
  40     instance, you cannot build a bipartitioning strategy and feed it to a
  41     mapping method (this is your second error).
  42
  43     To use the "mapCompute" routine, you must create a mapping strategy, not
  44     a bipartitioning one, and so use stratGraphMap() and not
  45     stratGraphBipart(). Your mapping strategy should however be based on the
  46     "recursive bipartitioning" method ("b"). For instance, a simple (and
  47     hence not very efficient) mapping strategy can be :
  48
  49     "b{sep=f}"
  50
  51     which computes mappings with the recursive bipartitioning method "b",
  52     this latter using the Fiduccia-Mattheyses method "f" to compute its
  53     separators.
  54
  55     If you want an exact partition (see your previous post), try
  56     "b{sep=fx}".
  57
  58     However, these strategies are not the most efficient, as they do not
  59     make use of the multi-level framework.
  60
  61     To use the multi-level framework, try for instance:
  62
  63     "b{sep=m{vert=100,low=h,asc=f}x}"
  64
  65     The current default mapping strategy in Scotch can be seen by using the
  66     "-vs" option of program gmap. It is, to date:
  67
  68     b
  69     {
  70         job=t,
  71         map=t,
  72         poli=S,
  73         sep=
  74         (
  75             m
  76             {
  77                 asc=b
  78                 {
  79                     bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
  80                     org=f{move=80,pass=-1,bal=0.005},
  81                     width=3
  82                 },
  83                 low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
  84                 type=h,
  85                 vert=80,
  86                 rat=0.8
  87             }
  88           | m
  89             {
  90                 asc=b
  91                 {
  92                     bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
  93                     org=f{move=80,pass=-1,bal=0.005},
  94                     width=3
  95                 },
  96                 low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
  97                 type=h,
  98                 vert=80,
  99                 rat=0.8
 100             }
 101         )
 102     }
 103
 104
 105 \*---------------------------------------------------------------------------*/
 106
 107 #include "scotchDecomp.H"
 108 #include "addToRunTimeSelectionTable.H"
 109 #include "floatScalar.H"
 110 #include "Time.H"
 111 #include "cyclicPolyPatch.H"
 112 #include "OFstream.H"
 113
 114 extern "C"
 115 {
 116 #include "scotch.h"
 117 }
 118
 119
 120 // Hack: scotch generates floating point errors so need to switch of error
 121 //       trapping!
 122 #if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
 123 #    define LINUX
 124 #endif
 125
 126 #if defined(LINUX) && defined(__GNUC__)
 127 #    define LINUX_GNUC
 128 #endif
 129
 130 #ifdef LINUX_GNUC
 131 #   ifndef __USE_GNU
 132 #       define __USE_GNU
 133 #   endif
 134 #   include <fenv.h>
 135 #endif
 136
 137
 138 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 139
 140 namespace Foam
 141 {
 142     defineTypeNameAndDebug(scotchDecomp, 0);
 143
 144     addToRunTimeSelectionTable
 145     (
 146         decompositionMethod,
 147         scotchDecomp,
 148         dictionaryMesh
 149     );
 150 }
 151
 152 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 153
 154 void Foam::scotchDecomp::check(const int retVal, const char* str)
 155 {
 156     if (retVal)
 157     {
 158         FatalErrorIn("scotchDecomp::decompose(..)")
 159             << "Call to scotch routine " << str << " failed."
 160             << exit(FatalError);
 161     }
 162 }
 163
 164
 165 // Call scotch with options from dictionary.
 166 Foam::label Foam::scotchDecomp::decompose
 167 (
 168     const List<int>& adjncy,
 169     const List<int>& xadj,
 170     const scalarField& cWeights,
 171
 172     List<int>& finalDecomp
 173 )
 174 {
 175     // Dump graph
 176     if (decompositionDict_.found("scotchCoeffs"))
 177     {
 178         const dictionary& scotchCoeffs =
 179             decompositionDict_.subDict("scotchCoeffs");
 180
 181         if (scotchCoeffs.found("writeGraph"))
 182         {
 183             Switch writeGraph(scotchCoeffs.lookup("writeGraph"));
 184
 185             if (writeGraph)
 186             {
 187                 OFstream str(mesh_.time().path() / mesh_.name() + ".grf");
 188
 189                 Info<< "Dumping Scotch graph file to " << str.name() << endl
 190                     << "Use this in combination with gpart." << endl;
 191
 192                 label version = 0;
 193                 str << version << nl;
 194                 // Numer of vertices
 195                 str << xadj.size()-1 << ' ' << adjncy.size() << nl;
 196                 // Numbering starts from 0
 197                 label baseval = 0;
 198                 // Has weights?
 199                 label hasEdgeWeights = 0;
 200                 label hasVertexWeights = 0;
 201                 label numericflag = 10*hasEdgeWeights+hasVertexWeights;
 202                 str << baseval << ' ' << numericflag << nl;
 203                 for (label cellI = 0; cellI < xadj.size()-1; cellI++)
 204                 {
 205                     label start = xadj[cellI];
 206                     label end = xadj[cellI+1];
 207                     str << end-start;
 208
 209                     for (label i = start; i < end; i++)
 210                     {
 211                         str << ' ' << adjncy[i];
 212                     }
 213                     str << nl;
 214                 }
 215             }
 216         }
 217     }
 218
 219
 220     // Strategy
 221     // ~~~~~~~~
 222
 223     // Default.
 224     SCOTCH_Strat stradat;
 225     check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
 226
 227     if (decompositionDict_.found("scotchCoeffs"))
 228     {
 229         const dictionary& scotchCoeffs =
 230             decompositionDict_.subDict("scotchCoeffs");
 231
 232
 233         string strategy;
 234         if (scotchCoeffs.readIfPresent("strategy", strategy))
 235         {
 236             if (debug)
 237             {
 238                 Info<< "scotchDecomp : Using strategy " << strategy << endl;
 239             }
 240             SCOTCH_stratGraphMap(&stradat, strategy.c_str());
 241             //fprintf(stdout, "S\tStrat=");
 242             //SCOTCH_stratSave(&stradat, stdout);
 243             //fprintf(stdout, "\n");
 244         }
 245     }
 246
 247
 248     // Graph
 249     // ~~~~~
 250
 251     List<int> velotab;
 252
 253
 254     // Check for externally provided cellweights and if so initialise weights
 255     scalar minWeights = gMin(cWeights);
 256     if (cWeights.size() > 0)
 257     {
 258         if (minWeights <= 0)
 259         {
 260             WarningIn
 261             (
 262                 "scotchDecomp::decompose"
 263                 "(const pointField&, const scalarField&)"
 264             )   << "Illegal minimum weight " << minWeights
 265                 << endl;
 266         }
 267
 268         if (cWeights.size() != xadj.size()-1)
 269         {
 270             FatalErrorIn
 271             (
 272                 "scotchDecomp::decompose"
 273                 "(const pointField&, const scalarField&)"
 274             )   << "Number of cell weights " << cWeights.size()
 275                 << " does not equal number of cells " << xadj.size()-1
 276                 << exit(FatalError);
 277         }
 278
 279         // Convert to integers.
 280         velotab.setSize(cWeights.size());
 281         forAll(velotab, i)
 282         {
 283             velotab[i] = int(cWeights[i]/minWeights);
 284         }
 285     }
 286
 287
 288
 289     SCOTCH_Graph grafdat;
 290     check(SCOTCH_graphInit(&grafdat), "SCOTCH_graphInit");
 291     check
 292     (
 293         SCOTCH_graphBuild
 294         (
 295             &grafdat,
 296             0,                      // baseval, c-style numbering
 297             xadj.size()-1,          // vertnbr, nCells
 298             xadj.begin(),           // verttab, start index per cell into adjncy
 299             &xadj[1],               // vendtab, end index  ,,
 300             velotab.begin(),        // velotab, vertex weights
 301             NULL,                   // vlbltab
 302             adjncy.size(),          // edgenbr, number of arcs
 303             adjncy.begin(),         // edgetab
 304             NULL                    // edlotab, edge weights
 305         ),
 306         "SCOTCH_graphBuild"
 307     );
 308     check(SCOTCH_graphCheck(&grafdat), "SCOTCH_graphCheck");
 309
 310
 311     // Architecture
 312     // ~~~~~~~~~~~~
 313     // (fully connected network topology since using switch)
 314
 315     SCOTCH_Arch archdat;
 316     check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
 317
 318     List<label> processorWeights;
 319     if (decompositionDict_.found("scotchCoeffs"))
 320     {
 321         const dictionary& scotchCoeffs =
 322             decompositionDict_.subDict("scotchCoeffs");
 323
 324         scotchCoeffs.readIfPresent("processorWeights", processorWeights);
 325     }
 326     if (processorWeights.size())
 327     {
 328         if (debug)
 329         {
 330             Info<< "scotchDecomp : Using procesor weights " << processorWeights
 331                 << endl;
 332         }
 333         check
 334         (
 335             SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
 336             "SCOTCH_archCmpltw"
 337         );
 338     }
 339     else
 340     {
 341         check
 342         (
 343             SCOTCH_archCmplt(&archdat, nProcessors_),
 344             "SCOTCH_archCmplt"
 345         );
 346     }
 347
 348
 349     //SCOTCH_Mapping mapdat;
 350     //SCOTCH_graphMapInit(&grafdat, &mapdat, &archdat, NULL);
 351     //SCOTCH_graphMapCompute(&grafdat, &mapdat, &stradat); /* Perform mapping */
 352     //SCOTCH_graphMapExit(&grafdat, &mapdat);
 353
 354
 355     // Hack:switch off fpu error trapping
 356 #   ifdef LINUX_GNUC
 357     int oldExcepts = fedisableexcept
 358     (
 359         FE_DIVBYZERO
 360       | FE_INVALID
 361       | FE_OVERFLOW
 362     );
 363 #   endif
 364
 365     finalDecomp.setSize(xadj.size()-1);
 366     finalDecomp = 0;
 367     check
 368     (
 369         SCOTCH_graphMap
 370         (
 371             &grafdat,
 372             &archdat,
 373             &stradat,           // const SCOTCH_Strat *
 374             finalDecomp.begin() // parttab
 375         ),
 376         "SCOTCH_graphMap"
 377     );
 378
 379 #   ifdef LINUX_GNUC
 380     feenableexcept(oldExcepts);
 381 #   endif
 382
 383
 384
 385     //finalDecomp.setSize(xadj.size()-1);
 386     //check
 387     //(
 388     //    SCOTCH_graphPart
 389     //    (
 390     //        &grafdat,
 391     //        nProcessors_,       // partnbr
 392     //        &stradat,           // const SCOTCH_Strat *
 393     //        finalDecomp.begin() // parttab
 394     //    ),
 395     //    "SCOTCH_graphPart"
 396     //);
 397
 398     // Release storage for graph
 399     SCOTCH_graphExit(&grafdat);
 400     // Release storage for strategy
 401     SCOTCH_stratExit(&stradat);
 402     // Release storage for network topology
 403     SCOTCH_archExit(&archdat);
 404
 405     return 0;
 406 }
 407
 408
 409 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 410
 411 Foam::scotchDecomp::scotchDecomp
 412 (
 413     const dictionary& decompositionDict,
 414     const polyMesh& mesh
 415 )
 416 :
 417     decompositionMethod(decompositionDict),
 418     mesh_(mesh)
 419 {}
 420
 421
 422 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 423
 424 Foam::labelList Foam::scotchDecomp::decompose
 425 (
 426     const pointField& points,
 427     const scalarField& pointWeights
 428 )
 429 {
 430     if (points.size() != mesh_.nCells())
 431     {
 432         FatalErrorIn
 433         (
 434             "scotchDecomp::decompose(const pointField&, const scalarField&)"
 435         )
 436             << "Can use this decomposition method only for the whole mesh"
 437             << endl
 438             << "and supply one coordinate (cellCentre) for every cell." << endl
 439             << "The number of coordinates " << points.size() << endl
 440             << "The number of cells in the mesh " << mesh_.nCells()
 441             << exit(FatalError);
 442     }
 443
 444     // Make Metis CSR (Compressed Storage Format) storage
 445     //   adjncy      : contains neighbours (= edges in graph)
 446     //   xadj(celli) : start of information in adjncy for celli
 447     List<int> adjncy;
 448     List<int> xadj;
 449     calcCSR(mesh_, adjncy, xadj);
 450
 451     // Decompose using default weights
 452     List<int> finalDecomp;
 453     decompose(adjncy, xadj, pointWeights, finalDecomp);
 454
 455     // Copy back to labelList
 456     labelList decomp(finalDecomp.size());
 457     forAll(decomp, i)
 458     {
 459         decomp[i] = finalDecomp[i];
 460     }
 461     return decomp;
 462 }
 463
 464
 465 Foam::labelList Foam::scotchDecomp::decompose
 466 (
 467     const labelList& agglom,
 468     const pointField& agglomPoints,
 469     const scalarField& pointWeights
 470 )
 471 {
 472     if (agglom.size() != mesh_.nCells())
 473     {
 474         FatalErrorIn
 475         (
 476             "parMetisDecomp::decompose(const labelList&, const pointField&)"
 477         )   << "Size of cell-to-coarse map " << agglom.size()
 478             << " differs from number of cells in mesh " << mesh_.nCells()
 479             << exit(FatalError);
 480     }
 481
 482     // Make Metis CSR (Compressed Storage Format) storage
 483     //   adjncy      : contains neighbours (= edges in graph)
 484     //   xadj(celli) : start of information in adjncy for celli
 485     List<int> adjncy;
 486     List<int> xadj;
 487     {
 488         // Get cellCells on coarse mesh.
 489         labelListList cellCells;
 490         calcCellCells
 491         (
 492             mesh_,
 493             agglom,
 494             agglomPoints.size(),
 495             cellCells
 496         );
 497
 498         calcCSR(cellCells, adjncy, xadj);
 499     }
 500
 501     // Decompose using weights
 502     List<int> finalDecomp;
 503     decompose(adjncy, xadj, pointWeights, finalDecomp);
 504
 505     // Rework back into decomposition for original mesh_
 506     labelList fineDistribution(agglom.size());
 507
 508     forAll(fineDistribution, i)
 509     {
 510         fineDistribution[i] = finalDecomp[agglom[i]];
 511     }
 512
 513     return fineDistribution;
 514 }
 515
 516
 517 Foam::labelList Foam::scotchDecomp::decompose
 518 (
 519     const labelListList& globalCellCells,
 520     const pointField& cellCentres,
 521     const scalarField& cWeights
 522 )
 523 {
 524     if (cellCentres.size() != globalCellCells.size())
 525     {
 526         FatalErrorIn
 527         (
 528             "scotchDecomp::decompose"
 529             "(const labelListList&, const pointField&, const scalarField&)"
 530         )   << "Inconsistent number of cells (" << globalCellCells.size()
 531             << ") and number of cell centres (" << cellCentres.size()
 532             << ")." << exit(FatalError);
 533     }
 534
 535
 536     // Make Metis CSR (Compressed Storage Format) storage
 537     //   adjncy      : contains neighbours (= edges in graph)
 538     //   xadj(celli) : start of information in adjncy for celli
 539
 540     List<int> adjncy;
 541     List<int> xadj;
 542     calcCSR(globalCellCells, adjncy, xadj);
 543
 544     // Decompose using weights
 545     List<int> finalDecomp;
 546     decompose(adjncy, xadj, cWeights, finalDecomp);
 547
 548     // Copy back to labelList
 549     labelList decomp(finalDecomp.size());
 550     forAll(decomp, i)
 551     {
 552         decomp[i] = finalDecomp[i];
 553     }
 554     return decomp;
 555 }
 556
 557
 558 void Foam::scotchDecomp::calcCSR
 559 (
 560     const polyMesh& mesh,
 561     List<int>& adjncy,
 562     List<int>& xadj
 563 )
 564 {
 565     // Make Metis CSR (Compressed Storage Format) storage
 566     //   adjncy      : contains neighbours (= edges in graph)
 567     //   xadj(celli) : start of information in adjncy for celli
 568
 569     xadj.setSize(mesh.nCells()+1);
 570
 571     // Initialise the number of internal faces of the cells to twice the
 572     // number of internal faces
 573     label nInternalFaces = 2*mesh.nInternalFaces();
 574
 575     // Check the boundary for coupled patches and add to the number of
 576     // internal faces
 577     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 578
 579     forAll(pbm, patchi)
 580     {
 581         if (isA<cyclicPolyPatch>(pbm[patchi]))
 582         {
 583             nInternalFaces += pbm[patchi].size();
 584         }
 585     }
 586
 587     // Create the adjncy array the size of the total number of internal and
 588     // coupled faces
 589     adjncy.setSize(nInternalFaces);
 590
 591     // Fill in xadj
 592     // ~~~~~~~~~~~~
 593     label freeAdj = 0;
 594
 595     for (label cellI = 0; cellI < mesh.nCells(); cellI++)
 596     {
 597         xadj[cellI] = freeAdj;
 598
 599         const labelList& cFaces = mesh.cells()[cellI];
 600
 601         forAll(cFaces, i)
 602         {
 603             label faceI = cFaces[i];
 604
 605             if
 606             (
 607                 mesh.isInternalFace(faceI)
 608              || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
 609             )
 610             {
 611                 freeAdj++;
 612             }
 613         }
 614     }
 615     xadj[mesh.nCells()] = freeAdj;
 616
 617
 618     // Fill in adjncy
 619     // ~~~~~~~~~~~~~~
 620
 621     labelList nFacesPerCell(mesh.nCells(), 0);
 622
 623     // Internal faces
 624     for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
 625     {
 626         label own = mesh.faceOwner()[faceI];
 627         label nei = mesh.faceNeighbour()[faceI];
 628
 629         adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
 630         adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
 631     }
 632
 633     // Coupled faces. Only cyclics done.
 634     forAll(pbm, patchi)
 635     {
 636         if (isA<cyclicPolyPatch>(pbm[patchi]))
 637         {
 638             const unallocLabelList& faceCells = pbm[patchi].faceCells();
 639
 640             label sizeby2 = faceCells.size()/2;
 641
 642             for (label facei=0; facei<sizeby2; facei++)
 643             {
 644                 label own = faceCells[facei];
 645                 label nei = faceCells[facei + sizeby2];
 646
 647                 adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
 648                 adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
 649             }
 650         }
 651     }
 652 }
 653
 654
 655 // From cell-cell connections to Metis format (like CompactListList)
 656 void Foam::scotchDecomp::calcCSR
 657 (
 658     const labelListList& cellCells,
 659     List<int>& adjncy,
 660     List<int>& xadj
 661 )
 662 {
 663     // Count number of internal faces
 664     label nConnections = 0;
 665
 666     forAll(cellCells, coarseI)
 667     {
 668         nConnections += cellCells[coarseI].size();
 669     }
 670
 671     // Create the adjncy array as twice the size of the total number of
 672     // internal faces
 673     adjncy.setSize(nConnections);
 674
 675     xadj.setSize(cellCells.size()+1);
 676
 677
 678     // Fill in xadj
 679     // ~~~~~~~~~~~~
 680     label freeAdj = 0;
 681
 682     forAll(cellCells, coarseI)
 683     {
 684         xadj[coarseI] = freeAdj;
 685
 686         const labelList& cCells = cellCells[coarseI];
 687
 688         forAll(cCells, i)
 689         {
 690             adjncy[freeAdj++] = cCells[i];
 691         }
 692     }
 693     xadj[cellCells.size()] = freeAdj;
 694 }
 695
 696
 697
 698 // ************************************************************************* //