FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

scotchDecomp.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024     From scotch forum:
00025     
00026     By: Francois PELLEGRINI RE: Graph mapping 'strategy' string [ reply ]  
00027     2008-08-22 10:09 Strategy handling in Scotch is a bit tricky. In order
00028     not to be confused, you must have a clear view of how they are built.
00029     Here are some rules:
00030 
00031     1- Strategies are made up of "methods" which are combined by means of
00032     "operators".
00033 
00034     2- A method is of the form "m{param=value,param=value,...}", where "m"
00035     is a single character (this is your first error: "f" is a method name,
00036     not a parameter name).
00037 
00038     3- There exist different sort of strategies : bipartitioning strategies,
00039     mapping strategies, ordering strategies, which cannot be mixed. For
00040     instance, you cannot build a bipartitioning strategy and feed it to a
00041     mapping method (this is your second error).
00042 
00043     To use the "mapCompute" routine, you must create a mapping strategy, not
00044     a bipartitioning one, and so use stratGraphMap() and not
00045     stratGraphBipart(). Your mapping strategy should however be based on the
00046     "recursive bipartitioning" method ("b"). For instance, a simple (and
00047     hence not very efficient) mapping strategy can be :
00048 
00049     "b{sep=f}"
00050 
00051     which computes mappings with the recursive bipartitioning method "b",
00052     this latter using the Fiduccia-Mattheyses method "f" to compute its
00053     separators.
00054 
00055     If you want an exact partition (see your previous post), try
00056     "b{sep=fx}".
00057 
00058     However, these strategies are not the most efficient, as they do not
00059     make use of the multi-level framework.
00060 
00061     To use the multi-level framework, try for instance:
00062 
00063     "b{sep=m{vert=100,low=h,asc=f}x}"
00064 
00065     The current default mapping strategy in Scotch can be seen by using the
00066     "-vs" option of program gmap. It is, to date:
00067 
00068     b
00069     {
00070         job=t,
00071         map=t,
00072         poli=S,
00073         sep=
00074         (
00075             m
00076             {
00077                 asc=b
00078                 {
00079                     bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
00080                     org=f{move=80,pass=-1,bal=0.005},
00081                     width=3
00082                 },
00083                 low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
00084                 type=h,
00085                 vert=80,
00086                 rat=0.8
00087             }
00088           | m
00089             {
00090                 asc=b
00091                 {
00092                     bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
00093                     org=f{move=80,pass=-1,bal=0.005},
00094                     width=3
00095                 },
00096                 low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
00097                 type=h,
00098                 vert=80,
00099                 rat=0.8
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 // Hack: scotch generates floating point errors so need to switch of error
00118 //       trapping!
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 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
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 // Call scotch with options from dictionary.
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     // Dump graph
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                 // Numer of vertices
00192                 str << xadj.size()-1 << ' ' << adjncy.size() << nl;
00193                 // Numbering starts from 0
00194                 label baseval = 0;
00195                 // Has weights?
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     // Strategy
00218     // ~~~~~~~~
00219 
00220     // Default.
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             //fprintf(stdout, "S\tStrat=");
00239             //SCOTCH_stratSave(&stradat, stdout);
00240             //fprintf(stdout, "\n");
00241         }
00242     }
00243 
00244 
00245     // Graph
00246     // ~~~~~
00247 
00248     List<int> velotab;
00249 
00250 
00251     // Check for externally provided cellweights and if so initialise weights
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         // Convert to integers.
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,                      // baseval, c-style numbering
00294             xadj.size()-1,          // vertnbr, nCells
00295             xadj.begin(),           // verttab, start index per cell into adjncy
00296             &xadj[1],               // vendtab, end index  ,,
00297             velotab.begin(),        // velotab, vertex weights
00298             NULL,                   // vlbltab
00299             adjncy.size(),          // edgenbr, number of arcs
00300             adjncy.begin(),         // edgetab
00301             NULL                    // edlotab, edge weights
00302         ),
00303         "SCOTCH_graphBuild"
00304     );
00305     check(SCOTCH_graphCheck(&grafdat), "SCOTCH_graphCheck");
00306 
00307 
00308     // Architecture
00309     // ~~~~~~~~~~~~
00310     // (fully connected network topology since using switch)
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     //SCOTCH_Mapping mapdat;
00347     //SCOTCH_graphMapInit(&grafdat, &mapdat, &archdat, NULL);
00348     //SCOTCH_graphMapCompute(&grafdat, &mapdat, &stradat); /* Perform mapping */
00349     //SCOTCH_graphMapExit(&grafdat, &mapdat);
00350 
00351 
00352     // Hack:switch off fpu error trapping
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,           // const SCOTCH_Strat *
00371             finalDecomp.begin() // parttab
00372         ),
00373         "SCOTCH_graphMap"
00374     );
00375 
00376 #   ifdef LINUX_GNUC
00377     feenableexcept(oldExcepts);
00378 #   endif
00379 
00380 
00381 
00382     //finalDecomp.setSize(xadj.size()-1);
00383     //check
00384     //(
00385     //    SCOTCH_graphPart
00386     //    (
00387     //        &grafdat,
00388     //        nProcessors_,       // partnbr
00389     //        &stradat,           // const SCOTCH_Strat *
00390     //        finalDecomp.begin() // parttab
00391     //    ),
00392     //    "SCOTCH_graphPart"
00393     //);
00394 
00395     // Release storage for graph
00396     SCOTCH_graphExit(&grafdat);
00397     // Release storage for strategy
00398     SCOTCH_stratExit(&stradat);
00399     // Release storage for network topology
00400     SCOTCH_archExit(&archdat);
00401 
00402     return 0;
00403 }
00404 
00405 
00406 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
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 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
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     // Make Metis CSR (Compressed Storage Format) storage
00442     //   adjncy      : contains neighbours (= edges in graph)
00443     //   xadj(celli) : start of information in adjncy for celli
00444     List<int> adjncy;
00445     List<int> xadj;
00446     calcCSR(mesh_, adjncy, xadj);
00447 
00448     // Decompose using default weights
00449     List<int> finalDecomp;
00450     decompose(adjncy, xadj, pointWeights, finalDecomp);
00451 
00452     // Copy back to labelList
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     // Make Metis CSR (Compressed Storage Format) storage
00480     //   adjncy      : contains neighbours (= edges in graph)
00481     //   xadj(celli) : start of information in adjncy for celli
00482     List<int> adjncy;
00483     List<int> xadj;
00484     {
00485         // Get cellCells on coarse mesh.
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     // Decompose using weights
00499     List<int> finalDecomp;
00500     decompose(adjncy, xadj, pointWeights, finalDecomp);
00501 
00502     // Rework back into decomposition for original mesh_
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     // Make Metis CSR (Compressed Storage Format) storage
00534     //   adjncy      : contains neighbours (= edges in graph)
00535     //   xadj(celli) : start of information in adjncy for celli
00536 
00537     List<int> adjncy;
00538     List<int> xadj;
00539     calcCSR(globalCellCells, adjncy, xadj);
00540 
00541     // Decompose using weights
00542     List<int> finalDecomp;
00543     decompose(adjncy, xadj, cWeights, finalDecomp);
00544 
00545     // Copy back to labelList
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     // Make Metis CSR (Compressed Storage Format) storage
00563     //   adjncy      : contains neighbours (= edges in graph)
00564     //   xadj(celli) : start of information in adjncy for celli
00565 
00566     xadj.setSize(mesh.nCells()+1);
00567 
00568     // Initialise the number of internal faces of the cells to twice the
00569     // number of internal faces
00570     label nInternalFaces = 2*mesh.nInternalFaces();
00571 
00572     // Check the boundary for coupled patches and add to the number of
00573     // internal faces
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     // Create the adjncy array the size of the total number of internal and
00585     // coupled faces
00586     adjncy.setSize(nInternalFaces);
00587 
00588     // Fill in xadj
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     // Fill in adjncy
00616     // ~~~~~~~~~~~~~~
00617 
00618     labelList nFacesPerCell(mesh.nCells(), 0);
00619 
00620     // Internal faces
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     // Coupled faces. Only cyclics done.
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 // From cell-cell connections to Metis format (like CompactListList)
00653 void Foam::scotchDecomp::calcCSR
00654 (
00655     const labelListList& cellCells,
00656     List<int>& adjncy,
00657     List<int>& xadj
00658 )
00659 {
00660     // Count number of internal faces
00661     label nConnections = 0;
00662 
00663     forAll(cellCells, coarseI)
00664     {
00665         nConnections += cellCells[coarseI].size();
00666     }
00667 
00668     // Create the adjncy array as twice the size of the total number of
00669     // internal faces
00670     adjncy.setSize(nConnections);
00671 
00672     xadj.setSize(cellCells.size()+1);
00673 
00674 
00675     // Fill in xadj
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 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines