src/parallel/decompose/scotchDecomp/scotchDecomp.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     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     r
  69     {
  70         job=t,
  71         map=t,
  72         poli=S,
  73         sep=
  74         (
  75             m
  76             {
  77                 asc=b
  78                 {
  79                     bnd=
  80                     (
  81                         d{pass=40,dif=1,rem=1}
  82                      |
  83                     )
  84                     f{move=80,pass=-1,bal=0.002491},
  85                     org=f{move=80,pass=-1,bal=0.002491},
  86                     width=3
  87                 },
  88                 low=h{pass=10}
  89                 f{move=80,pass=-1,bal=0.002491},
  90                 type=h,
  91                 vert=80,
  92                 rat=0.8
  93             }
  94           | m
  95             {
  96                 asc=b
  97                 {
  98                     bnd=
  99                     (
 100                         d{pass=40,dif=1,rem=1}
 101                       |
 102                     )
 103                     f{move=80,pass=-1,bal=0.002491},
 104                     org=f{move=80,pass=-1,bal=0.002491},
 105                     width=3
 106                 },
 107                 low=h{pass=10}
 108                 f{move=80,pass=-1,bal=0.002491},
 109                 type=h,
 110                 vert=80,
 111                 rat=0.8
 112             }
 113         )
 114     }
 115
 116
 117     Note: instead of gmap run gpart <nProcs> -vs <grfFile>
 118     where <grfFile> can be obtained by running with 'writeGraph=true'
 119
 120 \*---------------------------------------------------------------------------*/
 121
 122 #include "scotchDecomp.H"
 123 #include "addToRunTimeSelectionTable.H"
 124 #include "floatScalar.H"
 125 #include "Time.H"
 126 #include "OFstream.H"
 127 #include "globalIndex.H"
 128 #include "SubField.H"
 129
 130 extern "C"
 131 {
 132 #include "scotch.h"
 133 }
 134
 135
 136 // Hack: scotch generates floating point errors so need to switch of error
 137 //       trapping!
 138 #if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
 139 #    define LINUX
 140 #endif
 141
 142 #if defined(LINUX) && defined(__GNUC__)
 143 #    define LINUX_GNUC
 144 #endif
 145
 146 #ifdef LINUX_GNUC
 147 #   ifndef __USE_GNU
 148 #       define __USE_GNU
 149 #   endif
 150 #   include <fenv.h>
 151 #endif
 152
 153
 154 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 155
 156 namespace Foam
 157 {
 158     defineTypeNameAndDebug(scotchDecomp, 0);
 159
 160     addToRunTimeSelectionTable
 161     (
 162         decompositionMethod,
 163         scotchDecomp,
 164         dictionary
 165     );
 166 }
 167
 168 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 169
 170 void Foam::scotchDecomp::check(const int retVal, const char* str)
 171 {
 172     if (retVal)
 173     {
 174         FatalErrorIn("scotchDecomp::decompose(..)")
 175             << "Call to scotch routine " << str << " failed."
 176             << exit(FatalError);
 177     }
 178 }
 179
 180
 181 Foam::label Foam::scotchDecomp::decompose
 182 (
 183     const fileName& meshPath,
 184     const List<int>& adjncy,
 185     const List<int>& xadj,
 186     const scalarField& cWeights,
 187
 188     List<int>& finalDecomp
 189 )
 190 {
 191     if (!Pstream::parRun())
 192     {
 193         decomposeOneProc
 194         (
 195             meshPath,
 196             adjncy,
 197             xadj,
 198             cWeights,
 199             finalDecomp
 200         );
 201     }
 202     else
 203     {
 204         if (debug)
 205         {
 206             Info<< "scotchDecomp : running in parallel."
 207                 << " Decomposing all of graph on master processor." << endl;
 208         }
 209         globalIndex globalCells(xadj.size()-1);
 210         label nTotalConnections = returnReduce(adjncy.size(), sumOp<label>());
 211
 212         // Send all to master. Use scheduled to save some storage.
 213         if (Pstream::master())
 214         {
 215             Field<int> allAdjncy(nTotalConnections);
 216             Field<int> allXadj(globalCells.size()+1);
 217             scalarField allWeights(globalCells.size());
 218
 219             // Insert my own
 220             label nTotalCells = 0;
 221             forAll(cWeights, cellI)
 222             {
 223                 allXadj[nTotalCells] = xadj[cellI];
 224                 allWeights[nTotalCells++] = cWeights[cellI];
 225             }
 226             nTotalConnections = 0;
 227             forAll(adjncy, i)
 228             {
 229                 allAdjncy[nTotalConnections++] = adjncy[i];
 230             }
 231
 232             for (int slave=1; slave<Pstream::nProcs(); slave++)
 233             {
 234                 IPstream fromSlave(Pstream::scheduled, slave);
 235                 Field<int> nbrAdjncy(fromSlave);
 236                 Field<int> nbrXadj(fromSlave);
 237                 scalarField nbrWeights(fromSlave);
 238
 239                 // Append.
 240                 //label procStart = nTotalCells;
 241                 forAll(nbrXadj, cellI)
 242                 {
 243                     allXadj[nTotalCells] = nTotalConnections+nbrXadj[cellI];
 244                     allWeights[nTotalCells++] = nbrWeights[cellI];
 245                 }
 246                 // No need to renumber xadj since already global.
 247                 forAll(nbrAdjncy, i)
 248                 {
 249                     allAdjncy[nTotalConnections++] = nbrAdjncy[i];
 250                 }
 251             }
 252             allXadj[nTotalCells] = nTotalConnections;
 253
 254
 255             Field<int> allFinalDecomp;
 256             decomposeOneProc
 257             (
 258                 meshPath,
 259                 allAdjncy,
 260                 allXadj,
 261                 allWeights,
 262                 allFinalDecomp
 263             );
 264
 265
 266             // Send allFinalDecomp back
 267             for (int slave=1; slave<Pstream::nProcs(); slave++)
 268             {
 269                 OPstream toSlave(Pstream::scheduled, slave);
 270                 toSlave << SubField<int>
 271                 (
 272                     allFinalDecomp,
 273                     globalCells.localSize(slave),
 274                     globalCells.offset(slave)
 275                 );
 276             }
 277             // Get my own part (always first)
 278             finalDecomp = SubField<int>
 279             (
 280                 allFinalDecomp,
 281                 globalCells.localSize()
 282             );
 283         }
 284         else
 285         {
 286             // Send my part of the graph (already in global numbering)
 287             {
 288                 OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
 289                 toMaster<< adjncy << SubField<int>(xadj, xadj.size()-1)
 290                     << cWeights;
 291             }
 292
 293             // Receive back decomposition
 294             IPstream fromMaster(Pstream::scheduled, Pstream::masterNo());
 295             fromMaster >> finalDecomp;
 296         }
 297     }
 298     return 0;
 299 }
 300
 301
 302 // Call scotch with options from dictionary.
 303 Foam::label Foam::scotchDecomp::decomposeOneProc
 304 (
 305     const fileName& meshPath,
 306     const List<int>& adjncy,
 307     const List<int>& xadj,
 308     const scalarField& cWeights,
 309
 310     List<int>& finalDecomp
 311 )
 312 {
 313     // Dump graph
 314     if (decompositionDict_.found("scotchCoeffs"))
 315     {
 316         const dictionary& scotchCoeffs =
 317             decompositionDict_.subDict("scotchCoeffs");
 318
 319         if (scotchCoeffs.lookupOrDefault("writeGraph", false))
 320         {
 321             OFstream str(meshPath + ".grf");
 322
 323             Info<< "Dumping Scotch graph file to " << str.name() << endl
 324                 << "Use this in combination with gpart." << endl;
 325
 326             label version = 0;
 327             str << version << nl;
 328             // Numer of vertices
 329             str << xadj.size()-1 << ' ' << adjncy.size() << nl;
 330             // Numbering starts from 0
 331             label baseval = 0;
 332             // Has weights?
 333             label hasEdgeWeights = 0;
 334             label hasVertexWeights = 0;
 335             label numericflag = 10*hasEdgeWeights+hasVertexWeights;
 336             str << baseval << ' ' << numericflag << nl;
 337             for (label cellI = 0; cellI < xadj.size()-1; cellI++)
 338             {
 339                 label start = xadj[cellI];
 340                 label end = xadj[cellI+1];
 341                 str << end-start;
 342
 343                 for (label i = start; i < end; i++)
 344                 {
 345                     str << ' ' << adjncy[i];
 346                 }
 347                 str << nl;
 348             }
 349         }
 350     }
 351
 352
 353     // Strategy
 354     // ~~~~~~~~
 355
 356     // Default.
 357     SCOTCH_Strat stradat;
 358     check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
 359
 360     if (decompositionDict_.found("scotchCoeffs"))
 361     {
 362         const dictionary& scotchCoeffs =
 363             decompositionDict_.subDict("scotchCoeffs");
 364
 365         string strategy;
 366         if (scotchCoeffs.readIfPresent("strategy", strategy))
 367         {
 368             if (debug)
 369             {
 370                 Info<< "scotchDecomp : Using strategy " << strategy << endl;
 371             }
 372             SCOTCH_stratGraphMap(&stradat, strategy.c_str());
 373             //fprintf(stdout, "S\tStrat=");
 374             //SCOTCH_stratSave(&stradat, stdout);
 375             //fprintf(stdout, "\n");
 376         }
 377     }
 378
 379
 380     // Graph
 381     // ~~~~~
 382
 383     List<int> velotab;
 384
 385
 386     // Check for externally provided cellweights and if so initialise weights
 387     // Note: min, not gMin since routine runs on master only.
 388     scalar minWeights = min(cWeights);
 389     if (cWeights.size() > 0)
 390     {
 391         if (minWeights <= 0)
 392         {
 393             WarningIn
 394             (
 395                 "scotchDecomp::decompose"
 396                 "(const pointField&, const scalarField&)"
 397             )   << "Illegal minimum weight " << minWeights
 398                 << endl;
 399         }
 400
 401         if (cWeights.size() != xadj.size()-1)
 402         {
 403             FatalErrorIn
 404             (
 405                 "scotchDecomp::decompose"
 406                 "(const pointField&, const scalarField&)"
 407             )   << "Number of cell weights " << cWeights.size()
 408                 << " does not equal number of cells " << xadj.size()-1
 409                 << exit(FatalError);
 410         }
 411
 412         scalar velotabSum = sum(cWeights)/minWeights;
 413
 414         scalar rangeScale(1.0);
 415
 416         if (velotabSum > scalar(INT_MAX - 1))
 417         {
 418             // 0.9 factor of safety to avoid floating point round-off in
 419             // rangeScale tipping the subsequent sum over the integer limit.
 420             rangeScale = 0.9*scalar(INT_MAX - 1)/velotabSum;
 421
 422             WarningIn
 423             (
 424                 "scotchDecomp::decompose"
 425                 "(const pointField&, const scalarField&)"
 426             )   << "Sum of weights has overflowed integer: " << velotabSum
 427                 << ", compressing weight scale by a factor of " << rangeScale
 428                 << endl;
 429         }
 430
 431         // Convert to integers.
 432         velotab.setSize(cWeights.size());
 433
 434         forAll(velotab, i)
 435         {
 436             velotab[i] =
 437                 int((cWeights[i]/minWeights - 1)*rangeScale) + 1;
 438         }
 439     }
 440
 441
 442
 443     SCOTCH_Graph grafdat;
 444     check(SCOTCH_graphInit(&grafdat), "SCOTCH_graphInit");
 445     check
 446     (
 447         SCOTCH_graphBuild
 448         (
 449             &grafdat,
 450             0,                      // baseval, c-style numbering
 451             xadj.size()-1,          // vertnbr, nCells
 452             xadj.begin(),           // verttab, start index per cell into adjncy
 453             &xadj[1],               // vendtab, end index  ,,
 454             velotab.begin(),        // velotab, vertex weights
 455             NULL,                   // vlbltab
 456             adjncy.size(),          // edgenbr, number of arcs
 457             adjncy.begin(),         // edgetab
 458             NULL                    // edlotab, edge weights
 459         ),
 460         "SCOTCH_graphBuild"
 461     );
 462     check(SCOTCH_graphCheck(&grafdat), "SCOTCH_graphCheck");
 463
 464
 465     // Architecture
 466     // ~~~~~~~~~~~~
 467     // (fully connected network topology since using switch)
 468
 469     SCOTCH_Arch archdat;
 470     check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
 471
 472     List<label> processorWeights;
 473     if (decompositionDict_.found("scotchCoeffs"))
 474     {
 475         const dictionary& scotchCoeffs =
 476             decompositionDict_.subDict("scotchCoeffs");
 477
 478         scotchCoeffs.readIfPresent("processorWeights", processorWeights);
 479     }
 480     if (processorWeights.size())
 481     {
 482         if (debug)
 483         {
 484             Info<< "scotchDecomp : Using procesor weights " << processorWeights
 485                 << endl;
 486         }
 487         check
 488         (
 489             SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
 490             "SCOTCH_archCmpltw"
 491         );
 492     }
 493     else
 494     {
 495         check
 496         (
 497             SCOTCH_archCmplt(&archdat, nProcessors_),
 498             "SCOTCH_archCmplt"
 499         );
 500     }
 501
 502
 503     //SCOTCH_Mapping mapdat;
 504     //SCOTCH_graphMapInit(&grafdat, &mapdat, &archdat, NULL);
 505     //SCOTCH_graphMapCompute(&grafdat, &mapdat, &stradat); /* Perform mapping */
 506     //SCOTCH_graphMapExit(&grafdat, &mapdat);
 507
 508
 509     // Hack:switch off fpu error trapping
 510 #   ifdef LINUX_GNUC
 511     int oldExcepts = fedisableexcept
 512     (
 513         FE_DIVBYZERO
 514       | FE_INVALID
 515       | FE_OVERFLOW
 516     );
 517 #   endif
 518
 519     finalDecomp.setSize(xadj.size()-1);
 520     finalDecomp = 0;
 521     check
 522     (
 523         SCOTCH_graphMap
 524         (
 525             &grafdat,
 526             &archdat,
 527             &stradat,           // const SCOTCH_Strat *
 528             finalDecomp.begin() // parttab
 529         ),
 530         "SCOTCH_graphMap"
 531     );
 532
 533 #   ifdef LINUX_GNUC
 534     feenableexcept(oldExcepts);
 535 #   endif
 536
 537
 538
 539     //finalDecomp.setSize(xadj.size()-1);
 540     //check
 541     //(
 542     //    SCOTCH_graphPart
 543     //    (
 544     //        &grafdat,
 545     //        nProcessors_,       // partnbr
 546     //        &stradat,           // const SCOTCH_Strat *
 547     //        finalDecomp.begin() // parttab
 548     //    ),
 549     //    "SCOTCH_graphPart"
 550     //);
 551
 552     // Release storage for graph
 553     SCOTCH_graphExit(&grafdat);
 554     // Release storage for strategy
 555     SCOTCH_stratExit(&stradat);
 556     // Release storage for network topology
 557     SCOTCH_archExit(&archdat);
 558
 559     return 0;
 560 }
 561
 562
 563 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 564
 565 Foam::scotchDecomp::scotchDecomp(const dictionary& decompositionDict)
 566 :
 567     decompositionMethod(decompositionDict)
 568 {}
 569
 570
 571 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 572
 573 Foam::labelList Foam::scotchDecomp::decompose
 574 (
 575     const polyMesh& mesh,
 576     const pointField& points,
 577     const scalarField& pointWeights
 578 )
 579 {
 580     if (points.size() != mesh.nCells())
 581     {
 582         FatalErrorIn
 583         (
 584             "scotchDecomp::decompose(const polyMesh&, const pointField&"
 585             ", const scalarField&)"
 586         )   << "Can use this decomposition method only for the whole mesh"
 587             << endl
 588             << "and supply one coordinate (cellCentre) for every cell." << endl
 589             << "The number of coordinates " << points.size() << endl
 590             << "The number of cells in the mesh " << mesh.nCells()
 591             << exit(FatalError);
 592     }
 593
 594     // Calculate local or global (if Pstream::parRun()) connectivity
 595     CompactListList<label> cellCells;
 596     calcCellCells(mesh, identity(mesh.nCells()), mesh.nCells(), cellCells);
 597
 598     // Decompose using default weights
 599     List<int> finalDecomp;
 600     decompose
 601     (
 602         mesh.time().path()/mesh.name(),
 603         cellCells.m(),
 604         cellCells.offsets(),
 605         pointWeights,
 606         finalDecomp
 607     );
 608
 609     // Copy back to labelList
 610     labelList decomp(finalDecomp.size());
 611     forAll(decomp, i)
 612     {
 613         decomp[i] = finalDecomp[i];
 614     }
 615     return decomp;
 616 }
 617
 618
 619 Foam::labelList Foam::scotchDecomp::decompose
 620 (
 621     const polyMesh& mesh,
 622     const labelList& agglom,
 623     const pointField& agglomPoints,
 624     const scalarField& pointWeights
 625 )
 626 {
 627     if (agglom.size() != mesh.nCells())
 628     {
 629         FatalErrorIn
 630         (
 631             "scotchDecomp::decompose"
 632             "(const polyMesh&, const labelList&, const pointField&"
 633             ", const scalarField&)"
 634         )   << "Size of cell-to-coarse map " << agglom.size()
 635             << " differs from number of cells in mesh " << mesh.nCells()
 636             << exit(FatalError);
 637     }
 638
 639     // Calculate local or global (if Pstream::parRun()) connectivity
 640     CompactListList<label> cellCells;
 641     calcCellCells(mesh, agglom, agglomPoints.size(), cellCells);
 642
 643     // Decompose using weights
 644     List<int> finalDecomp;
 645     decompose
 646     (
 647         mesh.time().path()/mesh.name(),
 648         cellCells.m(),
 649         cellCells.offsets(),
 650         pointWeights,
 651         finalDecomp
 652     );
 653
 654     // Rework back into decomposition for original mesh_
 655     labelList fineDistribution(agglom.size());
 656
 657     forAll(fineDistribution, i)
 658     {
 659         fineDistribution[i] = finalDecomp[agglom[i]];
 660     }
 661
 662     return fineDistribution;
 663 }
 664
 665
 666 Foam::labelList Foam::scotchDecomp::decompose
 667 (
 668     const labelListList& globalCellCells,
 669     const pointField& cellCentres,
 670     const scalarField& cWeights
 671 )
 672 {
 673     if (cellCentres.size() != globalCellCells.size())
 674     {
 675         FatalErrorIn
 676         (
 677             "scotchDecomp::decompose"
 678             "(const labelListList&, const pointField&, const scalarField&)"
 679         )   << "Inconsistent number of cells (" << globalCellCells.size()
 680             << ") and number of cell centres (" << cellCentres.size()
 681             << ")." << exit(FatalError);
 682     }
 683
 684
 685     // Make Metis CSR (Compressed Storage Format) storage
 686     //   adjncy      : contains neighbours (= edges in graph)
 687     //   xadj(celli) : start of information in adjncy for celli
 688
 689     CompactListList<label> cellCells(globalCellCells);
 690
 691     // Decompose using weights
 692     List<int> finalDecomp;
 693     decompose
 694     (
 695         "scotch",
 696         cellCells.m(),
 697         cellCells.offsets(),
 698         cWeights,
 699         finalDecomp
 700     );
 701
 702     // Copy back to labelList
 703     labelList decomp(finalDecomp.size());
 704     forAll(decomp, i)
 705     {
 706         decomp[i] = finalDecomp[i];
 707     }
 708     return decomp;
 709 }
 710
 711
 712 // ************************************************************************* //