scotchDecomp.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 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 "scotchDecomp.H"
28 #include "floatScalar.H"
29 #include "Time.H"
30 #include "OFstream.H"
31 #include "globalIndex.H"
32 #include "SubField.H"
33 
34 extern "C"
35 {
36 #include "scotch.h"
37 }
38 
39 
40 // Hack: scotch generates floating point errors so need to switch of error
41 // trapping!
42 #ifdef __GLIBC__
43  #ifndef _GNU_SOURCE
44  #define _GNU_SOURCE
45  #endif
46  #include <fenv.h>
47 #endif
48 
49 
50 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
51 
52 namespace Foam
53 {
54  defineTypeNameAndDebug(scotchDecomp, 0);
55 
57  (
58  decompositionMethod,
59  scotchDecomp,
60  decomposer
61  );
62 
64  (
65  decompositionMethod,
66  scotchDecomp,
67  distributor
68  );
69 }
70 
71 
72 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
73 
74 void Foam::scotchDecomp::check(const int retVal, const char* str)
75 {
76  if (retVal)
77  {
79  << "Call to scotch routine " << str << " failed."
80  << exit(FatalError);
81  }
82 }
83 
84 
85 Foam::label Foam::scotchDecomp::decompose
86 (
87  const fileName& meshPath,
88  const List<label>& adjncy,
89  const List<label>& xadj,
90  const scalarField& cWeights,
91 
92  List<label>& finalDecomp
93 )
94 {
95  if (!Pstream::parRun())
96  {
97  decomposeOneProc
98  (
99  meshPath,
100  adjncy,
101  xadj,
102  cWeights,
103  finalDecomp
104  );
105  }
106  else
107  {
108  if (debug)
109  {
110  Info<< "scotchDecomp : running in parallel."
111  << " Decomposing all of graph on master processor." << endl;
112  }
113  globalIndex globalCells(xadj.size()-1);
114  label nTotalConnections = returnReduce(adjncy.size(), sumOp<label>());
115 
116  // Send all to master. Use scheduled to save some storage.
117  if (Pstream::master())
118  {
119  Field<label> allAdjncy(nTotalConnections);
120  Field<label> allXadj(globalCells.size()+1);
121  scalarField allWeights(globalCells.size());
122 
123  // Insert my own
124  label nTotalCells = 0;
125  forAll(cWeights, celli)
126  {
127  allXadj[nTotalCells] = xadj[celli];
128  allWeights[nTotalCells++] = cWeights[celli];
129  }
130  nTotalConnections = 0;
131  forAll(adjncy, i)
132  {
133  allAdjncy[nTotalConnections++] = adjncy[i];
134  }
135 
136  for (int slave=1; slave<Pstream::nProcs(); slave++)
137  {
138  IPstream fromSlave(Pstream::commsTypes::scheduled, slave);
139  Field<label> nbrAdjncy(fromSlave);
140  Field<label> nbrXadj(fromSlave);
141  scalarField nbrWeights(fromSlave);
142 
143  // Append.
144  // label procStart = nTotalCells;
145  forAll(nbrXadj, celli)
146  {
147  allXadj[nTotalCells] = nTotalConnections+nbrXadj[celli];
148  allWeights[nTotalCells++] = nbrWeights[celli];
149  }
150  // No need to renumber xadj since already global.
151  forAll(nbrAdjncy, i)
152  {
153  allAdjncy[nTotalConnections++] = nbrAdjncy[i];
154  }
155  }
156  allXadj[nTotalCells] = nTotalConnections;
157 
158 
159  Field<label> allFinalDecomp;
160  decomposeOneProc
161  (
162  meshPath,
163  allAdjncy,
164  allXadj,
165  allWeights,
166  allFinalDecomp
167  );
168 
169 
170  // Send allFinalDecomp back
171  for (int slave=1; slave<Pstream::nProcs(); slave++)
172  {
173  OPstream toSlave(Pstream::commsTypes::scheduled, slave);
174  toSlave << SubField<label>
175  (
176  allFinalDecomp,
177  globalCells.localSize(slave),
178  globalCells.offset(slave)
179  );
180  }
181  // Get my own part (always first)
182  finalDecomp = SubField<label>
183  (
184  allFinalDecomp,
185  globalCells.localSize()
186  );
187  }
188  else
189  {
190  // Send my part of the graph (already in global numbering)
191  {
192  OPstream toMaster
193  (
196  );
197  toMaster<< adjncy << SubField<label>(xadj, xadj.size()-1)
198  << cWeights;
199  }
200 
201  // Receive back decomposition
202  IPstream fromMaster
203  (
206  );
207  fromMaster >> finalDecomp;
208  }
209  }
210  return 0;
211 }
212 
213 
214 // Call scotch with options from dictionary.
215 Foam::label Foam::scotchDecomp::decomposeOneProc
216 (
217  const fileName& meshPath,
218  const List<label>& adjncy,
219  const List<label>& xadj,
220  const scalarField& cWeights,
221 
222  List<label>& finalDecomp
223 )
224 {
225  // Dump graph
226  if (decompositionDict_.found("scotchCoeffs"))
227  {
228  const dictionary& scotchCoeffs =
229  decompositionDict_.subDict("scotchCoeffs");
230 
231  if (scotchCoeffs.lookupOrDefault("writeGraph", false))
232  {
233  OFstream str(meshPath + ".grf");
234 
235  Info<< "Dumping Scotch graph file to " << str.name() << endl
236  << "Use this in combination with gpart." << endl;
237 
238  label version = 0;
239  str << version << nl;
240  // Number of vertices
241  str << xadj.size()-1 << ' ' << adjncy.size() << nl;
242  // Numbering starts from 0
243  label baseval = 0;
244  // Has weights?
245  label hasEdgeWeights = 0;
246  label hasVertexWeights = 0;
247  label numericflag = 10*hasEdgeWeights+hasVertexWeights;
248  str << baseval << ' ' << numericflag << nl;
249  for (label celli = 0; celli < xadj.size()-1; celli++)
250  {
251  label start = xadj[celli];
252  label end = xadj[celli+1];
253  str << end-start;
254 
255  for (label i = start; i < end; i++)
256  {
257  str << ' ' << adjncy[i];
258  }
259  str << nl;
260  }
261  }
262  }
263 
264  // Reset the seed of the pseudo-random generator used by the graph
265  // partitioning routines of the libScotch library. Two consecutive calls to
266  // the same libScotch partitioning routines, and separated by a call to
267  // SCOTCH randomReset, will always yield the same results, as if the
268  // equivalent standalone Scotch programs were used twice, independently,
269  SCOTCH_randomReset();
270 
271  // Strategy
272  // ~~~~~~~~
273 
274  // Default.
275  SCOTCH_Strat stradat;
276  check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
277 
278  if (decompositionDict_.found("scotchCoeffs"))
279  {
280  const dictionary& scotchCoeffs =
281  decompositionDict_.subDict("scotchCoeffs");
282 
283  string strategy;
284  if (scotchCoeffs.readIfPresent("strategy", strategy))
285  {
286  if (debug)
287  {
288  Info<< "scotchDecomp : Using strategy " << strategy << endl;
289  }
290  SCOTCH_stratGraphMap(&stradat, strategy.c_str());
291  }
292  }
293 
294 
295  // Graph
296  // ~~~~~
297 
298  List<label> velotab;
299 
300 
301  // Check for externally provided cellweights and if so initialise weights
302  // Note: min, not gMin since routine runs on master only.
303  scalar minWeights = min(cWeights);
304  if (!cWeights.empty())
305  {
306  if (minWeights <= 0)
307  {
309  << "Illegal minimum weight " << minWeights
310  << endl;
311  }
312 
313  if (cWeights.size() != xadj.size()-1)
314  {
316  << "Number of cell weights " << cWeights.size()
317  << " does not equal number of cells " << xadj.size()-1
318  << exit(FatalError);
319  }
320 
321  scalar velotabSum = sum(cWeights)/minWeights;
322 
323  scalar rangeScale(1.0);
324 
325  if (velotabSum > scalar(labelMax - 1))
326  {
327  // 0.9 factor of safety to avoid floating point round-off in
328  // rangeScale tipping the subsequent sum over the integer limit.
329  rangeScale = 0.9*scalar(labelMax - 1)/velotabSum;
330 
332  << "Sum of weights has overflowed integer: " << velotabSum
333  << ", compressing weight scale by a factor of " << rangeScale
334  << endl;
335  }
336 
337  // Convert to integers.
338  velotab.setSize(cWeights.size());
339 
340  forAll(velotab, i)
341  {
342  velotab[i] = int((cWeights[i]/minWeights - 1)*rangeScale) + 1;
343  }
344  }
345 
346 
347 
348  SCOTCH_Graph grafdat;
349  check(SCOTCH_graphInit(&grafdat), "SCOTCH_graphInit");
350  check
351  (
352  SCOTCH_graphBuild
353  (
354  &grafdat,
355  0, // baseval, c-style numbering
356  xadj.size()-1, // vertnbr, nCells
357  xadj.begin(), // verttab, start index per cell into adjncy
358  &xadj[1], // vendtab, end index ,,
359  velotab.begin(), // velotab, vertex weights
360  nullptr, // vlbltab
361  adjncy.size(), // edgenbr, number of arcs
362  adjncy.begin(), // edgetab
363  nullptr // edlotab, edge weights
364  ),
365  "SCOTCH_graphBuild"
366  );
367  check(SCOTCH_graphCheck(&grafdat), "SCOTCH_graphCheck");
368 
369 
370  // Architecture
371  // ~~~~~~~~~~~~
372  // (fully connected network topology since using switch)
373 
374  SCOTCH_Arch archdat;
375  check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
376 
377  List<label> processorWeights;
378  if (decompositionDict_.found("scotchCoeffs"))
379  {
380  const dictionary& scotchCoeffs =
381  decompositionDict_.subDict("scotchCoeffs");
382 
383  scotchCoeffs.readIfPresent("processorWeights", processorWeights);
384  }
385  if (processorWeights.size())
386  {
387  if (debug)
388  {
389  Info<< "scotchDecomp : Using processor weights " << processorWeights
390  << endl;
391  }
392  check
393  (
394  SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
395  "SCOTCH_archCmpltw"
396  );
397  }
398  else
399  {
400  check
401  (
402  SCOTCH_archCmplt(&archdat, nProcessors_),
403  "SCOTCH_archCmplt"
404  );
405  }
406 
407  // Hack:switch off fpu error trapping
408  #ifdef FE_NOMASK_ENV
409  int oldExcepts = fedisableexcept
410  (
411  FE_DIVBYZERO
412  | FE_INVALID
413  | FE_OVERFLOW
414  );
415  #endif
416 
417  finalDecomp.setSize(xadj.size()-1);
418  finalDecomp = 0;
419  check
420  (
421  SCOTCH_graphMap
422  (
423  &grafdat,
424  &archdat,
425  &stradat, // const SCOTCH_Strat *
426  finalDecomp.begin() // parttab
427  ),
428  "SCOTCH_graphMap"
429  );
430 
431  #ifdef FE_NOMASK_ENV
432  feenableexcept(oldExcepts);
433  #endif
434 
435  // Release storage for graph
436  SCOTCH_graphExit(&grafdat);
437 
438  // Release storage for strategy
439  SCOTCH_stratExit(&stradat);
440 
441  // Release storage for network topology
442  SCOTCH_archExit(&archdat);
443 
444  return 0;
445 }
446 
447 
448 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
449 
450 Foam::scotchDecomp::scotchDecomp(const dictionary& decompositionDict)
451 :
452  decompositionMethod(decompositionDict)
453 {}
454 
455 
456 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
457 
458 Foam::labelList Foam::scotchDecomp::decompose
459 (
460  const polyMesh& mesh,
461  const pointField& points,
462  const scalarField& pointWeights
463 )
464 {
465  if (points.size() != mesh.nCells())
466  {
468  << "Can use this decomposition method only for the whole mesh"
469  << endl
470  << "and supply one coordinate (cellCentre) for every cell." << endl
471  << "The number of coordinates " << points.size() << endl
472  << "The number of cells in the mesh " << mesh.nCells()
473  << exit(FatalError);
474  }
475 
476  // Calculate local or global (if Pstream::parRun()) connectivity
477  CompactListList<label> cellCells;
479  (
480  mesh,
481  identity(mesh.nCells()),
482  mesh.nCells(),
483  true,
484  cellCells
485  );
486 
487  // Decompose using default weights
488  List<label> finalDecomp;
489  decompose
490  (
491  mesh.time().path()/mesh.name(),
492  cellCells.m(),
493  cellCells.offsets(),
494  pointWeights,
495  finalDecomp
496  );
497 
498  // Copy back to labelList
499  labelList decomp(finalDecomp.size());
500  forAll(decomp, i)
501  {
502  decomp[i] = finalDecomp[i];
503  }
504  return decomp;
505 }
506 
507 
508 Foam::labelList Foam::scotchDecomp::decompose
509 (
510  const polyMesh& mesh,
511  const labelList& agglom,
512  const pointField& agglomPoints,
513  const scalarField& pointWeights
514 )
515 {
516  if (agglom.size() != mesh.nCells())
517  {
519  << "Size of cell-to-coarse map " << agglom.size()
520  << " differs from number of cells in mesh " << mesh.nCells()
521  << exit(FatalError);
522  }
523 
524  // Calculate local or global (if Pstream::parRun()) connectivity
525  CompactListList<label> cellCells;
527  (
528  mesh,
529  agglom,
530  agglomPoints.size(),
531  true,
532  cellCells
533  );
534 
535  // Decompose using weights
536  List<label> finalDecomp;
537  decompose
538  (
539  mesh.time().path()/mesh.name(),
540  cellCells.m(),
541  cellCells.offsets(),
542  pointWeights,
543  finalDecomp
544  );
545 
546  // Rework back into decomposition for original mesh_
547  labelList fineDistribution(agglom.size());
548 
549  forAll(fineDistribution, i)
550  {
551  fineDistribution[i] = finalDecomp[agglom[i]];
552  }
553 
554  return fineDistribution;
555 }
556 
557 
558 Foam::labelList Foam::scotchDecomp::decompose
559 (
560  const labelListList& globalCellCells,
561  const pointField& cellCentres,
562  const scalarField& cWeights
563 )
564 {
565  if (cellCentres.size() != globalCellCells.size())
566  {
568  << "Inconsistent number of cells (" << globalCellCells.size()
569  << ") and number of cell centres (" << cellCentres.size()
570  << ")." << exit(FatalError);
571  }
572 
573 
574  // Make Metis CSR (Compressed Storage Format) storage
575  // adjncy : contains neighbours (= edges in graph)
576  // xadj(celli) : start of information in adjncy for celli
577 
578  CompactListList<label> cellCells(globalCellCells);
579 
580  // Decompose using weights
581  List<label> finalDecomp;
582  decompose
583  (
584  "scotch",
585  cellCells.m(),
586  cellCells.offsets(),
587  cWeights,
588  finalDecomp
589  );
590 
591  // Copy back to labelList
592  labelList decomp(finalDecomp.size());
593  forAll(decomp, i)
594  {
595  decomp[i] = finalDecomp[i];
596  }
597  return decomp;
598 }
599 
600 
601 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
static int masterNo()
Process index of the master.
Definition: UPstream.H:417
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:306
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
Holds information (coordinate and normal) regarding nearest wall point.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
static void calcCellCells(const polyMesh &mesh, const labelList &agglom, const label nLocalCoarse, const bool global, CompactListList< label > &cellCells)
Helper: determine (local or global) cellCells from mesh.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
static const label labelMax
Definition: label.H:62
List< label > labelList
A List of labels.
Definition: labelList.H:56
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
static const char nl
Definition: Ostream.H:260
defineTypeNameAndDebug(combustionModel, 0)
scotchDecomp(const dictionary &decompositionDict)
Construct given the decomposition dictionary and mesh.
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
#define WarningInFunction
Report a warning using Foam::Warning.
messageStream Info
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Namespace for OpenFOAM.