src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.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     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     Note: writes out .dgr files for debugging. Run with e.g.
 107
 108         mpirun -np 4 dgpart 2 'processor%r.grf'
 109
 110     - %r gets replaced by current processor rank
 111     - decompose into 2 domains
 112
 113 \*---------------------------------------------------------------------------*/
 114
 115 #include "ptscotchDecomp.H"
 116 #include "addToRunTimeSelectionTable.H"
 117 #include "Time.H"
 118 #include "OFstream.H"
 119 #include "globalIndex.H"
 120 #include "SubField.H"
 121
 122 extern "C"
 123 {
 124 #include <stdio.h>
 125 #include "mpi.h"
 126 #include "ptscotch.h"
 127 }
 128
 129
 130 // Hack: scotch generates floating point errors so need to switch of error
 131 //       trapping!
 132 #if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
 133 #    define LINUX
 134 #endif
 135
 136 #if defined(LINUX) && defined(__GNUC__)
 137 #    define LINUX_GNUC
 138 #endif
 139
 140 #ifdef LINUX_GNUC
 141 #   ifndef __USE_GNU
 142 #       define __USE_GNU
 143 #   endif
 144 #   include <fenv.h>
 145 #endif
 146
 147
 148 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 149
 150 namespace Foam
 151 {
 152     defineTypeNameAndDebug(ptscotchDecomp, 0);
 153
 154     addToRunTimeSelectionTable(decompositionMethod, ptscotchDecomp, dictionary);
 155 }
 156
 157 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 158
 159 void Foam::ptscotchDecomp::check(const int retVal, const char* str)
 160 {
 161     if (retVal)
 162     {
 163         FatalErrorIn("ptscotchDecomp::decompose(..)")
 164             << "Call to scotch routine " << str << " failed."
 165             << exit(FatalError);
 166     }
 167 }
 168
 169
 170 //- Does prevention of 0 cell domains and calls ptscotch.
 171 Foam::label Foam::ptscotchDecomp::decomposeZeroDomains
 172 (
 173     const fileName& meshPath,
 174     const List<int>& initadjncy,
 175     const List<int>& initxadj,
 176     const scalarField& initcWeights,
 177
 178     List<int>& finalDecomp
 179 ) const
 180 {
 181     globalIndex globalCells(initxadj.size()-1);
 182
 183     bool hasZeroDomain = false;
 184     for (label procI = 0; procI < Pstream::nProcs(); procI++)
 185     {
 186         if (globalCells.localSize(procI) == 0)
 187         {
 188             hasZeroDomain = true;
 189             break;
 190         }
 191     }
 192
 193     if (!hasZeroDomain)
 194     {
 195         return decompose
 196         (
 197             meshPath,
 198             initadjncy,
 199             initxadj,
 200             initcWeights,
 201             finalDecomp
 202         );
 203     }
 204
 205
 206     if (debug)
 207     {
 208         Info<< "ptscotchDecomp : have graphs with locally 0 cells."
 209             << " trickling down." << endl;
 210     }
 211
 212     // Make sure every domain has at least one cell
 213     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 214     // (scotch does not like zero sized domains)
 215     // Trickle cells from processors that have them up to those that
 216     // don't.
 217
 218
 219     // Number of cells to send to the next processor
 220     // (is same as number of cells next processor has to receive)
 221     List<int> nSendCells(Pstream::nProcs(), 0);
 222
 223     for (label procI = nSendCells.size()-1; procI >=1; procI--)
 224     {
 225         label nLocalCells = globalCells.localSize(procI);
 226         if (nLocalCells-nSendCells[procI] < 1)
 227         {
 228             nSendCells[procI-1] = nSendCells[procI]-nLocalCells+1;
 229         }
 230     }
 231
 232     // First receive (so increasing the sizes of all arrays)
 233
 234     Field<int> xadj(initxadj);
 235     Field<int> adjncy(initadjncy);
 236     scalarField cWeights(initcWeights);
 237
 238     if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
 239     {
 240         // Receive cells from previous processor
 241         IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
 242
 243         Field<int> prevXadj(fromPrevProc);
 244         Field<int> prevAdjncy(fromPrevProc);
 245         scalarField prevCellWeights(fromPrevProc);
 246
 247         if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1])
 248         {
 249             FatalErrorIn("ptscotchDecomp::decompose(..)")
 250                 << "Expected from processor " << Pstream::myProcNo()-1
 251                 << " connectivity for " << nSendCells[Pstream::myProcNo()-1]
 252                 << " nCells but only received " << prevXadj.size()
 253                 << abort(FatalError);
 254         }
 255
 256         // Insert adjncy
 257         prepend(prevAdjncy, adjncy);
 258         // Adapt offsets and prepend xadj
 259         xadj += prevAdjncy.size();
 260         prepend(prevXadj, xadj);
 261         // Weights
 262         prepend(prevCellWeights, cWeights);
 263     }
 264
 265
 266     // Send to my next processor
 267
 268     if (nSendCells[Pstream::myProcNo()] > 0)
 269     {
 270         // Send cells to next processor
 271         OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1);
 272
 273         int nCells = nSendCells[Pstream::myProcNo()];
 274         int startCell = xadj.size()-1 - nCells;
 275         int startFace = xadj[startCell];
 276         int nFaces = adjncy.size()-startFace;
 277
 278         // Send for all cell data: last nCells elements
 279         // Send for all face data: last nFaces elements
 280         toNextProc
 281             << Field<int>::subField(xadj, nCells, startCell)-startFace
 282             << Field<int>::subField(adjncy, nFaces, startFace)
 283             <<
 284             (
 285                 cWeights.size()
 286               ? static_cast<const scalarField&>
 287                 (
 288                     scalarField::subField(cWeights, nCells, startCell)
 289                 )
 290               : scalarField(0)
 291             );
 292
 293         // Remove data that has been sent
 294         if (cWeights.size())
 295         {
 296             cWeights.setSize(cWeights.size()-nCells);
 297         }
 298         adjncy.setSize(adjncy.size()-nFaces);
 299         xadj.setSize(xadj.size() - nCells);
 300     }
 301
 302
 303     // Do decomposition as normal. Sets finalDecomp.
 304     label result = decompose(meshPath, adjncy, xadj, cWeights, finalDecomp);
 305
 306
 307     if (debug)
 308     {
 309         Info<< "ptscotchDecomp : have graphs with locally 0 cells."
 310             << " trickling up." << endl;
 311     }
 312
 313
 314     // If we sent cells across make sure we undo it
 315     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 316
 317     // Receive back from next processor if I sent something
 318     if (nSendCells[Pstream::myProcNo()] > 0)
 319     {
 320         IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1);
 321
 322         List<int> nextFinalDecomp(fromNextProc);
 323
 324         if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()])
 325         {
 326             FatalErrorIn("parMetisDecomp::decompose(..)")
 327                 << "Expected from processor " << Pstream::myProcNo()+1
 328                 << " decomposition for " << nSendCells[Pstream::myProcNo()]
 329                 << " nCells but only received " << nextFinalDecomp.size()
 330                 << abort(FatalError);
 331         }
 332
 333         append(nextFinalDecomp, finalDecomp);
 334     }
 335
 336     // Send back to previous processor.
 337     if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
 338     {
 339         OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
 340
 341         int nToPrevious = nSendCells[Pstream::myProcNo()-1];
 342
 343         toPrevProc <<
 344             SubList<int>
 345             (
 346                 finalDecomp,
 347                 nToPrevious,
 348                 finalDecomp.size()-nToPrevious
 349             );
 350
 351         // Remove locally what has been sent
 352         finalDecomp.setSize(finalDecomp.size()-nToPrevious);
 353     }
 354     return result;
 355 }
 356
 357
 358 // Call scotch with options from dictionary.
 359 Foam::label Foam::ptscotchDecomp::decompose
 360 (
 361     const fileName& meshPath,
 362     const List<int>& adjncy,
 363     const List<int>& xadj,
 364     const scalarField& cWeights,
 365
 366     List<int>& finalDecomp
 367 ) const
 368 {
 369     if (debug)
 370     {
 371         Pout<< "ptscotchDecomp : entering with xadj:" << xadj.size() << endl;
 372     }
 373
 374     // Dump graph
 375     if (decompositionDict_.found("ptscotchCoeffs"))
 376     {
 377         const dictionary& scotchCoeffs =
 378             decompositionDict_.subDict("ptscotchCoeffs");
 379
 380         if (scotchCoeffs.lookupOrDefault("writeGraph", false))
 381         {
 382             OFstream str
 383             (
 384                meshPath + "_" + Foam::name(Pstream::myProcNo()) + ".dgr"
 385             );
 386
 387             Pout<< "Dumping Scotch graph file to " << str.name() << endl
 388                 << "Use this in combination with dgpart." << endl;
 389
 390             globalIndex globalCells(xadj.size()-1);
 391
 392             // Distributed graph file (.grf)
 393             label version = 2;
 394             str << version << nl;
 395             // Number of files (procglbnbr)
 396             str << Pstream::nProcs();
 397             // My file number (procloc)
 398             str << ' ' << Pstream::myProcNo() << nl;
 399
 400             // Total number of vertices (vertglbnbr)
 401             str << globalCells.size();
 402             // Total number of connections (edgeglbnbr)
 403             str << ' ' << returnReduce(xadj[xadj.size()-1], sumOp<label>())
 404                 << nl;
 405             // Local number of vertices (vertlocnbr)
 406             str << xadj.size()-1;
 407             // Local number of connections (edgelocnbr)
 408             str << ' ' << xadj[xadj.size()-1] << nl;
 409             // Numbering starts from 0
 410             label baseval = 0;
 411             // 100*hasVertlabels+10*hasEdgeWeights+1*hasVertWeighs
 412             str << baseval << ' ' << "000" << nl;
 413             for (label cellI = 0; cellI < xadj.size()-1; cellI++)
 414             {
 415                 label start = xadj[cellI];
 416                 label end = xadj[cellI+1];
 417                 str << end-start;
 418
 419                 for (label i = start; i < end; i++)
 420                 {
 421                     str << ' ' << adjncy[i];
 422                 }
 423                 str << nl;
 424             }
 425         }
 426     }
 427
 428     // Strategy
 429     // ~~~~~~~~
 430
 431     // Default.
 432     SCOTCH_Strat stradat;
 433     check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
 434
 435     if (decompositionDict_.found("scotchCoeffs"))
 436     {
 437         const dictionary& scotchCoeffs =
 438             decompositionDict_.subDict("scotchCoeffs");
 439
 440
 441         string strategy;
 442         if (scotchCoeffs.readIfPresent("strategy", strategy))
 443         {
 444             if (debug)
 445             {
 446                 Info<< "ptscotchDecomp : Using strategy " << strategy << endl;
 447             }
 448             SCOTCH_stratDgraphMap(&stradat, strategy.c_str());
 449             //fprintf(stdout, "S\tStrat=");
 450             //SCOTCH_stratSave(&stradat, stdout);
 451             //fprintf(stdout, "\n");
 452         }
 453     }
 454
 455
 456     // Graph
 457     // ~~~~~
 458
 459     List<int> velotab;
 460
 461
 462     // Check for externally provided cellweights and if so initialise weights
 463     scalar minWeights = gMin(cWeights);
 464     if (cWeights.size() > 0)
 465     {
 466         if (minWeights <= 0)
 467         {
 468             WarningIn
 469             (
 470                 "ptscotchDecomp::decompose(..)"
 471             )   << "Illegal minimum weight " << minWeights
 472                 << endl;
 473         }
 474
 475         if (cWeights.size() != xadj.size()-1)
 476         {
 477             FatalErrorIn
 478             (
 479                 "ptscotchDecomp::decompose(..)"
 480             )   << "Number of cell weights " << cWeights.size()
 481                 << " does not equal number of cells " << xadj.size()-1
 482                 << exit(FatalError);
 483         }
 484
 485         // Convert to integers.
 486         velotab.setSize(cWeights.size());
 487         forAll(velotab, i)
 488         {
 489             velotab[i] = int(cWeights[i]/minWeights);
 490         }
 491     }
 492
 493
 494
 495     if (debug)
 496     {
 497         Pout<< "SCOTCH_dgraphInit" << endl;
 498     }
 499     SCOTCH_Dgraph grafdat;
 500     check(SCOTCH_dgraphInit(&grafdat, MPI_COMM_WORLD), "SCOTCH_dgraphInit");
 501
 502
 503     if (debug)
 504     {
 505         Pout<< "SCOTCH_dgraphBuild with:" << nl
 506             << "xadj.size()-1   : " << xadj.size()-1 << nl
 507             << "xadj            : " << long(xadj.begin()) << nl
 508             << "velotab         : " << long(velotab.begin()) << nl
 509             << "adjncy.size()   : " << adjncy.size() << nl
 510             << "adjncy          : " << long(adjncy.begin()) << nl
 511             << endl;
 512     }
 513
 514     check
 515     (
 516         SCOTCH_dgraphBuild
 517         (
 518             &grafdat,               // grafdat
 519             0,                      // baseval, c-style numbering
 520             xadj.size()-1,          // vertlocnbr, nCells
 521             xadj.size()-1,          // vertlocmax
 522             const_cast<SCOTCH_Num*>(xadj.begin()),
 523                                     // vertloctab, start index per cell into
 524                                     // adjncy
 525             const_cast<SCOTCH_Num*>(&xadj[1]),// vendloctab, end index  ,,
 526
 527             const_cast<SCOTCH_Num*>(velotab.begin()),// veloloctab, vtx weights
 528             NULL,                   // vlblloctab
 529
 530             adjncy.size(),          // edgelocnbr, number of arcs
 531             adjncy.size(),          // edgelocsiz
 532             const_cast<SCOTCH_Num*>(adjncy.begin()),         // edgeloctab
 533             NULL,                   // edgegsttab
 534             NULL                    // edlotab, edge weights
 535         ),
 536         "SCOTCH_dgraphBuild"
 537     );
 538
 539
 540     if (debug)
 541     {
 542         Pout<< "SCOTCH_dgraphCheck" << endl;
 543     }
 544     check(SCOTCH_dgraphCheck(&grafdat), "SCOTCH_dgraphCheck");
 545
 546
 547     // Architecture
 548     // ~~~~~~~~~~~~
 549     // (fully connected network topology since using switch)
 550
 551     if (debug)
 552     {
 553         Pout<< "SCOTCH_archInit" << endl;
 554     }
 555     SCOTCH_Arch archdat;
 556     check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
 557
 558     List<label> processorWeights;
 559     if (decompositionDict_.found("scotchCoeffs"))
 560     {
 561         const dictionary& scotchCoeffs =
 562             decompositionDict_.subDict("scotchCoeffs");
 563
 564         scotchCoeffs.readIfPresent("processorWeights", processorWeights);
 565     }
 566     if (processorWeights.size())
 567     {
 568         if (debug)
 569         {
 570             Info<< "ptscotchDecomp : Using procesor weights "
 571                 << processorWeights
 572                 << endl;
 573         }
 574         check
 575         (
 576             SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
 577             "SCOTCH_archCmpltw"
 578         );
 579     }
 580     else
 581     {
 582         if (debug)
 583         {
 584             Pout<< "SCOTCH_archCmplt" << endl;
 585         }
 586         check
 587         (
 588             SCOTCH_archCmplt(&archdat, nProcessors_),
 589             "SCOTCH_archCmplt"
 590         );
 591     }
 592
 593
 594     //SCOTCH_Mapping mapdat;
 595     //SCOTCH_dgraphMapInit(&grafdat, &mapdat, &archdat, NULL);
 596     //SCOTCH_dgraphMapCompute(&grafdat, &mapdat, &stradat); /*Perform mapping*/
 597     //SCOTCHdgraphMapExit(&grafdat, &mapdat);
 598
 599
 600     // Hack:switch off fpu error trapping
 601 #   ifdef LINUX_GNUC
 602     int oldExcepts = fedisableexcept
 603     (
 604         FE_DIVBYZERO
 605       | FE_INVALID
 606       | FE_OVERFLOW
 607     );
 608 #   endif
 609
 610     if (debug)
 611     {
 612         Pout<< "SCOTCH_dgraphMap" << endl;
 613     }
 614     finalDecomp.setSize(xadj.size()-1);
 615     finalDecomp = 0;
 616     check
 617     (
 618         SCOTCH_dgraphMap
 619         (
 620             &grafdat,
 621             &archdat,
 622             &stradat,           // const SCOTCH_Strat *
 623             finalDecomp.begin() // parttab
 624         ),
 625         "SCOTCH_graphMap"
 626     );
 627
 628 #   ifdef LINUX_GNUC
 629     feenableexcept(oldExcepts);
 630 #   endif
 631
 632
 633
 634     //finalDecomp.setSize(xadj.size()-1);
 635     //check
 636     //(
 637     //    SCOTCH_dgraphPart
 638     //    (
 639     //        &grafdat,
 640     //        nProcessors_,       // partnbr
 641     //        &stradat,           // const SCOTCH_Strat *
 642     //        finalDecomp.begin() // parttab
 643     //    ),
 644     //    "SCOTCH_graphPart"
 645     //);
 646
 647     if (debug)
 648     {
 649         Pout<< "SCOTCH_dgraphExit" << endl;
 650     }
 651     // Release storage for graph
 652     SCOTCH_dgraphExit(&grafdat);
 653     // Release storage for strategy
 654     SCOTCH_stratExit(&stradat);
 655     // Release storage for network topology
 656     SCOTCH_archExit(&archdat);
 657
 658     return 0;
 659 }
 660
 661
 662 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 663
 664 Foam::ptscotchDecomp::ptscotchDecomp(const dictionary& decompositionDict)
 665 :
 666     decompositionMethod(decompositionDict)
 667 {}
 668
 669
 670 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 671
 672 Foam::labelList Foam::ptscotchDecomp::decompose
 673 (
 674     const polyMesh& mesh,
 675     const pointField& points,
 676     const scalarField& pointWeights
 677 )
 678 {
 679     if (points.size() != mesh.nCells())
 680     {
 681         FatalErrorIn
 682         (
 683             "ptscotchDecomp::decompose(const pointField&, const scalarField&)"
 684         )
 685             << "Can use this decomposition method only for the whole mesh"
 686             << endl
 687             << "and supply one coordinate (cellCentre) for every cell." << endl
 688             << "The number of coordinates " << points.size() << endl
 689             << "The number of cells in the mesh " << mesh.nCells()
 690             << exit(FatalError);
 691     }
 692
 693 //    // For running sequential ...
 694 //    if (Pstream::nProcs() <= 1)
 695 //    {
 696 //        return scotchDecomp(decompositionDict_, mesh_)
 697 //            .decompose(points, pointWeights);
 698 //    }
 699
 700     // Make Metis CSR (Compressed Storage Format) storage
 701     //   adjncy      : contains neighbours (= edges in graph)
 702     //   xadj(celli) : start of information in adjncy for celli
 703
 704
 705     CompactListList<label> cellCells;
 706     calcCellCells(mesh, identity(mesh.nCells()), mesh.nCells(), cellCells);
 707
 708     // Decompose using default weights
 709     List<int> finalDecomp;
 710     decomposeZeroDomains
 711     (
 712         mesh.time().path()/mesh.name(),
 713         cellCells.m(),
 714         cellCells.offsets(),
 715         pointWeights,
 716         finalDecomp
 717     );
 718
 719     // Copy back to labelList
 720     labelList decomp(finalDecomp.size());
 721     forAll(decomp, i)
 722     {
 723         decomp[i] = finalDecomp[i];
 724     }
 725     return decomp;
 726 }
 727
 728
 729 Foam::labelList Foam::ptscotchDecomp::decompose
 730 (
 731     const polyMesh& mesh,
 732     const labelList& agglom,
 733     const pointField& agglomPoints,
 734     const scalarField& pointWeights
 735 )
 736 {
 737     if (agglom.size() != mesh.nCells())
 738     {
 739         FatalErrorIn
 740         (
 741             "ptscotchDecomp::decompose(const labelList&, const pointField&)"
 742         )   << "Size of cell-to-coarse map " << agglom.size()
 743             << " differs from number of cells in mesh " << mesh.nCells()
 744             << exit(FatalError);
 745     }
 746
 747 //    // For running sequential ...
 748 //    if (Pstream::nProcs() <= 1)
 749 //    {
 750 //        return scotchDecomp(decompositionDict_, mesh)
 751 //            .decompose(agglom, agglomPoints, pointWeights);
 752 //    }
 753
 754     // Make Metis CSR (Compressed Storage Format) storage
 755     //   adjncy      : contains neighbours (= edges in graph)
 756     //   xadj(celli) : start of information in adjncy for celli
 757     CompactListList<label> cellCells;
 758     calcCellCells(mesh, agglom, agglomPoints.size(), cellCells);
 759
 760     // Decompose using weights
 761     List<int> finalDecomp;
 762     decomposeZeroDomains
 763     (
 764         mesh.time().path()/mesh.name(),
 765         cellCells.m(),
 766         cellCells.offsets(),
 767         pointWeights,
 768         finalDecomp
 769     );
 770
 771     // Rework back into decomposition for original mesh
 772     labelList fineDistribution(agglom.size());
 773
 774     forAll(fineDistribution, i)
 775     {
 776         fineDistribution[i] = finalDecomp[agglom[i]];
 777     }
 778
 779     return fineDistribution;
 780 }
 781
 782
 783 Foam::labelList Foam::ptscotchDecomp::decompose
 784 (
 785     const labelListList& globalCellCells,
 786     const pointField& cellCentres,
 787     const scalarField& cWeights
 788 )
 789 {
 790     if (cellCentres.size() != globalCellCells.size())
 791     {
 792         FatalErrorIn
 793         (
 794             "ptscotchDecomp::decompose(const pointField&, const labelListList&)"
 795         )   << "Inconsistent number of cells (" << globalCellCells.size()
 796             << ") and number of cell centres (" << cellCentres.size()
 797             << ")." << exit(FatalError);
 798     }
 799
 800 //    // For running sequential ...
 801 //    if (Pstream::nProcs() <= 1)
 802 //    {
 803 //        return scotchDecomp(decompositionDict_, mesh)
 804 //            .decompose(globalCellCells, cellCentres, cWeights);
 805 //    }
 806
 807
 808     // Make Metis CSR (Compressed Storage Format) storage
 809     //   adjncy      : contains neighbours (= edges in graph)
 810     //   xadj(celli) : start of information in adjncy for celli
 811
 812     CompactListList<label> cellCells(globalCellCells);
 813
 814     // Decompose using weights
 815     List<int> finalDecomp;
 816     decomposeZeroDomains
 817     (
 818         "ptscotch",
 819         cellCells.m(),
 820         cellCells.offsets(),
 821         cWeights,
 822         finalDecomp
 823     );
 824
 825     // Copy back to labelList
 826     labelList decomp(finalDecomp.size());
 827     forAll(decomp, i)
 828     {
 829         decomp[i] = finalDecomp[i];
 830     }
 831     return decomp;
 832 }
 833
 834
 835 // ************************************************************************* //