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) 2011-2016 OpenFOAM Foundation
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 "Time.H"
29 
30 extern "C"
31 {
32  #define OMPI_SKIP_MPICXX
33  #include "metis.h"
34 }
35 
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41  defineTypeNameAndDebug(metisDecomp, 0);
42  addToRunTimeSelectionTable(decompositionMethod, metisDecomp, dictionary);
43 }
44 
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 Foam::label Foam::metisDecomp::decompose
49 (
50  const List<label>& adjncy,
51  const List<label>& xadj,
52  const scalarField& cWeights,
53 
54  List<label>& finalDecomp
55 )
56 {
57  // Method of decomposition
58  // recursive: multi-level recursive bisection (default)
59  // k-way: multi-level k-way
60  word method("recursive");
61 
62  label numCells = xadj.size()-1;
63 
64  // Decomposition options
65  List<label> options(METIS_NOPTIONS);
66  METIS_SetDefaultOptions(options.begin());
67 
68  // Processor weights initialised with no size, only used if specified in
69  // a file
70  Field<scalar> processorWeights;
71 
72  // Cell weights (so on the vertices of the dual)
73  List<label> cellWeights;
74 
75  // Face weights (so on the edges of the dual)
76  List<label> faceWeights;
77 
78 
79  // Check for externally provided cellweights and if so initialise weights
80  scalar minWeights = gMin(cWeights);
81  if (cWeights.size() > 0)
82  {
83  if (minWeights <= 0)
84  {
86  << "Illegal minimum weight " << minWeights
87  << endl;
88  }
89 
90  if (cWeights.size() != numCells)
91  {
93  << "Number of cell weights " << cWeights.size()
94  << " does not equal number of cells " << numCells
95  << exit(FatalError);
96  }
97 
98  // Convert to integers.
99  cellWeights.setSize(cWeights.size());
100  forAll(cellWeights, i)
101  {
102  cellWeights[i] = int(cWeights[i]/minWeights);
103  }
104  }
105 
106 
107  // Check for user supplied weights and decomp options
108  if (decompositionDict_.found("metisCoeffs"))
109  {
110  const dictionary& metisCoeffs =
111  decompositionDict_.subDict("metisCoeffs");
112 
113  word weightsFile;
114 
115  if (metisCoeffs.readIfPresent("method", method))
116  {
117  if (method != "recursive" && method != "k-way")
118  {
120  << "Method " << method << " in metisCoeffs in dictionary : "
121  << decompositionDict_.name()
122  << " should be 'recursive' or 'k-way'"
123  << exit(FatalError);
124  }
125 
126  Info<< "metisDecomp : Using Metis method " << method
127  << nl << endl;
128  }
129 
130  if (metisCoeffs.readIfPresent("options", options))
131  {
132  if (options.size() != METIS_NOPTIONS)
133  {
135  << "Number of options in metisCoeffs in dictionary : "
136  << decompositionDict_.name()
137  << " should be " << METIS_NOPTIONS
138  << exit(FatalError);
139  }
140 
141  Info<< "metisDecomp : Using Metis options " << options
142  << nl << endl;
143  }
144 
145  if (metisCoeffs.readIfPresent("processorWeights", processorWeights))
146  {
147  processorWeights /= sum(processorWeights);
148 
149  if (processorWeights.size() != nProcessors_)
150  {
152  << "Number of processor weights "
153  << processorWeights.size()
154  << " does not equal number of domains " << nProcessors_
155  << exit(FatalError);
156  }
157  }
158  }
159 
160  label ncon = 1;
161  label nProcs = nProcessors_;
162 
163  // Output: cell -> processor addressing
164  finalDecomp.setSize(numCells);
165 
166  // Output: number of cut edges
167  label edgeCut = 0;
168 
169  if (method == "recursive")
170  {
171  METIS_PartGraphRecursive
172  (
173  &numCells, // num vertices in graph
174  &ncon, // num balancing constraints
175  const_cast<List<label>&>(xadj).begin(), // indexing into adjncy
176  const_cast<List<label>&>(adjncy).begin(), // neighbour info
177  cellWeights.begin(),// vertexweights
178  NULL, // vsize: total communication vol
179  faceWeights.begin(),// edgeweights
180  &nProcs, // nParts
181  processorWeights.begin(), // tpwgts
182  NULL, // ubvec: processor imbalance (default)
183  options.begin(),
184  &edgeCut,
185  finalDecomp.begin()
186  );
187  }
188  else
189  {
190  METIS_PartGraphKway
191  (
192  &numCells, // num vertices in graph
193  &ncon, // num balancing constraints
194  const_cast<List<label>&>(xadj).begin(), // indexing into adjncy
195  const_cast<List<label>&>(adjncy).begin(), // neighbour info
196  cellWeights.begin(),// vertexweights
197  NULL, // vsize: total communication vol
198  faceWeights.begin(),// edgeweights
199  &nProcs, // nParts
200  processorWeights.begin(), // tpwgts
201  NULL, // ubvec: processor imbalance (default)
202  options.begin(),
203  &edgeCut,
204  finalDecomp.begin()
205  );
206  }
207 
208  return edgeCut;
209 }
210 
211 
212 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
213 
214 Foam::metisDecomp::metisDecomp(const dictionary& decompositionDict)
215 :
216  decompositionMethod(decompositionDict)
217 {}
218 
219 
220 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
221 
222 Foam::labelList Foam::metisDecomp::decompose
223 (
224  const polyMesh& mesh,
225  const pointField& points,
226  const scalarField& pointWeights
227 )
228 {
229  if (points.size() != mesh.nCells())
230  {
232  << "Can use this decomposition method only for the whole mesh"
233  << endl
234  << "and supply one coordinate (cellCentre) for every cell." << endl
235  << "The number of coordinates " << points.size() << endl
236  << "The number of cells in the mesh " << mesh.nCells()
237  << exit(FatalError);
238  }
239 
240  CompactListList<label> cellCells;
241  calcCellCells
242  (
243  mesh,
244  identity(mesh.nCells()),
245  mesh.nCells(),
246  false,
247  cellCells
248  );
249 
250  // Decompose using default weights
251  labelList decomp;
252  decompose(cellCells.m(), cellCells.offsets(), pointWeights, decomp);
253 
254  return decomp;
255 }
256 
257 
258 Foam::labelList Foam::metisDecomp::decompose
259 (
260  const polyMesh& mesh,
261  const labelList& agglom,
262  const pointField& agglomPoints,
263  const scalarField& agglomWeights
264 )
265 {
266  if (agglom.size() != mesh.nCells())
267  {
269  << "Size of cell-to-coarse map " << agglom.size()
270  << " differs from number of cells in mesh " << mesh.nCells()
271  << exit(FatalError);
272  }
273 
274  // Make Metis CSR (Compressed Storage Format) storage
275  // adjncy : contains neighbours (= edges in graph)
276  // xadj(celli) : start of information in adjncy for celli
277 
278  CompactListList<label> cellCells;
279  calcCellCells(mesh, agglom, agglomPoints.size(), false, cellCells);
280 
281  // Decompose using default weights
282  labelList finalDecomp;
283  decompose(cellCells.m(), cellCells.offsets(), agglomWeights, finalDecomp);
284 
285 
286  // Rework back into decomposition for original mesh
287  labelList fineDistribution(agglom.size());
288 
289  forAll(fineDistribution, i)
290  {
291  fineDistribution[i] = finalDecomp[agglom[i]];
292  }
293 
294  return finalDecomp;
295 }
296 
297 
298 Foam::labelList Foam::metisDecomp::decompose
299 (
300  const labelListList& globalCellCells,
301  const pointField& cellCentres,
302  const scalarField& cellWeights
303 )
304 {
305  if (cellCentres.size() != globalCellCells.size())
306  {
308  << "Inconsistent number of cells (" << globalCellCells.size()
309  << ") and number of cell centres (" << cellCentres.size()
310  << ")." << exit(FatalError);
311  }
312 
313 
314  // Make Metis CSR (Compressed Storage Format) storage
315  // adjncy : contains neighbours (= edges in graph)
316  // xadj(celli) : start of information in adjncy for celli
317 
318  CompactListList<label> cellCells(globalCellCells);
319 
320  // Decompose using default weights
321  labelList decomp;
322  decompose(cellCells.m(), cellCells.offsets(), cellWeights, decomp);
323 
324  return decomp;
325 }
326 
327 
328 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Type gMin(const FieldField< Field, Type > &f)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
label nCells() const
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
dynamicFvMesh & mesh
const pointField & points
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< label > labelList
A List of labels.
Definition: labelList.H:56
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
static const char nl
Definition: Ostream.H:262
defineTypeNameAndDebug(combustionModel, 0)
#define WarningInFunction
Report a warning using Foam::Warning.
messageStream Info
Namespace for OpenFOAM.