00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107 #include "scotchDecomp.H"
00108 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00109 #include <OpenFOAM/floatScalar.H>
00110 #include <OpenFOAM/Time.H>
00111 #include <OpenFOAM/cyclicPolyPatch.H>
00112 #include <OpenFOAM/OFstream.H>
00113
00114 #include "scotch.h"
00115
00116
00117
00118
00119 #if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
00120 # define LINUX
00121 #endif
00122
00123 #if defined(LINUX) && defined(__GNUC__)
00124 # define LINUX_GNUC
00125 #endif
00126
00127 #ifdef LINUX_GNUC
00128 # ifndef __USE_GNU
00129 # define __USE_GNU
00130 # endif
00131 # include <fenv.h>
00132 #endif
00133
00134
00135
00136
00137 namespace Foam
00138 {
00139 defineTypeNameAndDebug(scotchDecomp, 0);
00140
00141 addToRunTimeSelectionTable
00142 (
00143 decompositionMethod,
00144 scotchDecomp,
00145 dictionaryMesh
00146 );
00147 }
00148
00149
00150
00151 void Foam::scotchDecomp::check(const int retVal, const char* str)
00152 {
00153 if (retVal)
00154 {
00155 FatalErrorIn("scotchDecomp::decompose(..)")
00156 << "Call to scotch routine " << str << " failed."
00157 << exit(FatalError);
00158 }
00159 }
00160
00161
00162
00163 Foam::label Foam::scotchDecomp::decompose
00164 (
00165 const List<int>& adjncy,
00166 const List<int>& xadj,
00167 const scalarField& cWeights,
00168
00169 List<int>& finalDecomp
00170 )
00171 {
00172
00173 if (decompositionDict_.found("scotchCoeffs"))
00174 {
00175 const dictionary& scotchCoeffs =
00176 decompositionDict_.subDict("scotchCoeffs");
00177
00178 if (scotchCoeffs.found("writeGraph"))
00179 {
00180 Switch writeGraph(scotchCoeffs.lookup("writeGraph"));
00181
00182 if (writeGraph)
00183 {
00184 OFstream str(mesh_.time().path() / mesh_.name() + ".grf");
00185
00186 Info<< "Dumping Scotch graph file to " << str.name() << endl
00187 << "Use this in combination with gpart." << endl;
00188
00189 label version = 0;
00190 str << version << nl;
00191
00192 str << xadj.size()-1 << ' ' << adjncy.size() << nl;
00193
00194 label baseval = 0;
00195
00196 label hasEdgeWeights = 0;
00197 label hasVertexWeights = 0;
00198 label numericflag = 10*hasEdgeWeights+hasVertexWeights;
00199 str << baseval << ' ' << numericflag << nl;
00200 for (label cellI = 0; cellI < xadj.size()-1; cellI++)
00201 {
00202 label start = xadj[cellI];
00203 label end = xadj[cellI+1];
00204 str << end-start;
00205
00206 for (label i = start; i < end; i++)
00207 {
00208 str << ' ' << adjncy[i];
00209 }
00210 str << nl;
00211 }
00212 }
00213 }
00214 }
00215
00216
00217
00218
00219
00220
00221 SCOTCH_Strat stradat;
00222 check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
00223
00224 if (decompositionDict_.found("scotchCoeffs"))
00225 {
00226 const dictionary& scotchCoeffs =
00227 decompositionDict_.subDict("scotchCoeffs");
00228
00229
00230 string strategy;
00231 if (scotchCoeffs.readIfPresent("strategy", strategy))
00232 {
00233 if (debug)
00234 {
00235 Info<< "scotchDecomp : Using strategy " << strategy << endl;
00236 }
00237 SCOTCH_stratGraphMap(&stradat, strategy.c_str());
00238
00239
00240
00241 }
00242 }
00243
00244
00245
00246
00247
00248 List<int> velotab;
00249
00250
00251
00252 scalar minWeights = gMin(cWeights);
00253 if (cWeights.size() > 0)
00254 {
00255 if (minWeights <= 0)
00256 {
00257 WarningIn
00258 (
00259 "scotchDecomp::decompose"
00260 "(const pointField&, const scalarField&)"
00261 ) << "Illegal minimum weight " << minWeights
00262 << endl;
00263 }
00264
00265 if (cWeights.size() != xadj.size()-1)
00266 {
00267 FatalErrorIn
00268 (
00269 "scotchDecomp::decompose"
00270 "(const pointField&, const scalarField&)"
00271 ) << "Number of cell weights " << cWeights.size()
00272 << " does not equal number of cells " << xadj.size()-1
00273 << exit(FatalError);
00274 }
00275
00276
00277 velotab.setSize(cWeights.size());
00278 forAll(velotab, i)
00279 {
00280 velotab[i] = int(cWeights[i]/minWeights);
00281 }
00282 }
00283
00284
00285
00286 SCOTCH_Graph grafdat;
00287 check(SCOTCH_graphInit(&grafdat), "SCOTCH_graphInit");
00288 check
00289 (
00290 SCOTCH_graphBuild
00291 (
00292 &grafdat,
00293 0,
00294 xadj.size()-1,
00295 xadj.begin(),
00296 &xadj[1],
00297 velotab.begin(),
00298 NULL,
00299 adjncy.size(),
00300 adjncy.begin(),
00301 NULL
00302 ),
00303 "SCOTCH_graphBuild"
00304 );
00305 check(SCOTCH_graphCheck(&grafdat), "SCOTCH_graphCheck");
00306
00307
00308
00309
00310
00311
00312 SCOTCH_Arch archdat;
00313 check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
00314
00315 List<label> processorWeights;
00316 if (decompositionDict_.found("scotchCoeffs"))
00317 {
00318 const dictionary& scotchCoeffs =
00319 decompositionDict_.subDict("scotchCoeffs");
00320
00321 scotchCoeffs.readIfPresent("processorWeights", processorWeights);
00322 }
00323 if (processorWeights.size())
00324 {
00325 if (debug)
00326 {
00327 Info<< "scotchDecomp : Using procesor weights " << processorWeights
00328 << endl;
00329 }
00330 check
00331 (
00332 SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
00333 "SCOTCH_archCmpltw"
00334 );
00335 }
00336 else
00337 {
00338 check
00339 (
00340 SCOTCH_archCmplt(&archdat, nProcessors_),
00341 "SCOTCH_archCmplt"
00342 );
00343 }
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353 # ifdef LINUX_GNUC
00354 int oldExcepts = fedisableexcept
00355 (
00356 FE_DIVBYZERO
00357 | FE_INVALID
00358 | FE_OVERFLOW
00359 );
00360 # endif
00361
00362 finalDecomp.setSize(xadj.size()-1);
00363 finalDecomp = 0;
00364 check
00365 (
00366 SCOTCH_graphMap
00367 (
00368 &grafdat,
00369 &archdat,
00370 &stradat,
00371 finalDecomp.begin()
00372 ),
00373 "SCOTCH_graphMap"
00374 );
00375
00376 # ifdef LINUX_GNUC
00377 feenableexcept(oldExcepts);
00378 # endif
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396 SCOTCH_graphExit(&grafdat);
00397
00398 SCOTCH_stratExit(&stradat);
00399
00400 SCOTCH_archExit(&archdat);
00401
00402 return 0;
00403 }
00404
00405
00406
00407
00408 Foam::scotchDecomp::scotchDecomp
00409 (
00410 const dictionary& decompositionDict,
00411 const polyMesh& mesh
00412 )
00413 :
00414 decompositionMethod(decompositionDict),
00415 mesh_(mesh)
00416 {}
00417
00418
00419
00420
00421 Foam::labelList Foam::scotchDecomp::decompose
00422 (
00423 const pointField& points,
00424 const scalarField& pointWeights
00425 )
00426 {
00427 if (points.size() != mesh_.nCells())
00428 {
00429 FatalErrorIn
00430 (
00431 "scotchDecomp::decompose(const pointField&, const scalarField&)"
00432 )
00433 << "Can use this decomposition method only for the whole mesh"
00434 << endl
00435 << "and supply one coordinate (cellCentre) for every cell." << endl
00436 << "The number of coordinates " << points.size() << endl
00437 << "The number of cells in the mesh " << mesh_.nCells()
00438 << exit(FatalError);
00439 }
00440
00441
00442
00443
00444 List<int> adjncy;
00445 List<int> xadj;
00446 calcCSR(mesh_, adjncy, xadj);
00447
00448
00449 List<int> finalDecomp;
00450 decompose(adjncy, xadj, pointWeights, finalDecomp);
00451
00452
00453 labelList decomp(finalDecomp.size());
00454 forAll(decomp, i)
00455 {
00456 decomp[i] = finalDecomp[i];
00457 }
00458 return decomp;
00459 }
00460
00461
00462 Foam::labelList Foam::scotchDecomp::decompose
00463 (
00464 const labelList& agglom,
00465 const pointField& agglomPoints,
00466 const scalarField& pointWeights
00467 )
00468 {
00469 if (agglom.size() != mesh_.nCells())
00470 {
00471 FatalErrorIn
00472 (
00473 "parMetisDecomp::decompose(const labelList&, const pointField&)"
00474 ) << "Size of cell-to-coarse map " << agglom.size()
00475 << " differs from number of cells in mesh " << mesh_.nCells()
00476 << exit(FatalError);
00477 }
00478
00479
00480
00481
00482 List<int> adjncy;
00483 List<int> xadj;
00484 {
00485
00486 labelListList cellCells;
00487 calcCellCells
00488 (
00489 mesh_,
00490 agglom,
00491 agglomPoints.size(),
00492 cellCells
00493 );
00494
00495 calcCSR(cellCells, adjncy, xadj);
00496 }
00497
00498
00499 List<int> finalDecomp;
00500 decompose(adjncy, xadj, pointWeights, finalDecomp);
00501
00502
00503 labelList fineDistribution(agglom.size());
00504
00505 forAll(fineDistribution, i)
00506 {
00507 fineDistribution[i] = finalDecomp[agglom[i]];
00508 }
00509
00510 return fineDistribution;
00511 }
00512
00513
00514 Foam::labelList Foam::scotchDecomp::decompose
00515 (
00516 const labelListList& globalCellCells,
00517 const pointField& cellCentres,
00518 const scalarField& cWeights
00519 )
00520 {
00521 if (cellCentres.size() != globalCellCells.size())
00522 {
00523 FatalErrorIn
00524 (
00525 "scotchDecomp::decompose"
00526 "(const labelListList&, const pointField&, const scalarField&)"
00527 ) << "Inconsistent number of cells (" << globalCellCells.size()
00528 << ") and number of cell centres (" << cellCentres.size()
00529 << ")." << exit(FatalError);
00530 }
00531
00532
00533
00534
00535
00536
00537 List<int> adjncy;
00538 List<int> xadj;
00539 calcCSR(globalCellCells, adjncy, xadj);
00540
00541
00542 List<int> finalDecomp;
00543 decompose(adjncy, xadj, cWeights, finalDecomp);
00544
00545
00546 labelList decomp(finalDecomp.size());
00547 forAll(decomp, i)
00548 {
00549 decomp[i] = finalDecomp[i];
00550 }
00551 return decomp;
00552 }
00553
00554
00555 void Foam::scotchDecomp::calcCSR
00556 (
00557 const polyMesh& mesh,
00558 List<int>& adjncy,
00559 List<int>& xadj
00560 )
00561 {
00562
00563
00564
00565
00566 xadj.setSize(mesh.nCells()+1);
00567
00568
00569
00570 label nInternalFaces = 2*mesh.nInternalFaces();
00571
00572
00573
00574 const polyBoundaryMesh& pbm = mesh.boundaryMesh();
00575
00576 forAll(pbm, patchi)
00577 {
00578 if (isA<cyclicPolyPatch>(pbm[patchi]))
00579 {
00580 nInternalFaces += pbm[patchi].size();
00581 }
00582 }
00583
00584
00585
00586 adjncy.setSize(nInternalFaces);
00587
00588
00589
00590 label freeAdj = 0;
00591
00592 for (label cellI = 0; cellI < mesh.nCells(); cellI++)
00593 {
00594 xadj[cellI] = freeAdj;
00595
00596 const labelList& cFaces = mesh.cells()[cellI];
00597
00598 forAll(cFaces, i)
00599 {
00600 label faceI = cFaces[i];
00601
00602 if
00603 (
00604 mesh.isInternalFace(faceI)
00605 || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
00606 )
00607 {
00608 freeAdj++;
00609 }
00610 }
00611 }
00612 xadj[mesh.nCells()] = freeAdj;
00613
00614
00615
00616
00617
00618 labelList nFacesPerCell(mesh.nCells(), 0);
00619
00620
00621 for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
00622 {
00623 label own = mesh.faceOwner()[faceI];
00624 label nei = mesh.faceNeighbour()[faceI];
00625
00626 adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
00627 adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
00628 }
00629
00630
00631 forAll(pbm, patchi)
00632 {
00633 if (isA<cyclicPolyPatch>(pbm[patchi]))
00634 {
00635 const unallocLabelList& faceCells = pbm[patchi].faceCells();
00636
00637 label sizeby2 = faceCells.size()/2;
00638
00639 for (label facei=0; facei<sizeby2; facei++)
00640 {
00641 label own = faceCells[facei];
00642 label nei = faceCells[facei + sizeby2];
00643
00644 adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
00645 adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
00646 }
00647 }
00648 }
00649 }
00650
00651
00652
00653 void Foam::scotchDecomp::calcCSR
00654 (
00655 const labelListList& cellCells,
00656 List<int>& adjncy,
00657 List<int>& xadj
00658 )
00659 {
00660
00661 label nConnections = 0;
00662
00663 forAll(cellCells, coarseI)
00664 {
00665 nConnections += cellCells[coarseI].size();
00666 }
00667
00668
00669
00670 adjncy.setSize(nConnections);
00671
00672 xadj.setSize(cellCells.size()+1);
00673
00674
00675
00676
00677 label freeAdj = 0;
00678
00679 forAll(cellCells, coarseI)
00680 {
00681 xadj[coarseI] = freeAdj;
00682
00683 const labelList& cCells = cellCells[coarseI];
00684
00685 forAll(cCells, i)
00686 {
00687 adjncy[freeAdj++] = cCells[i];
00688 }
00689 }
00690 xadj[cellCells.size()] = freeAdj;
00691 }
00692
00693
00694
00695