scotchDecomp.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-2017 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  dictionary
61  );
62 }
63 
64 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
65 
66 void Foam::scotchDecomp::check(const int retVal, const char* str)
67 {
68  if (retVal)
69  {
71  << "Call to scotch routine " << str << " failed."
72  << exit(FatalError);
73  }
74 }
75 
76 
77 Foam::label Foam::scotchDecomp::decompose
78 (
79  const fileName& meshPath,
80  const List<label>& adjncy,
81  const List<label>& xadj,
82  const scalarField& cWeights,
83 
84  List<label>& finalDecomp
85 )
86 {
87  if (!Pstream::parRun())
88  {
89  decomposeOneProc
90  (
91  meshPath,
92  adjncy,
93  xadj,
94  cWeights,
95  finalDecomp
96  );
97  }
98  else
99  {
100  if (debug)
101  {
102  Info<< "scotchDecomp : running in parallel."
103  << " Decomposing all of graph on master processor." << endl;
104  }
105  globalIndex globalCells(xadj.size()-1);
106  label nTotalConnections = returnReduce(adjncy.size(), sumOp<label>());
107 
108  // Send all to master. Use scheduled to save some storage.
109  if (Pstream::master())
110  {
111  Field<label> allAdjncy(nTotalConnections);
112  Field<label> allXadj(globalCells.size()+1);
113  scalarField allWeights(globalCells.size());
114 
115  // Insert my own
116  label nTotalCells = 0;
117  forAll(cWeights, celli)
118  {
119  allXadj[nTotalCells] = xadj[celli];
120  allWeights[nTotalCells++] = cWeights[celli];
121  }
122  nTotalConnections = 0;
123  forAll(adjncy, i)
124  {
125  allAdjncy[nTotalConnections++] = adjncy[i];
126  }
127 
128  for (int slave=1; slave<Pstream::nProcs(); slave++)
129  {
130  IPstream fromSlave(Pstream::commsTypes::scheduled, slave);
131  Field<label> nbrAdjncy(fromSlave);
132  Field<label> nbrXadj(fromSlave);
133  scalarField nbrWeights(fromSlave);
134 
135  // Append.
136  //label procStart = nTotalCells;
137  forAll(nbrXadj, celli)
138  {
139  allXadj[nTotalCells] = nTotalConnections+nbrXadj[celli];
140  allWeights[nTotalCells++] = nbrWeights[celli];
141  }
142  // No need to renumber xadj since already global.
143  forAll(nbrAdjncy, i)
144  {
145  allAdjncy[nTotalConnections++] = nbrAdjncy[i];
146  }
147  }
148  allXadj[nTotalCells] = nTotalConnections;
149 
150 
151  Field<label> allFinalDecomp;
152  decomposeOneProc
153  (
154  meshPath,
155  allAdjncy,
156  allXadj,
157  allWeights,
158  allFinalDecomp
159  );
160 
161 
162  // Send allFinalDecomp back
163  for (int slave=1; slave<Pstream::nProcs(); slave++)
164  {
165  OPstream toSlave(Pstream::commsTypes::scheduled, slave);
166  toSlave << SubField<label>
167  (
168  allFinalDecomp,
169  globalCells.localSize(slave),
170  globalCells.offset(slave)
171  );
172  }
173  // Get my own part (always first)
174  finalDecomp = SubField<label>
175  (
176  allFinalDecomp,
177  globalCells.localSize()
178  );
179  }
180  else
181  {
182  // Send my part of the graph (already in global numbering)
183  {
184  OPstream toMaster
185  (
188  );
189  toMaster<< adjncy << SubField<label>(xadj, xadj.size()-1)
190  << cWeights;
191  }
192 
193  // Receive back decomposition
194  IPstream fromMaster
195  (
198  );
199  fromMaster >> finalDecomp;
200  }
201  }
202  return 0;
203 }
204 
205 
206 // Call scotch with options from dictionary.
207 Foam::label Foam::scotchDecomp::decomposeOneProc
208 (
209  const fileName& meshPath,
210  const List<label>& adjncy,
211  const List<label>& xadj,
212  const scalarField& cWeights,
213 
214  List<label>& finalDecomp
215 )
216 {
217  // Dump graph
218  if (decompositionDict_.found("scotchCoeffs"))
219  {
220  const dictionary& scotchCoeffs =
221  decompositionDict_.subDict("scotchCoeffs");
222 
223  if (scotchCoeffs.lookupOrDefault("writeGraph", false))
224  {
225  OFstream str(meshPath + ".grf");
226 
227  Info<< "Dumping Scotch graph file to " << str.name() << endl
228  << "Use this in combination with gpart." << endl;
229 
230  label version = 0;
231  str << version << nl;
232  // Numer of vertices
233  str << xadj.size()-1 << ' ' << adjncy.size() << nl;
234  // Numbering starts from 0
235  label baseval = 0;
236  // Has weights?
237  label hasEdgeWeights = 0;
238  label hasVertexWeights = 0;
239  label numericflag = 10*hasEdgeWeights+hasVertexWeights;
240  str << baseval << ' ' << numericflag << nl;
241  for (label celli = 0; celli < xadj.size()-1; celli++)
242  {
243  label start = xadj[celli];
244  label end = xadj[celli+1];
245  str << end-start;
246 
247  for (label i = start; i < end; i++)
248  {
249  str << ' ' << adjncy[i];
250  }
251  str << nl;
252  }
253  }
254  }
255 
256 
257  // Strategy
258  // ~~~~~~~~
259 
260  // Default.
261  SCOTCH_Strat stradat;
262  check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
263 
264  if (decompositionDict_.found("scotchCoeffs"))
265  {
266  const dictionary& scotchCoeffs =
267  decompositionDict_.subDict("scotchCoeffs");
268 
269  string strategy;
270  if (scotchCoeffs.readIfPresent("strategy", strategy))
271  {
272  if (debug)
273  {
274  Info<< "scotchDecomp : Using strategy " << strategy << endl;
275  }
276  SCOTCH_stratGraphMap(&stradat, strategy.c_str());
277  //fprintf(stdout, "S\tStrat=");
278  //SCOTCH_stratSave(&stradat, stdout);
279  //fprintf(stdout, "\n");
280  }
281  }
282 
283 
284  // Graph
285  // ~~~~~
286 
287  List<label> velotab;
288 
289 
290  // Check for externally provided cellweights and if so initialise weights
291  // Note: min, not gMin since routine runs on master only.
292  scalar minWeights = min(cWeights);
293  if (!cWeights.empty())
294  {
295  if (minWeights <= 0)
296  {
298  << "Illegal minimum weight " << minWeights
299  << endl;
300  }
301 
302  if (cWeights.size() != xadj.size()-1)
303  {
305  << "Number of cell weights " << cWeights.size()
306  << " does not equal number of cells " << xadj.size()-1
307  << exit(FatalError);
308  }
309 
310  scalar velotabSum = sum(cWeights)/minWeights;
311 
312  scalar rangeScale(1.0);
313 
314  if (velotabSum > scalar(labelMax - 1))
315  {
316  // 0.9 factor of safety to avoid floating point round-off in
317  // rangeScale tipping the subsequent sum over the integer limit.
318  rangeScale = 0.9*scalar(labelMax - 1)/velotabSum;
319 
321  << "Sum of weights has overflowed integer: " << velotabSum
322  << ", compressing weight scale by a factor of " << rangeScale
323  << endl;
324  }
325 
326  // Convert to integers.
327  velotab.setSize(cWeights.size());
328 
329  forAll(velotab, i)
330  {
331  velotab[i] = int((cWeights[i]/minWeights - 1)*rangeScale) + 1;
332  }
333  }
334 
335 
336 
337  SCOTCH_Graph grafdat;
338  check(SCOTCH_graphInit(&grafdat), "SCOTCH_graphInit");
339  check
340  (
341  SCOTCH_graphBuild
342  (
343  &grafdat,
344  0, // baseval, c-style numbering
345  xadj.size()-1, // vertnbr, nCells
346  xadj.begin(), // verttab, start index per cell into adjncy
347  &xadj[1], // vendtab, end index ,,
348  velotab.begin(), // velotab, vertex weights
349  nullptr, // vlbltab
350  adjncy.size(), // edgenbr, number of arcs
351  adjncy.begin(), // edgetab
352  nullptr // edlotab, edge weights
353  ),
354  "SCOTCH_graphBuild"
355  );
356  check(SCOTCH_graphCheck(&grafdat), "SCOTCH_graphCheck");
357 
358 
359  // Architecture
360  // ~~~~~~~~~~~~
361  // (fully connected network topology since using switch)
362 
363  SCOTCH_Arch archdat;
364  check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
365 
366  List<label> processorWeights;
367  if (decompositionDict_.found("scotchCoeffs"))
368  {
369  const dictionary& scotchCoeffs =
370  decompositionDict_.subDict("scotchCoeffs");
371 
372  scotchCoeffs.readIfPresent("processorWeights", processorWeights);
373  }
374  if (processorWeights.size())
375  {
376  if (debug)
377  {
378  Info<< "scotchDecomp : Using procesor weights " << processorWeights
379  << endl;
380  }
381  check
382  (
383  SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
384  "SCOTCH_archCmpltw"
385  );
386  }
387  else
388  {
389  check
390  (
391  SCOTCH_archCmplt(&archdat, nProcessors_),
392  "SCOTCH_archCmplt"
393  );
394 
395 
396  //- Hack to test clustering. Note that finalDecomp is non-compact
397  // numbers!
398  //
400  //check
401  //(
402  // SCOTCH_archVcmplt(&archdat),
403  // "SCOTCH_archVcmplt"
404  //);
405  //
407  //SCOTCH_Num straval = 0;
410  //
413  //SCOTCH_Num agglomSize = 3;
414  //
416  //check
417  //(
418  // SCOTCH_stratGraphClusterBuild
419  // (
420  // &stradat, // strategy to build
421  // straval, // strategy flags
422  // agglomSize, // cells per cluster
423  // 1.0, // weight?
424  // 0.01 // max load imbalance
425  // ),
426  // "SCOTCH_stratGraphClusterBuild"
427  //);
428  }
429 
430 
431  //SCOTCH_Mapping mapdat;
432  //SCOTCH_graphMapInit(&grafdat, &mapdat, &archdat, nullptr);
433  //SCOTCH_graphMapCompute(&grafdat, &mapdat, &stradat); /* Perform mapping */
434  //SCOTCH_graphMapExit(&grafdat, &mapdat);
435 
436 
437  // Hack:switch off fpu error trapping
438  #ifdef FE_NOMASK_ENV
439  int oldExcepts = fedisableexcept
440  (
441  FE_DIVBYZERO
442  | FE_INVALID
443  | FE_OVERFLOW
444  );
445  #endif
446 
447  finalDecomp.setSize(xadj.size()-1);
448  finalDecomp = 0;
449  check
450  (
451  SCOTCH_graphMap
452  (
453  &grafdat,
454  &archdat,
455  &stradat, // const SCOTCH_Strat *
456  finalDecomp.begin() // parttab
457  ),
458  "SCOTCH_graphMap"
459  );
460 
461  #ifdef FE_NOMASK_ENV
462  feenableexcept(oldExcepts);
463  #endif
464 
465 
466 
467  //finalDecomp.setSize(xadj.size()-1);
468  //check
469  //(
470  // SCOTCH_graphPart
471  // (
472  // &grafdat,
473  // nProcessors_, // partnbr
474  // &stradat, // const SCOTCH_Strat *
475  // finalDecomp.begin() // parttab
476  // ),
477  // "SCOTCH_graphPart"
478  //);
479 
480  // Release storage for graph
481  SCOTCH_graphExit(&grafdat);
482  // Release storage for strategy
483  SCOTCH_stratExit(&stradat);
484  // Release storage for network topology
485  SCOTCH_archExit(&archdat);
486 
487  return 0;
488 }
489 
490 
491 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
492 
493 Foam::scotchDecomp::scotchDecomp(const dictionary& decompositionDict)
494 :
495  decompositionMethod(decompositionDict)
496 {}
497 
498 
499 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
500 
501 Foam::labelList Foam::scotchDecomp::decompose
502 (
503  const polyMesh& mesh,
504  const pointField& points,
505  const scalarField& pointWeights
506 )
507 {
508  if (points.size() != mesh.nCells())
509  {
511  << "Can use this decomposition method only for the whole mesh"
512  << endl
513  << "and supply one coordinate (cellCentre) for every cell." << endl
514  << "The number of coordinates " << points.size() << endl
515  << "The number of cells in the mesh " << mesh.nCells()
516  << exit(FatalError);
517  }
518 
519  // Calculate local or global (if Pstream::parRun()) connectivity
520  CompactListList<label> cellCells;
522  (
523  mesh,
524  identity(mesh.nCells()),
525  mesh.nCells(),
526  true,
527  cellCells
528  );
529 
530  // Decompose using default weights
531  List<label> finalDecomp;
532  decompose
533  (
534  mesh.time().path()/mesh.name(),
535  cellCells.m(),
536  cellCells.offsets(),
537  pointWeights,
538  finalDecomp
539  );
540 
541  // Copy back to labelList
542  labelList decomp(finalDecomp.size());
543  forAll(decomp, i)
544  {
545  decomp[i] = finalDecomp[i];
546  }
547  return decomp;
548 }
549 
550 
551 Foam::labelList Foam::scotchDecomp::decompose
552 (
553  const polyMesh& mesh,
554  const labelList& agglom,
555  const pointField& agglomPoints,
556  const scalarField& pointWeights
557 )
558 {
559  if (agglom.size() != mesh.nCells())
560  {
562  << "Size of cell-to-coarse map " << agglom.size()
563  << " differs from number of cells in mesh " << mesh.nCells()
564  << exit(FatalError);
565  }
566 
567  // Calculate local or global (if Pstream::parRun()) connectivity
568  CompactListList<label> cellCells;
570  (
571  mesh,
572  agglom,
573  agglomPoints.size(),
574  true,
575  cellCells
576  );
577 
578  // Decompose using weights
579  List<label> finalDecomp;
580  decompose
581  (
582  mesh.time().path()/mesh.name(),
583  cellCells.m(),
584  cellCells.offsets(),
585  pointWeights,
586  finalDecomp
587  );
588 
589  // Rework back into decomposition for original mesh_
590  labelList fineDistribution(agglom.size());
591 
592  forAll(fineDistribution, i)
593  {
594  fineDistribution[i] = finalDecomp[agglom[i]];
595  }
596 
597  return fineDistribution;
598 }
599 
600 
601 Foam::labelList Foam::scotchDecomp::decompose
602 (
603  const labelListList& globalCellCells,
604  const pointField& cellCentres,
605  const scalarField& cWeights
606 )
607 {
608  if (cellCentres.size() != globalCellCells.size())
609  {
611  << "Inconsistent number of cells (" << globalCellCells.size()
612  << ") and number of cell centres (" << cellCentres.size()
613  << ")." << exit(FatalError);
614  }
615 
616 
617  // Make Metis CSR (Compressed Storage Format) storage
618  // adjncy : contains neighbours (= edges in graph)
619  // xadj(celli) : start of information in adjncy for celli
620 
621  CompactListList<label> cellCells(globalCellCells);
622 
623  // Decompose using weights
624  List<label> finalDecomp;
625  decompose
626  (
627  "scotch",
628  cellCells.m(),
629  cellCells.offsets(),
630  cWeights,
631  finalDecomp
632  );
633 
634  // Copy back to labelList
635  labelList decomp(finalDecomp.size());
636  forAll(decomp, i)
637  {
638  decomp[i] = finalDecomp[i];
639  }
640  return decomp;
641 }
642 
643 
644 // ************************************************************************* //
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
static int masterNo()
Process index of the master.
Definition: UPstream.H:406
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:412
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.
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:262
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:394
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:400
#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.