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 #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
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 \
00052 \
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
00084 \
00085 switch (options_p[1]) \
00086 { \
00087 case 1: opts[METIS_OPTION_CTYPE] = METIS_CTYPE_RM; break; \
00088 case 2: \
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 \
00097 METIS_ ## BASE_NAME \
00098 ( \
00099 &nvtxs, \
00100 &ncon, \
00101 xadj.begin(), \
00102 adjncy.begin(), \
00103 vwgt_p ? vwgt.begin() : NULL, \
00104 NULL, \
00105 adjwgt_p ? adjwgt.begin() : NULL, \
00106 &nparts, \
00107 tpwgts_p ? tpwgts.begin() : NULL, \
00108 NULL, \
00109 opts.begin(), \
00110 &edgecut, \
00111 part.begin() \
00112 ); \
00113 \
00114 \
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 , \
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
00164
00165
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
00176 int numFlag = 0;
00177
00178
00179
00180
00181 word method("k-way");
00182
00183 int numCells = xadj.size()-1;
00184
00185
00186 List<int> options(5, 0);
00187
00188
00189
00190 Field<floatScalar> processorWeights;
00191
00192
00193 List<int> cellWeights;
00194
00195
00196 List<int> faceWeights;
00197
00198
00199
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
00224 cellWeights.setSize(cWeights.size());
00225 forAll(cellWeights, i)
00226 {
00227 cellWeights[i] = int(cWeights[i]/minWeights);
00228 }
00229 }
00230
00231
00232
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
00313 finalDecomp.setSize(numCells);
00314
00315
00316 int edgeCut = 0;
00317
00318
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;
00327 }
00328 if (faceWeights.size())
00329 {
00330 adjwgtPtr = faceWeights.begin();
00331 wgtFlag += 1;
00332 }
00333
00334 if (method == "recursive")
00335 {
00336 if (processorWeights.size())
00337 {
00338 FOAM_WPartGraphRecursive
00339 (
00340 &numCells,
00341 const_cast<List<int>&>(xadj).begin(),
00342 const_cast<List<int>&>(adjncy).begin(),
00343 vwgtPtr,
00344 adjwgtPtr,
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,
00359 const_cast<List<int>&>(xadj).begin(),
00360 const_cast<List<int>&>(adjncy).begin(),
00361 vwgtPtr,
00362 adjwgtPtr,
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,
00379 const_cast<List<int>&>(xadj).begin(),
00380 const_cast<List<int>&>(adjncy).begin(),
00381 vwgtPtr,
00382 adjwgtPtr,
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,
00397 const_cast<List<int>&>(xadj).begin(),
00398 const_cast<List<int>&>(adjncy).begin(),
00399 vwgtPtr,
00400 adjwgtPtr,
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
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
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
00459 List<int> finalDecomp;
00460 decompose(adjncy, xadj, pointWeights, finalDecomp);
00461
00462
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
00491
00492
00493 List<int> adjncy;
00494 List<int> xadj;
00495 {
00496
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
00510 List<int> finalDecomp;
00511 decompose(adjncy, xadj, agglomWeights, finalDecomp);
00512
00513
00514
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
00546
00547
00548
00549 List<int> adjncy;
00550 List<int> xadj;
00551 scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);
00552
00553
00554
00555 List<int> finalDecomp;
00556 decompose(adjncy, xadj, cellWeights, finalDecomp);
00557
00558
00559 labelList decomp(finalDecomp.size());
00560 forAll(decomp, i)
00561 {
00562 decomp[i] = finalDecomp[i];
00563 }
00564 return decomp;
00565 }
00566
00567
00568