src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.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     Default strategy:
  26         Strat=b
  27         {
  28             job=t,
  29             map=t,
  30             poli=S,
  31             sep=
  32             (
  33                 m
  34                 {
  35                     asc=b
  36                     {
  37                         bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
  38                         org=f{move=80,pass=-1,bal=0.005},width=3
  39                     },
  40                     low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
  41                     type=h,
  42                     vert=80,
  43                     rat=0.8
  44                 }
  45               | m
  46                 {
  47                     asc=b
  48                     {
  49                         bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
  50                         org=f{move=80,pass=-1,bal=0.005},
  51                         width=3},
  52                         low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
  53                         type=h,
  54                         vert=80,
  55                         rat=0.8
  56                 }
  57             )
  58         }
  59 \*---------------------------------------------------------------------------*/
  60
  61 #include "scotchDecomp.H"
  62 #include "addToRunTimeSelectionTable.H"
  63 #include "floatScalar.H"
  64 #include "IFstream.H"
  65 #include "Time.H"
  66 #include "cyclicPolyPatch.H"
  67 #include "OFstream.H"
  68
  69 extern "C"
  70 {
  71 #define OMPI_SKIP_MPICXX
  72 #include "module.h"
  73 #include "common.h"
  74 #include "scotch.h"
  75 }
  76
  77
  78 // Hack: scotch generates floating point errors so need to switch of error
  79 //       trapping!
  80 #if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
  81 #    define LINUX
  82 #endif
  83
  84 #if defined(LINUX) && defined(__GNUC__)
  85 #    define LINUX_GNUC
  86 #endif
  87
  88 #ifdef LINUX_GNUC
  89 #   ifndef __USE_GNU
  90 #       define __USE_GNU
  91 #   endif
  92 #   include <fenv.h>
  93 #endif
  94
  95
  96 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  97
  98 namespace Foam
  99 {
 100     defineTypeNameAndDebug(scotchDecomp, 0);
 101
 102     addToRunTimeSelectionTable
 103     (
 104         decompositionMethod,
 105         scotchDecomp,
 106         dictionaryMesh
 107     );
 108 }
 109
 110 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 111
 112 void Foam::scotchDecomp::check(const int retVal, const char* str)
 113 {
 114     if (retVal)
 115     {
 116         FatalErrorIn("scotchDecomp::decompose(..)")
 117             << "Called to scotch routine " << str << " failed."
 118             << exit(FatalError);
 119     }
 120 }
 121
 122
 123 // Call Metis with options from dictionary.
 124 Foam::label Foam::scotchDecomp::decompose
 125 (
 126     const List<int>& adjncy,
 127     const List<int>& xadj,
 128
 129     List<int>& finalDecomp
 130 )
 131 {
 132     // Strategy
 133     // ~~~~~~~~
 134
 135     // Default.
 136     SCOTCH_Strat stradat;
 137     check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
 138
 139     if (decompositionDict_.found("scotchCoeffs"))
 140     {
 141         const dictionary& scotchCoeffs =
 142             decompositionDict_.subDict("scotchCoeffs");
 143
 144
 145         string strategy;
 146         if (scotchCoeffs.readIfPresent("strategy", strategy))
 147         {
 148             if (debug)
 149             {
 150                 Info<< "scotchDecomp : Using strategy " << strategy << endl;
 151             }
 152             SCOTCH_stratGraphMap(&stradat, strategy.c_str());
 153             //fprintf(stdout, "S\tStrat=");
 154             //SCOTCH_stratSave(&stradat, stdout);
 155             //fprintf(stdout, "\n");
 156         }
 157     }
 158
 159
 160     // Graph
 161     // ~~~~~
 162
 163     SCOTCH_Graph grafdat;
 164     check(SCOTCH_graphInit(&grafdat), "SCOTCH_graphInit");
 165     check
 166     (
 167         SCOTCH_graphBuild
 168         (
 169             &grafdat,
 170             0,                      // baseval, c-style numbering
 171             xadj.size()-1,          // vertnbr, nCells
 172             xadj.begin(),           // verttab, start index per cell into adjncy
 173             &xadj[1],               // vendtab, end index  ,,
 174             NULL,                   // velotab, vertex weights
 175             NULL,                   // vlbltab
 176             adjncy.size(),          // edgenbr, number of arcs
 177             adjncy.begin(),         // edgetab
 178             NULL                    // edlotab, edge weights
 179         ),
 180         "SCOTCH_graphBuild"
 181     );
 182     check(SCOTCH_graphCheck(&grafdat), "SCOTCH_graphCheck");
 183
 184
 185     // Dump graph
 186     if (decompositionDict_.found("scotchCoeffs"))
 187     {
 188         const dictionary& scotchCoeffs =
 189             decompositionDict_.subDict("scotchCoeffs");
 190
 191         if (scotchCoeffs.found("writeGraph"))
 192         {
 193             Switch writeGraph(scotchCoeffs.lookup("writeGraph"));
 194
 195             if (writeGraph)
 196             {
 197                 OFstream str(mesh_.time().path() / mesh_.name() + ".grf");
 198
 199                 Info<< "Dumping Scotch graph file to " << str.name() << endl
 200                     << "Use this in combination with gpart." << endl;
 201
 202                 label version = 0;
 203                 str << version << nl;
 204                 // Numer of vertices
 205                 str << xadj.size()-1 << ' ' << adjncy.size() << nl;
 206                 // Numbering starts from 0
 207                 label baseval = 0;
 208                 // Has weights?
 209                 label hasEdgeWeights = 0;
 210                 label hasVertexWeights = 0;
 211                 label numericflag = 10*hasEdgeWeights+hasVertexWeights;
 212                 str << baseval << ' ' << numericflag << nl;
 213                 for (label cellI = 0; cellI < xadj.size()-1; cellI++)
 214                 {
 215                     label start = xadj[cellI];
 216                     label end = xadj[cellI+1];
 217                     str << end-start;
 218
 219                     for (label i = start; i < end; i++)
 220                     {
 221                         str << ' ' << adjncy[i];
 222                     }
 223                     str << nl;
 224                 }
 225             }
 226         }
 227     }
 228
 229
 230     // Architecture
 231     // ~~~~~~~~~~~~
 232     // (fully connected network topology since using switch)
 233
 234     SCOTCH_Arch archdat;
 235     check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
 236
 237     List<label> processorWeights;
 238     if (decompositionDict_.found("scotchCoeffs"))
 239     {
 240         const dictionary& scotchCoeffs =
 241             decompositionDict_.subDict("scotchCoeffs");
 242
 243         scotchCoeffs.readIfPresent("processorWeights", processorWeights);
 244     }
 245     if (processorWeights.size())
 246     {
 247         if (debug)
 248         {
 249             Info<< "scotchDecomp : Using procesor weights " << processorWeights
 250                 << endl;
 251         }
 252         check
 253         (
 254             SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
 255             "SCOTCH_archCmpltw"
 256         );
 257     }
 258     else
 259     {
 260         check
 261         (
 262             SCOTCH_archCmplt(&archdat, nProcessors_),
 263             "SCOTCH_archCmplt"
 264         );
 265     }
 266
 267
 268     //SCOTCH_Mapping mapdat;
 269     //SCOTCH_graphMapInit(&grafdat, &mapdat, &archdat, NULL);
 270     //SCOTCH_graphMapCompute(&grafdat, &mapdat, &stradat); /* Perform mapping */
 271     //SCOTCH_graphMapExit(&grafdat, &mapdat);
 272
 273
 274     // Hack:switch off fpu error trapping
 275 #   ifdef LINUX_GNUC
 276     int oldExcepts = fedisableexcept
 277     (
 278         FE_DIVBYZERO
 279       | FE_INVALID
 280       | FE_OVERFLOW
 281     );
 282 #   endif
 283
 284     finalDecomp.setSize(xadj.size()-1);
 285     finalDecomp = 0;
 286     check
 287     (
 288         SCOTCH_graphMap
 289         (
 290             &grafdat,
 291             &archdat,
 292             &stradat,           // const SCOTCH_Strat *
 293             finalDecomp.begin() // parttab
 294         ),
 295         "SCOTCH_graphMap"
 296     );
 297
 298 #   ifdef LINUX_GNUC
 299     feenableexcept(oldExcepts);
 300 #   endif
 301
 302
 303
 304     //finalDecomp.setSize(xadj.size()-1);
 305     //check
 306     //(
 307     //    SCOTCH_graphPart
 308     //    (
 309     //        &grafdat,
 310     //        nProcessors_,       // partnbr
 311     //        &stradat,           // const SCOTCH_Strat *
 312     //        finalDecomp.begin() // parttab
 313     //    ),
 314     //    "SCOTCH_graphPart"
 315     //);
 316
 317     // Release storage for graph
 318     SCOTCH_graphExit(&grafdat);
 319     // Release storage for strategy
 320     SCOTCH_stratExit(&stradat);
 321     // Release storage for network topology
 322     SCOTCH_archExit(&archdat);
 323
 324     return 0;
 325 }
 326
 327
 328 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 329
 330 Foam::scotchDecomp::scotchDecomp
 331 (
 332     const dictionary& decompositionDict,
 333     const polyMesh& mesh
 334 )
 335 :
 336     decompositionMethod(decompositionDict),
 337     mesh_(mesh)
 338 {}
 339
 340
 341 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 342
 343 Foam::labelList Foam::scotchDecomp::decompose(const pointField& points)
 344 {
 345     if (points.size() != mesh_.nCells())
 346     {
 347         FatalErrorIn("scotchDecomp::decompose(const pointField&)")
 348             << "Can use this decomposition method only for the whole mesh"
 349             << endl
 350             << "and supply one coordinate (cellCentre) for every cell." << endl
 351             << "The number of coordinates " << points.size() << endl
 352             << "The number of cells in the mesh " << mesh_.nCells()
 353             << exit(FatalError);
 354     }
 355
 356     // Make Metis CSR (Compressed Storage Format) storage
 357     //   adjncy      : contains neighbours (= edges in graph)
 358     //   xadj(celli) : start of information in adjncy for celli
 359
 360     List<int> xadj(mesh_.nCells()+1);
 361
 362     // Initialise the number of internal faces of the cells to twice the
 363     // number of internal faces
 364     label nInternalFaces = 2*mesh_.nInternalFaces();
 365
 366     // Check the boundary for coupled patches and add to the number of
 367     // internal faces
 368     const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
 369
 370     forAll(pbm, patchi)
 371     {
 372         if (isA<cyclicPolyPatch>(pbm[patchi]))
 373         {
 374             nInternalFaces += pbm[patchi].size();
 375         }
 376     }
 377
 378     // Create the adjncy array the size of the total number of internal and
 379     // coupled faces
 380     List<int> adjncy(nInternalFaces);
 381
 382     // Fill in xadj
 383     // ~~~~~~~~~~~~
 384     label freeAdj = 0;
 385
 386     for (label cellI = 0; cellI < mesh_.nCells(); cellI++)
 387     {
 388         xadj[cellI] = freeAdj;
 389
 390         const labelList& cFaces = mesh_.cells()[cellI];
 391
 392         forAll(cFaces, i)
 393         {
 394             label faceI = cFaces[i];
 395
 396             if
 397             (
 398                 mesh_.isInternalFace(faceI)
 399              || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
 400             )
 401             {
 402                 freeAdj++;
 403             }
 404         }
 405     }
 406     xadj[mesh_.nCells()] = freeAdj;
 407
 408
 409     // Fill in adjncy
 410     // ~~~~~~~~~~~~~~
 411
 412     labelList nFacesPerCell(mesh_.nCells(), 0);
 413
 414     // Internal faces
 415     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 416     {
 417         label own = mesh_.faceOwner()[faceI];
 418         label nei = mesh_.faceNeighbour()[faceI];
 419
 420         adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
 421         adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
 422     }
 423
 424     // Coupled faces. Only cyclics done.
 425     forAll(pbm, patchi)
 426     {
 427         if (isA<cyclicPolyPatch>(pbm[patchi]))
 428         {
 429             const unallocLabelList& faceCells = pbm[patchi].faceCells();
 430
 431             label sizeby2 = faceCells.size()/2;
 432
 433             for (label facei=0; facei<sizeby2; facei++)
 434             {
 435                 label own = faceCells[facei];
 436                 label nei = faceCells[facei + sizeby2];
 437
 438                 adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
 439                 adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
 440             }
 441         }
 442     }
 443
 444     // Decompose using default weights
 445     List<int> finalDecomp;
 446     decompose(adjncy, xadj, finalDecomp);
 447
 448     // Copy back to labelList
 449     labelList decomp(finalDecomp.size());
 450     forAll(decomp, i)
 451     {
 452         decomp[i] = finalDecomp[i];
 453     }
 454     return decomp;
 455 }
 456
 457
 458 // From cell-cell connections to Metis format (like CompactListList)
 459 void Foam::scotchDecomp::calcMetisCSR
 460 (
 461     const labelListList& cellCells,
 462     List<int>& adjncy,
 463     List<int>& xadj
 464 )
 465 {
 466     // Count number of internal faces
 467     label nConnections = 0;
 468
 469     forAll(cellCells, coarseI)
 470     {
 471         nConnections += cellCells[coarseI].size();
 472     }
 473
 474     // Create the adjncy array as twice the size of the total number of
 475     // internal faces
 476     adjncy.setSize(nConnections);
 477
 478     xadj.setSize(cellCells.size()+1);
 479
 480
 481     // Fill in xadj
 482     // ~~~~~~~~~~~~
 483     label freeAdj = 0;
 484
 485     forAll(cellCells, coarseI)
 486     {
 487         xadj[coarseI] = freeAdj;
 488
 489         const labelList& cCells = cellCells[coarseI];
 490
 491         forAll(cCells, i)
 492         {
 493             adjncy[freeAdj++] = cCells[i];
 494         }
 495     }
 496     xadj[cellCells.size()] = freeAdj;
 497 }
 498
 499
 500 Foam::labelList Foam::scotchDecomp::decompose
 501 (
 502     const labelList& agglom,
 503     const pointField& agglomPoints
 504 )
 505 {
 506     if (agglom.size() != mesh_.nCells())
 507     {
 508         FatalErrorIn
 509         (
 510             "parMetisDecomp::decompose(const labelList&, const pointField&)"
 511         )   << "Size of cell-to-coarse map " << agglom.size()
 512             << " differs from number of cells in mesh " << mesh_.nCells()
 513             << exit(FatalError);
 514     }
 515
 516     // Make Metis CSR (Compressed Storage Format) storage
 517     //   adjncy      : contains neighbours (= edges in graph)
 518     //   xadj(celli) : start of information in adjncy for celli
 519     List<int> adjncy;
 520     List<int> xadj;
 521     {
 522         // Get cellCells on coarse mesh.
 523         labelListList cellCells;
 524         calcCellCells
 525         (
 526             mesh_,
 527             agglom,
 528             agglomPoints.size(),
 529             cellCells
 530         );
 531
 532         calcMetisCSR(cellCells, adjncy, xadj);
 533     }
 534
 535     // Decompose using default weights
 536     List<int> finalDecomp;
 537     decompose(adjncy, xadj, finalDecomp);
 538
 539
 540     // Rework back into decomposition for original mesh_
 541     labelList fineDistribution(agglom.size());
 542
 543     forAll(fineDistribution, i)
 544     {
 545         fineDistribution[i] = finalDecomp[agglom[i]];
 546     }
 547
 548     return fineDistribution;
 549 }
 550
 551
 552 Foam::labelList Foam::scotchDecomp::decompose
 553 (
 554     const labelListList& globalCellCells,
 555     const pointField& cellCentres
 556 )
 557 {
 558     if (cellCentres.size() != globalCellCells.size())
 559     {
 560         FatalErrorIn
 561         (
 562             "scotchDecomp::decompose(const pointField&, const labelListList&)"
 563         )   << "Inconsistent number of cells (" << globalCellCells.size()
 564             << ") and number of cell centres (" << cellCentres.size()
 565             << ")." << exit(FatalError);
 566     }
 567
 568
 569     // Make Metis CSR (Compressed Storage Format) storage
 570     //   adjncy      : contains neighbours (= edges in graph)
 571     //   xadj(celli) : start of information in adjncy for celli
 572
 573     List<int> adjncy;
 574     List<int> xadj;
 575     calcMetisCSR(globalCellCells, adjncy, xadj);
 576
 577
 578     // Decompose using default weights
 579     List<int> finalDecomp;
 580     decompose(adjncy, xadj, finalDecomp);
 581
 582     // Copy back to labelList
 583     labelList decomp(finalDecomp.size());
 584     forAll(decomp, i)
 585     {
 586         decomp[i] = finalDecomp[i];
 587     }
 588     return decomp;
 589 }
 590
 591
 592 // ************************************************************************* //