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

metisDecomp.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 \*---------------------------------------------------------------------------*/
00025 
00026 #include "metisDecomp.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 #include <OpenFOAM/floatScalar.H>
00029 #include <OpenFOAM/Time.H>
00030 #include <scotchDecomp/scotchDecomp.H>
00031 
00032 extern "C"
00033 {
00034 #define OMPI_SKIP_MPICXX
00035 #   include <metis.h>
00036 }
00037 
00038 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00039 
00040 #if defined(METIS_NEW_API) || (defined(METIS_API) && defined(METIS_NOPTIONS))
00041 
00042 // Compatibility functions implementing the old API
00043 #define DEFINE_METIS_COMPAT_FUNC(BASE_NAME)                                   \
00044     void FOAM_W ## BASE_NAME(int *nvtxs_p, int *xadj_p, int *adjncy_p,        \
00045                                   int *vwgt_p, int *adjwgt_p, int *wgtflag_p, \
00046                                   int *numflag_p, int *nparts_p,              \
00047                                   float *tpwgts_p, int *options_p,            \
00048                                   int *edgecut_p, int *part_p)                \
00049 {                                                                             \
00050     using namespace Foam;                                                     \
00051     /* ignored: wgtflag_p */                                                  \
00052     /* convert arguments to correct type */                                   \
00053     idx_t nvtxs = *nvtxs_p;                                                   \
00054     idx_t ncon = 1;                                                           \
00055     idx_t madjncy = xadj_p[nvtxs];                                            \
00056     List<idx_t> xadj(nvtxs+1);                                                \
00057     std::copy(xadj_p, xadj_p+nvtxs+1, xadj.begin());                          \
00058     List<idx_t> adjncy(madjncy);                                              \
00059     std::copy(adjncy_p, adjncy_p+madjncy, adjncy.begin());                    \
00060     List<idx_t> vwgt;                                                         \
00061     if (vwgt_p)                                                               \
00062     {                                                                         \
00063         vwgt.setSize(nvtxs);                                                  \
00064         std::copy(vwgt_p, vwgt_p+nvtxs, vwgt.begin());                        \
00065     }                                                                         \
00066     List<idx_t> adjwgt;                                                       \
00067     if (adjwgt_p)                                                             \
00068     {                                                                         \
00069         adjwgt.setSize(madjncy);                                              \
00070         std::copy(adjwgt_p, adjwgt_p+madjncy, adjwgt.begin());                \
00071     }                                                                         \
00072     idx_t nparts = *nparts_p;                                                 \
00073     List<real_t> tpwgts;                                                      \
00074     if (tpwgts_p)                                                             \
00075     {                                                                         \
00076         tpwgts.setSize(nparts);                                               \
00077         std::copy(tpwgts_p, tpwgts_p+nparts, tpwgts.begin());                 \
00078     }                                                                         \
00079     List<idx_t> opts(METIS_NOPTIONS);                                         \
00080     METIS_SetDefaultOptions(opts.begin());                                    \
00081     if (options_p && options_p[0])                                            \
00082     {                                                                         \
00083         /* only options_p[1] has multiple, non-default options which has      \
00084          * equivalent in new API. Not sure about refinement in Kway-case... */\
00085         switch (options_p[1])                                                 \
00086         {                                                                     \
00087             case 1: opts[METIS_OPTION_CTYPE] = METIS_CTYPE_RM; break;         \
00088             case 2: /* fall-through */                                        \
00089             case 3: opts[METIS_OPTION_CTYPE] = METIS_CTYPE_SHEM; break;       \
00090         }                                                                     \
00091     }                                                                         \
00092     opts[METIS_OPTION_NUMBERING] = *numflag_p;                                \
00093     idx_t edgecut;                                                            \
00094     List<idx_t> part(nvtxs);                                                  \
00095                                                                               \
00096     /* call new API function */                                               \
00097     METIS_ ## BASE_NAME                                                       \
00098     (                                                                         \
00099         &nvtxs,                                                               \
00100         &ncon,                                                                \
00101         xadj.begin(),                                                         \
00102         adjncy.begin(),                                                       \
00103         vwgt_p ? vwgt.begin() : NULL,                                         \
00104         NULL,   /* vsize */                                                   \
00105         adjwgt_p ? adjwgt.begin() : NULL,                                     \
00106         &nparts,                                                              \
00107         tpwgts_p ? tpwgts.begin() : NULL,                                     \
00108         NULL,   /* ubvec */                                                   \
00109         opts.begin(),                                                         \
00110         &edgecut,                                                             \
00111         part.begin()                                                          \
00112     );                                                                        \
00113                                                                               \
00114     /* pass output values back */                                             \
00115     *edgecut_p = edgecut;                                                     \
00116     forAll(part, partI)                                                       \
00117     {                                                                         \
00118         part_p[partI] = part[partI];                                          \
00119     }                                                                         \
00120 }                                                                             \
00121                                                                               \
00122                                                                               \
00123 inline                                                                        \
00124 void FOAM_ ## BASE_NAME(int *nvtxs_p, int *xadj_p, int *adjncy_p, int *vwgt_p,\
00125                         int *adjwgt_p, int *wgtflag_p, int *numflag_p,        \
00126                         int *nparts_p, int *options_p, int *edgecut_p,        \
00127                         int *part_p)                                          \
00128 {                                                                             \
00129     FOAM_W ## BASE_NAME(nvtxs_p, xadj_p, adjncy_p, vwgt_p, adjwgt_p,          \
00130                         wgtflag_p, numflag_p, nparts_p, NULL /* tpwgts */,    \
00131                         options_p, edgecut_p, part_p);                        \
00132 }
00133 
00134 namespace {
00135     DEFINE_METIS_COMPAT_FUNC(PartGraphRecursive)
00136     DEFINE_METIS_COMPAT_FUNC(PartGraphKway)
00137 }
00138 
00139 #else
00140 
00141 #define FOAM_PartGraphRecursive  METIS_PartGraphRecursive
00142 #define FOAM_WPartGraphRecursive METIS_WPartGraphRecursive
00143 #define FOAM_PartGraphKway       METIS_PartGraphKway
00144 #define FOAM_WPartGraphKway      METIS_WPartGraphKway
00145 
00146 #endif
00147 
00148 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00149 
00150 namespace Foam
00151 {
00152     defineTypeNameAndDebug(metisDecomp, 0);
00153 
00154     addToRunTimeSelectionTable
00155     (
00156         decompositionMethod,
00157         metisDecomp,
00158         dictionaryMesh
00159     );
00160 }
00161 
00162 
00163 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00164 
00165 // Call Metis with options from dictionary.
00166 Foam::label Foam::metisDecomp::decompose
00167 (
00168     const List<int>& adjncy,
00169     const List<int>& xadj,
00170     const scalarField& cWeights,
00171 
00172     List<int>& finalDecomp
00173 )
00174 {
00175     // C style numbering
00176     int numFlag = 0;
00177 
00178     // Method of decomposition
00179     // recursive: multi-level recursive bisection (default)
00180     // k-way: multi-level k-way
00181     word method("k-way");
00182 
00183     int numCells = xadj.size()-1;
00184 
00185     // decomposition options. 0 = use defaults
00186     List<int> options(5, 0);
00187 
00188     // processor weights initialised with no size, only used if specified in
00189     // a file
00190     Field<floatScalar> processorWeights;
00191 
00192     // cell weights (so on the vertices of the dual)
00193     List<int> cellWeights;
00194 
00195     // face weights (so on the edges of the dual)
00196     List<int> faceWeights;
00197 
00198 
00199     // Check for externally provided cellweights and if so initialise weights
00200     scalar minWeights = gMin(cWeights);
00201     if (cWeights.size() > 0)
00202     {
00203         if (minWeights <= 0)
00204         {
00205             WarningIn
00206             (
00207                 "metisDecomp::decompose"
00208                 "(const pointField&, const scalarField&)"
00209             )   << "Illegal minimum weight " << minWeights
00210                 << endl;
00211         }
00212 
00213         if (cWeights.size() != numCells)
00214         {
00215             FatalErrorIn
00216             (
00217                 "metisDecomp::decompose"
00218                 "(const pointField&, const scalarField&)"
00219             )   << "Number of cell weights " << cWeights.size()
00220                 << " does not equal number of cells " << numCells
00221                 << exit(FatalError);
00222         }
00223         // Convert to integers.
00224         cellWeights.setSize(cWeights.size());
00225         forAll(cellWeights, i)
00226         {
00227             cellWeights[i] = int(cWeights[i]/minWeights);
00228         }
00229     }
00230 
00231 
00232     // Check for user supplied weights and decomp options
00233     if (decompositionDict_.found("metisCoeffs"))
00234     {
00235         const dictionary& metisCoeffs =
00236             decompositionDict_.subDict("metisCoeffs");
00237         word weightsFile;
00238 
00239         if (metisCoeffs.readIfPresent("method", method))
00240         {
00241             if (method != "recursive" && method != "k-way")
00242             {
00243                 FatalErrorIn("metisDecomp::decompose()")
00244                     << "Method " << method << " in metisCoeffs in dictionary : "
00245                     << decompositionDict_.name()
00246                     << " should be 'recursive' or 'k-way'"
00247                     << exit(FatalError);
00248             }
00249 
00250             Info<< "metisDecomp : Using Metis method     " << method
00251                 << nl << endl;
00252         }
00253 
00254         if (metisCoeffs.readIfPresent("options", options))
00255         {
00256             if (options.size() != 5)
00257             {
00258                 FatalErrorIn("metisDecomp::decompose()")
00259                     << "Number of options in metisCoeffs in dictionary : "
00260                     << decompositionDict_.name()
00261                     << " should be 5"
00262                     << exit(FatalError);
00263             }
00264 
00265             Info<< "metisDecomp : Using Metis options     " << options
00266                 << nl << endl;
00267         }
00268 
00269         if (metisCoeffs.readIfPresent("processorWeights", processorWeights))
00270         {
00271             processorWeights /= sum(processorWeights);
00272 
00273             if (processorWeights.size() != nProcessors_)
00274             {
00275                 FatalErrorIn("metisDecomp::decompose(const pointField&)")
00276                     << "Number of processor weights "
00277                     << processorWeights.size()
00278                     << " does not equal number of domains " << nProcessors_
00279                     << exit(FatalError);
00280             }
00281         }
00282 
00283         if (metisCoeffs.readIfPresent("cellWeightsFile", weightsFile))
00284         {
00285             Info<< "metisDecomp : Using cell-based weights." << endl;
00286 
00287             IOList<int> cellIOWeights
00288             (
00289                 IOobject
00290                 (
00291                     weightsFile,
00292                     mesh_.time().timeName(),
00293                     mesh_,
00294                     IOobject::MUST_READ,
00295                     IOobject::AUTO_WRITE
00296                 )
00297             );
00298             cellWeights.transfer(cellIOWeights);
00299 
00300             if (cellWeights.size() != xadj.size()-1)
00301             {
00302                 FatalErrorIn("metisDecomp::decompose(const pointField&)")
00303                     << "Number of cell weights " << cellWeights.size()
00304                     << " does not equal number of cells " << xadj.size()-1
00305                     << exit(FatalError);
00306             }
00307         }
00308     }
00309 
00310     int nProcs = nProcessors_;
00311 
00312     // output: cell -> processor addressing
00313     finalDecomp.setSize(numCells);
00314 
00315     // output: number of cut edges
00316     int edgeCut = 0;
00317 
00318     // Vertex weight info
00319     int wgtFlag = 0;
00320     int* vwgtPtr = NULL;
00321     int* adjwgtPtr = NULL;
00322 
00323     if (cellWeights.size())
00324     {
00325         vwgtPtr = cellWeights.begin();
00326         wgtFlag += 2;       // Weights on vertices
00327     }
00328     if (faceWeights.size())
00329     {
00330         adjwgtPtr = faceWeights.begin();
00331         wgtFlag += 1;       // Weights on edges
00332     }
00333 
00334     if (method == "recursive")
00335     {
00336         if (processorWeights.size())
00337         {
00338             FOAM_WPartGraphRecursive
00339             (
00340                 &numCells,         // num vertices in graph
00341                 const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
00342                 const_cast<List<int>&>(adjncy).begin(), // neighbour info
00343                 vwgtPtr,           // vertexweights
00344                 adjwgtPtr,         // no edgeweights
00345                 &wgtFlag,
00346                 &numFlag,
00347                 &nProcs,
00348                 processorWeights.begin(),
00349                 options.begin(),
00350                 &edgeCut,
00351                 finalDecomp.begin()
00352             );
00353         }
00354         else
00355         {
00356             FOAM_PartGraphRecursive
00357             (
00358                 &numCells,         // num vertices in graph
00359                 const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
00360                 const_cast<List<int>&>(adjncy).begin(), // neighbour info
00361                 vwgtPtr,           // vertexweights
00362                 adjwgtPtr,         // no edgeweights
00363                 &wgtFlag,
00364                 &numFlag,
00365                 &nProcs,
00366                 options.begin(),
00367                 &edgeCut,
00368                 finalDecomp.begin()
00369             );
00370         }
00371     }
00372     else
00373     {
00374         if (processorWeights.size())
00375         {
00376             FOAM_WPartGraphKway
00377             (
00378                 &numCells,         // num vertices in graph
00379                 const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
00380                 const_cast<List<int>&>(adjncy).begin(), // neighbour info
00381                 vwgtPtr,           // vertexweights
00382                 adjwgtPtr,         // no edgeweights
00383                 &wgtFlag,
00384                 &numFlag,
00385                 &nProcs,
00386                 processorWeights.begin(),
00387                 options.begin(),
00388                 &edgeCut,
00389                 finalDecomp.begin()
00390             );
00391         }
00392         else
00393         {
00394             FOAM_PartGraphKway
00395             (
00396                 &numCells,         // num vertices in graph
00397                 const_cast<List<int>&>(xadj).begin(),   // indexing into adjncy
00398                 const_cast<List<int>&>(adjncy).begin(), // neighbour info
00399                 vwgtPtr,           // vertexweights
00400                 adjwgtPtr,         // no edgeweights
00401                 &wgtFlag,
00402                 &numFlag,
00403                 &nProcs,
00404                 options.begin(),
00405                 &edgeCut,
00406                 finalDecomp.begin()
00407             );
00408         }
00409     }
00410 
00411     return edgeCut;
00412 }
00413 
00414 
00415 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00416 
00417 Foam::metisDecomp::metisDecomp
00418 (
00419     const dictionary& decompositionDict,
00420     const polyMesh& mesh
00421 )
00422 :
00423     decompositionMethod(decompositionDict),
00424     mesh_(mesh)
00425 {}
00426 
00427 
00428 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00429 
00430 Foam::labelList Foam::metisDecomp::decompose
00431 (
00432     const pointField& points,
00433     const scalarField& pointWeights
00434 )
00435 {
00436     if (points.size() != mesh_.nCells())
00437     {
00438         FatalErrorIn
00439         (
00440             "metisDecomp::decompose(const pointField&,const scalarField&)"
00441         )   << "Can use this decomposition method only for the whole mesh"
00442             << endl
00443             << "and supply one coordinate (cellCentre) for every cell." << endl
00444             << "The number of coordinates " << points.size() << endl
00445             << "The number of cells in the mesh " << mesh_.nCells()
00446             << exit(FatalError);
00447     }
00448 
00449     List<int> adjncy;
00450     List<int> xadj;
00451     scotchDecomp::calcCSR
00452     (
00453         mesh_,
00454         adjncy,
00455         xadj
00456     );
00457 
00458     // Decompose using default weights
00459     List<int> finalDecomp;
00460     decompose(adjncy, xadj, pointWeights, finalDecomp);
00461 
00462     // Copy back to labelList
00463     labelList decomp(finalDecomp.size());
00464     forAll(decomp, i)
00465     {
00466         decomp[i] = finalDecomp[i];
00467     }
00468     return decomp;
00469 }
00470 
00471 
00472 Foam::labelList Foam::metisDecomp::decompose
00473 (
00474     const labelList& agglom,
00475     const pointField& agglomPoints,
00476     const scalarField& agglomWeights
00477 )
00478 {
00479     if (agglom.size() != mesh_.nCells())
00480     {
00481         FatalErrorIn
00482         (
00483             "metisDecomp::decompose"
00484             "(const labelList&, const pointField&, const scalarField&)"
00485         )   << "Size of cell-to-coarse map " << agglom.size()
00486             << " differs from number of cells in mesh " << mesh_.nCells()
00487             << exit(FatalError);
00488     }
00489 
00490     // Make Metis CSR (Compressed Storage Format) storage
00491     //   adjncy      : contains neighbours (= edges in graph)
00492     //   xadj(celli) : start of information in adjncy for celli
00493     List<int> adjncy;
00494     List<int> xadj;
00495     {
00496         // Get cellCells on coarse mesh.
00497         labelListList cellCells;
00498         calcCellCells
00499         (
00500             mesh_,
00501             agglom,
00502             agglomPoints.size(),
00503             cellCells
00504         );
00505 
00506         scotchDecomp::calcCSR(cellCells, adjncy, xadj);
00507     }
00508 
00509     // Decompose using default weights
00510     List<int> finalDecomp;
00511     decompose(adjncy, xadj, agglomWeights, finalDecomp);
00512 
00513 
00514     // Rework back into decomposition for original mesh_
00515     labelList fineDistribution(agglom.size());
00516 
00517     forAll(fineDistribution, i)
00518     {
00519         fineDistribution[i] = finalDecomp[agglom[i]];
00520     }
00521 
00522     return fineDistribution;
00523 }
00524 
00525 
00526 Foam::labelList Foam::metisDecomp::decompose
00527 (
00528     const labelListList& globalCellCells,
00529     const pointField& cellCentres,
00530     const scalarField& cellWeights
00531 )
00532 {
00533     if (cellCentres.size() != globalCellCells.size())
00534     {
00535         FatalErrorIn
00536         (
00537             "metisDecomp::decompose"
00538             "(const pointField&, const labelListList&, const scalarField&)"
00539         )   << "Inconsistent number of cells (" << globalCellCells.size()
00540             << ") and number of cell centres (" << cellCentres.size()
00541             << ")." << exit(FatalError);
00542     }
00543 
00544 
00545     // Make Metis CSR (Compressed Storage Format) storage
00546     //   adjncy      : contains neighbours (= edges in graph)
00547     //   xadj(celli) : start of information in adjncy for celli
00548 
00549     List<int> adjncy;
00550     List<int> xadj;
00551     scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);
00552 
00553 
00554     // Decompose using default weights
00555     List<int> finalDecomp;
00556     decompose(adjncy, xadj, cellWeights, finalDecomp);
00557 
00558     // Copy back to labelList
00559     labelList decomp(finalDecomp.size());
00560     forAll(decomp, i)
00561     {
00562         decomp[i] = finalDecomp[i];
00563     }
00564     return decomp;
00565 }
00566 
00567 
00568 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines