FreeFOAM The Cross-Platform CFD Toolkit
metisDecomp.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
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 \*---------------------------------------------------------------------------*/
25 
26 #include "metisDecomp.H"
28 #include <OpenFOAM/floatScalar.H>
29 #include <OpenFOAM/Time.H>
31 
32 extern "C"
33 {
34 #define OMPI_SKIP_MPICXX
35 # include <metis.h>
36 }
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 #if defined(METIS_NEW_API) || (defined(METIS_API) && defined(METIS_NOPTIONS))
41 
42 // Compatibility functions implementing the old API
43 #define DEFINE_METIS_COMPAT_FUNC(BASE_NAME) \
44  void FOAM_W ## BASE_NAME(int *nvtxs_p, int *xadj_p, int *adjncy_p, \
45  int *vwgt_p, int *adjwgt_p, int *wgtflag_p, \
46  int *numflag_p, int *nparts_p, \
47  float *tpwgts_p, int *options_p, \
48  int *edgecut_p, int *part_p) \
49 { \
50  using namespace Foam; \
51  /* ignored: wgtflag_p */ \
52  /* convert arguments to correct type */ \
53  idx_t nvtxs = *nvtxs_p; \
54  idx_t ncon = 1; \
55  idx_t madjncy = xadj_p[nvtxs]; \
56  List<idx_t> xadj(nvtxs+1); \
57  std::copy(xadj_p, xadj_p+nvtxs+1, xadj.begin()); \
58  List<idx_t> adjncy(madjncy); \
59  std::copy(adjncy_p, adjncy_p+madjncy, adjncy.begin()); \
60  List<idx_t> vwgt; \
61  if (vwgt_p) \
62  { \
63  vwgt.setSize(nvtxs); \
64  std::copy(vwgt_p, vwgt_p+nvtxs, vwgt.begin()); \
65  } \
66  List<idx_t> adjwgt; \
67  if (adjwgt_p) \
68  { \
69  adjwgt.setSize(madjncy); \
70  std::copy(adjwgt_p, adjwgt_p+madjncy, adjwgt.begin()); \
71  } \
72  idx_t nparts = *nparts_p; \
73  List<real_t> tpwgts; \
74  if (tpwgts_p) \
75  { \
76  tpwgts.setSize(nparts); \
77  std::copy(tpwgts_p, tpwgts_p+nparts, tpwgts.begin()); \
78  } \
79  List<idx_t> opts(METIS_NOPTIONS); \
80  METIS_SetDefaultOptions(opts.begin()); \
81  if (options_p && options_p[0]) \
82  { \
83  /* only options_p[1] has multiple, non-default options which has \
84  * equivalent in new API. Not sure about refinement in Kway-case... */\
85  switch (options_p[1]) \
86  { \
87  case 1: opts[METIS_OPTION_CTYPE] = METIS_CTYPE_RM; break; \
88  case 2: /* fall-through */ \
89  case 3: opts[METIS_OPTION_CTYPE] = METIS_CTYPE_SHEM; break; \
90  } \
91  } \
92  opts[METIS_OPTION_NUMBERING] = *numflag_p; \
93  idx_t edgecut; \
94  List<idx_t> part(nvtxs); \
95  \
96  /* call new API function */ \
97  METIS_ ## BASE_NAME \
98  ( \
99  &nvtxs, \
100  &ncon, \
101  xadj.begin(), \
102  adjncy.begin(), \
103  vwgt_p ? vwgt.begin() : NULL, \
104  NULL, /* vsize */ \
105  adjwgt_p ? adjwgt.begin() : NULL, \
106  &nparts, \
107  tpwgts_p ? tpwgts.begin() : NULL, \
108  NULL, /* ubvec */ \
109  opts.begin(), \
110  &edgecut, \
111  part.begin() \
112  ); \
113  \
114  /* pass output values back */ \
115  *edgecut_p = edgecut; \
116  forAll(part, partI) \
117  { \
118  part_p[partI] = part[partI]; \
119  } \
120 } \
121  \
122  \
123 inline \
124 void FOAM_ ## BASE_NAME(int *nvtxs_p, int *xadj_p, int *adjncy_p, int *vwgt_p,\
125  int *adjwgt_p, int *wgtflag_p, int *numflag_p, \
126  int *nparts_p, int *options_p, int *edgecut_p, \
127  int *part_p) \
128 { \
129  FOAM_W ## BASE_NAME(nvtxs_p, xadj_p, adjncy_p, vwgt_p, adjwgt_p, \
130  wgtflag_p, numflag_p, nparts_p, NULL /* tpwgts */, \
131  options_p, edgecut_p, part_p); \
132 }
133 
134 namespace {
135  DEFINE_METIS_COMPAT_FUNC(PartGraphRecursive)
136  DEFINE_METIS_COMPAT_FUNC(PartGraphKway)
137 }
138 
139 #else
140 
141 #define FOAM_PartGraphRecursive METIS_PartGraphRecursive
142 #define FOAM_WPartGraphRecursive METIS_WPartGraphRecursive
143 #define FOAM_PartGraphKway METIS_PartGraphKway
144 #define FOAM_WPartGraphKway METIS_WPartGraphKway
145 
146 #endif
147 
148 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
149 
150 namespace Foam
151 {
152  defineTypeNameAndDebug(metisDecomp, 0);
153 
155  (
156  decompositionMethod,
157  metisDecomp,
158  dictionaryMesh
159  );
160 }
161 
162 
163 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
164 
165 // Call Metis with options from dictionary.
166 Foam::label Foam::metisDecomp::decompose
167 (
168  const List<int>& adjncy,
169  const List<int>& xadj,
170  const scalarField& cWeights,
171 
172  List<int>& finalDecomp
173 )
174 {
175  // C style numbering
176  int numFlag = 0;
177 
178  // Method of decomposition
179  // recursive: multi-level recursive bisection (default)
180  // k-way: multi-level k-way
181  word method("k-way");
182 
183  int numCells = xadj.size()-1;
184 
185  // decomposition options. 0 = use defaults
186  List<int> options(5, 0);
187 
188  // processor weights initialised with no size, only used if specified in
189  // a file
190  Field<floatScalar> processorWeights;
191 
192  // cell weights (so on the vertices of the dual)
193  List<int> cellWeights;
194 
195  // face weights (so on the edges of the dual)
196  List<int> faceWeights;
197 
198 
199  // Check for externally provided cellweights and if so initialise weights
200  scalar minWeights = gMin(cWeights);
201  if (cWeights.size() > 0)
202  {
203  if (minWeights <= 0)
204  {
205  WarningIn
206  (
207  "metisDecomp::decompose"
208  "(const pointField&, const scalarField&)"
209  ) << "Illegal minimum weight " << minWeights
210  << endl;
211  }
212 
213  if (cWeights.size() != numCells)
214  {
216  (
217  "metisDecomp::decompose"
218  "(const pointField&, const scalarField&)"
219  ) << "Number of cell weights " << cWeights.size()
220  << " does not equal number of cells " << numCells
221  << exit(FatalError);
222  }
223  // Convert to integers.
224  cellWeights.setSize(cWeights.size());
225  forAll(cellWeights, i)
226  {
227  cellWeights[i] = int(cWeights[i]/minWeights);
228  }
229  }
230 
231 
232  // Check for user supplied weights and decomp options
233  if (decompositionDict_.found("metisCoeffs"))
234  {
235  const dictionary& metisCoeffs =
236  decompositionDict_.subDict("metisCoeffs");
237  word weightsFile;
238 
239  if (metisCoeffs.readIfPresent("method", method))
240  {
241  if (method != "recursive" && method != "k-way")
242  {
243  FatalErrorIn("metisDecomp::decompose()")
244  << "Method " << method << " in metisCoeffs in dictionary : "
245  << decompositionDict_.name()
246  << " should be 'recursive' or 'k-way'"
247  << exit(FatalError);
248  }
249 
250  Info<< "metisDecomp : Using Metis method " << method
251  << nl << endl;
252  }
253 
254  if (metisCoeffs.readIfPresent("options", options))
255  {
256  if (options.size() != 5)
257  {
258  FatalErrorIn("metisDecomp::decompose()")
259  << "Number of options in metisCoeffs in dictionary : "
260  << decompositionDict_.name()
261  << " should be 5"
262  << exit(FatalError);
263  }
264 
265  Info<< "metisDecomp : Using Metis options " << options
266  << nl << endl;
267  }
268 
269  if (metisCoeffs.readIfPresent("processorWeights", processorWeights))
270  {
271  processorWeights /= sum(processorWeights);
272 
273  if (processorWeights.size() != nProcessors_)
274  {
275  FatalErrorIn("metisDecomp::decompose(const pointField&)")
276  << "Number of processor weights "
277  << processorWeights.size()
278  << " does not equal number of domains " << nProcessors_
279  << exit(FatalError);
280  }
281  }
282 
283  if (metisCoeffs.readIfPresent("cellWeightsFile", weightsFile))
284  {
285  Info<< "metisDecomp : Using cell-based weights." << endl;
286 
287  IOList<int> cellIOWeights
288  (
289  IOobject
290  (
291  weightsFile,
292  mesh_.time().timeName(),
293  mesh_,
294  IOobject::MUST_READ,
295  IOobject::AUTO_WRITE
296  )
297  );
298  cellWeights.transfer(cellIOWeights);
299 
300  if (cellWeights.size() != xadj.size()-1)
301  {
302  FatalErrorIn("metisDecomp::decompose(const pointField&)")
303  << "Number of cell weights " << cellWeights.size()
304  << " does not equal number of cells " << xadj.size()-1
305  << exit(FatalError);
306  }
307  }
308  }
309 
310  int nProcs = nProcessors_;
311 
312  // output: cell -> processor addressing
313  finalDecomp.setSize(numCells);
314 
315  // output: number of cut edges
316  int edgeCut = 0;
317 
318  // Vertex weight info
319  int wgtFlag = 0;
320  int* vwgtPtr = NULL;
321  int* adjwgtPtr = NULL;
322 
323  if (cellWeights.size())
324  {
325  vwgtPtr = cellWeights.begin();
326  wgtFlag += 2; // Weights on vertices
327  }
328  if (faceWeights.size())
329  {
330  adjwgtPtr = faceWeights.begin();
331  wgtFlag += 1; // Weights on edges
332  }
333 
334  if (method == "recursive")
335  {
336  if (processorWeights.size())
337  {
338  FOAM_WPartGraphRecursive
339  (
340  &numCells, // num vertices in graph
341  const_cast<List<int>&>(xadj).begin(), // indexing into adjncy
342  const_cast<List<int>&>(adjncy).begin(), // neighbour info
343  vwgtPtr, // vertexweights
344  adjwgtPtr, // no edgeweights
345  &wgtFlag,
346  &numFlag,
347  &nProcs,
348  processorWeights.begin(),
349  options.begin(),
350  &edgeCut,
351  finalDecomp.begin()
352  );
353  }
354  else
355  {
356  FOAM_PartGraphRecursive
357  (
358  &numCells, // num vertices in graph
359  const_cast<List<int>&>(xadj).begin(), // indexing into adjncy
360  const_cast<List<int>&>(adjncy).begin(), // neighbour info
361  vwgtPtr, // vertexweights
362  adjwgtPtr, // no edgeweights
363  &wgtFlag,
364  &numFlag,
365  &nProcs,
366  options.begin(),
367  &edgeCut,
368  finalDecomp.begin()
369  );
370  }
371  }
372  else
373  {
374  if (processorWeights.size())
375  {
376  FOAM_WPartGraphKway
377  (
378  &numCells, // num vertices in graph
379  const_cast<List<int>&>(xadj).begin(), // indexing into adjncy
380  const_cast<List<int>&>(adjncy).begin(), // neighbour info
381  vwgtPtr, // vertexweights
382  adjwgtPtr, // no edgeweights
383  &wgtFlag,
384  &numFlag,
385  &nProcs,
386  processorWeights.begin(),
387  options.begin(),
388  &edgeCut,
389  finalDecomp.begin()
390  );
391  }
392  else
393  {
394  FOAM_PartGraphKway
395  (
396  &numCells, // num vertices in graph
397  const_cast<List<int>&>(xadj).begin(), // indexing into adjncy
398  const_cast<List<int>&>(adjncy).begin(), // neighbour info
399  vwgtPtr, // vertexweights
400  adjwgtPtr, // no edgeweights
401  &wgtFlag,
402  &numFlag,
403  &nProcs,
404  options.begin(),
405  &edgeCut,
406  finalDecomp.begin()
407  );
408  }
409  }
410 
411  return edgeCut;
412 }
413 
414 
415 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
416 
417 Foam::metisDecomp::metisDecomp
418 (
419  const dictionary& decompositionDict,
420  const polyMesh& mesh
421 )
422 :
423  decompositionMethod(decompositionDict),
424  mesh_(mesh)
425 {}
426 
427 
428 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
429 
430 Foam::labelList Foam::metisDecomp::decompose
431 (
432  const pointField& points,
433  const scalarField& pointWeights
434 )
435 {
436  if (points.size() != mesh_.nCells())
437  {
439  (
440  "metisDecomp::decompose(const pointField&,const scalarField&)"
441  ) << "Can use this decomposition method only for the whole mesh"
442  << endl
443  << "and supply one coordinate (cellCentre) for every cell." << endl
444  << "The number of coordinates " << points.size() << endl
445  << "The number of cells in the mesh " << mesh_.nCells()
446  << exit(FatalError);
447  }
448 
449  List<int> adjncy;
450  List<int> xadj;
452  (
453  mesh_,
454  adjncy,
455  xadj
456  );
457 
458  // Decompose using default weights
459  List<int> finalDecomp;
460  decompose(adjncy, xadj, pointWeights, finalDecomp);
461 
462  // Copy back to labelList
463  labelList decomp(finalDecomp.size());
464  forAll(decomp, i)
465  {
466  decomp[i] = finalDecomp[i];
467  }
468  return decomp;
469 }
470 
471 
472 Foam::labelList Foam::metisDecomp::decompose
473 (
474  const labelList& agglom,
475  const pointField& agglomPoints,
476  const scalarField& agglomWeights
477 )
478 {
479  if (agglom.size() != mesh_.nCells())
480  {
482  (
483  "metisDecomp::decompose"
484  "(const labelList&, const pointField&, const scalarField&)"
485  ) << "Size of cell-to-coarse map " << agglom.size()
486  << " differs from number of cells in mesh " << mesh_.nCells()
487  << exit(FatalError);
488  }
489 
490  // Make Metis CSR (Compressed Storage Format) storage
491  // adjncy : contains neighbours (= edges in graph)
492  // xadj(celli) : start of information in adjncy for celli
493  List<int> adjncy;
494  List<int> xadj;
495  {
496  // Get cellCells on coarse mesh.
497  labelListList cellCells;
498  calcCellCells
499  (
500  mesh_,
501  agglom,
502  agglomPoints.size(),
503  cellCells
504  );
505 
506  scotchDecomp::calcCSR(cellCells, adjncy, xadj);
507  }
508 
509  // Decompose using default weights
510  List<int> finalDecomp;
511  decompose(adjncy, xadj, agglomWeights, finalDecomp);
512 
513 
514  // Rework back into decomposition for original mesh_
515  labelList fineDistribution(agglom.size());
516 
517  forAll(fineDistribution, i)
518  {
519  fineDistribution[i] = finalDecomp[agglom[i]];
520  }
521 
522  return fineDistribution;
523 }
524 
525 
526 Foam::labelList Foam::metisDecomp::decompose
527 (
528  const labelListList& globalCellCells,
529  const pointField& cellCentres,
530  const scalarField& cellWeights
531 )
532 {
533  if (cellCentres.size() != globalCellCells.size())
534  {
536  (
537  "metisDecomp::decompose"
538  "(const pointField&, const labelListList&, const scalarField&)"
539  ) << "Inconsistent number of cells (" << globalCellCells.size()
540  << ") and number of cell centres (" << cellCentres.size()
541  << ")." << exit(FatalError);
542  }
543 
544 
545  // Make Metis CSR (Compressed Storage Format) storage
546  // adjncy : contains neighbours (= edges in graph)
547  // xadj(celli) : start of information in adjncy for celli
548 
549  List<int> adjncy;
550  List<int> xadj;
551  scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);
552 
553 
554  // Decompose using default weights
555  List<int> finalDecomp;
556  decompose(adjncy, xadj, cellWeights, finalDecomp);
557 
558  // Copy back to labelList
559  labelList decomp(finalDecomp.size());
560  forAll(decomp, i)
561  {
562  decomp[i] = finalDecomp[i];
563  }
564  return decomp;
565 }
566 
567 
568 // ************************ vim: set sw=4 sts=4 et: ************************ //