ptscotchDecomp.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2015-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "ptscotchDecomp.H"
31 #include "floatScalar.H"
32 #include "Time.H"
33 #include "PrecisionAdaptor.H"
34 #include "OFstream.H"
35 #include <limits>
36 
37 // Avoid too many warnings from mpi.h
38 #pragma GCC diagnostic ignored "-Wold-style-cast"
39 
40 #include <cstdio>
41 #include <mpi.h>
42 #include "ptscotch.h"
43 
44 // Hack: scotch generates floating point errors so need to switch off error
45 // trapping!
46 #ifdef __GLIBC__
47  #ifndef _GNU_SOURCE
48  #define _GNU_SOURCE
49  #endif
50  #include <fenv.h>
51 #endif
52 
53 // Error if we attempt narrowing
54 static_assert
55 (
56  sizeof(Foam::label) <= sizeof(SCOTCH_Num),
57  "SCOTCH_Num is too small for Foam::label, check your scotch headers"
58 );
59 
60 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
61 
62 namespace Foam
63 {
64  defineTypeNameAndDebug(ptscotchDecomp, 0);
66  (
67  decompositionMethod,
68  ptscotchDecomp,
69  dictionary
70  );
71 }
72 
73 
74 // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
75 
76 namespace Foam
77 {
78 
79 // Check and print error message
80 static inline void check(const int retVal, const char* what)
81 {
82  if (retVal)
83  {
85  << "Call to scotch routine " << what
86  << " failed (" << retVal << ")\n"
87  << exit(FatalError);
88  }
89 }
90 
91 // The mesh-relative graph path/name (without extension)
93 {
94  return mesh.time().path()/mesh.name();
95 }
96 
97 
98 } // End namespace Foam
99 
100 
101 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
102 
103 Foam::label Foam::ptscotchDecomp::decompose
104 (
105  const labelList& adjncy,
106  const labelList& xadj,
107  const List<scalar>& cWeights,
108  labelList& decomp
109 ) const
110 {
111  const SCOTCH_Num numCells = max(0, (xadj.size()-1));
112 
113  // Addressing
114  ConstPrecisionAdaptor<SCOTCH_Num, label, List> adjncy_param(adjncy);
115  ConstPrecisionAdaptor<SCOTCH_Num, label, List> xadj_param(xadj);
116 
117  // Output: cell -> processor addressing
118  decomp.resize(numCells);
119  decomp = 0;
120  PrecisionAdaptor<SCOTCH_Num, label, List> decomp_param(decomp, false);
121 
122  // Avoid potential nullptr issues with zero-sized arrays
123  labelList adjncy_dummy, xadj_dummy, decomp_dummy;
124  if (!numCells)
125  {
126  adjncy_dummy.resize(1, 0);
127  adjncy_param.set(adjncy_dummy);
128 
129  xadj_dummy.resize(2, 0);
130  xadj_param.set(xadj_dummy);
131 
132  decomp_dummy.resize(1, 0);
133  decomp_param.clear(); // Avoid propagating spurious values
134  decomp_param.set(decomp_dummy);
135  }
136 
137 
138  if (debug & 2)
139  {
140  Pout<< "ptscotchDecomp : " << numCells << " cells" << endl;
141  }
142 
143  // Dump graph
144  if (coeffsDict_.getOrDefault("writeGraph", false))
145  {
146  OFstream str
147  (
148  graphPath_ + "_" + Foam::name(Pstream::myProcNo()) + ".dgr"
149  );
150 
151  Pout<< "Dumping Scotch graph file to " << str.name() << endl
152  << "Use this in combination with dgpart." << endl;
153 
154  const label numConnect = adjncy.size();
155  const label nTotCells = returnReduce(numCells, sumOp<label>());
156  const label nTotConnect = returnReduce(numConnect, sumOp<label>());
157 
158  // Version 2 = Distributed graph file (.dgrf)
159  str << "2" << nl;
160 
161  // Number of files (procglbnbr), my file number (procloc)
162  str << Pstream::nProcs() << ' ' << Pstream::myProcNo() << nl;
163 
164  // Total number of vertices (vertglbnbr),
165  // Total number of connections (edgeglbnbr)
166  str << nTotCells << ' ' << nTotConnect << nl;
167 
168  // Local number of vertices (vertlocnbr),
169  // Local number of connections (edgelocnbr)
170  str << numCells << ' ' << numConnect << nl;
171 
172  // Numbering starts from 0
173  // 100*hasVertlabels+10*hasEdgeWeights+1*hasVertWeights
174  str << "0 000" << nl;
175 
176  for (label celli = 0; celli < numCells; ++celli)
177  {
178  const label beg = xadj[celli];
179  const label end = xadj[celli+1];
180 
181  str << (end-beg); // size
182 
183  for (label i = beg; i < end; ++i)
184  {
185  str << ' ' << adjncy[i];
186  }
187  str << nl;
188  }
189  }
190 
191 
192  // Make repeatable
193  SCOTCH_randomReset();
194 
195  // Strategy
196  // ~~~~~~~~
197 
198  // Default.
199  SCOTCH_Strat stradat;
200  check
201  (
202  SCOTCH_stratInit(&stradat),
203  "SCOTCH_stratInit"
204  );
205 
206  string strategy;
207  if (coeffsDict_.readIfPresent("strategy", strategy))
208  {
209  DebugInfo
210  << "ptscotchDecomp : Using strategy " << strategy << endl;
211 
212  SCOTCH_stratDgraphMap(&stradat, strategy.c_str());
213  //fprintf(stdout, "S\tStrat=");
214  //SCOTCH_stratSave(&stradat, stdout);
215  //fprintf(stdout, "\n");
216  }
217 
218 
219  // Graph
220  // ~~~~~
221 
222  // Check for externally provided cellweights and if so initialise weights
223 
224  bool hasWeights = returnReduce(!cWeights.empty(), orOp<bool>());
225 
226  const scalar minWeights = hasWeights ? gMin(cWeights) : scalar(1);
227 
228  if (minWeights <= 0)
229  {
230  hasWeights = false;
232  << "Illegal minimum weight " << minWeights
233  << " ... ignoring"
234  << endl;
235  }
236  else if (hasWeights && (cWeights.size() != numCells))
237  {
239  << "Number of cell weights " << cWeights.size()
240  << " does not equal number of cells " << numCells
241  << exit(FatalError);
242  }
243 
244 
245  List<SCOTCH_Num> velotab;
246 
247  if (hasWeights)
248  {
249  scalar rangeScale(1);
250 
251  const scalar velotabSum = gSum(cWeights)/minWeights;
252 
253  const scalar upperRange = static_cast<scalar>
254  (
256  );
257 
258  if (velotabSum > upperRange)
259  {
260  // 0.9 factor of safety to avoid floating point round-off in
261  // rangeScale tipping the subsequent sum over the integer limit.
262  rangeScale = 0.9*upperRange/velotabSum;
263 
265  << "Sum of weights overflows SCOTCH_Num: " << velotabSum
266  << ", compressing by factor " << rangeScale << endl;
267  }
268 
269  if (cWeights.size())
270  {
271  // Convert to integers.
272  velotab.resize(cWeights.size());
273 
274  forAll(velotab, i)
275  {
276  velotab[i] = static_cast<SCOTCH_Num>
277  (
278  ((cWeights[i]/minWeights - 1)*rangeScale) + 1
279  );
280  }
281  }
282  else
283  {
284  // Locally zero cells but not globally.
285  // Provide some size to avoid null pointer.
286  velotab.resize(1, 1);
287  }
288  }
289 
290 
291  //
292  // Decomposition graph
293  //
294 
295  if (debug & 2)
296  {
297  Pout<< "SCOTCH_dgraphInit" << endl;
298  }
299  SCOTCH_Dgraph grafdat;
300  check
301  (
302  SCOTCH_dgraphInit(&grafdat, MPI_COMM_WORLD),
303  "SCOTCH_dgraphInit"
304  );
305 
306  if (debug & 2)
307  {
308  Pout<< "SCOTCH_dgraphBuild with:" << nl
309  << "numCells : " << numCells << nl
310  << "xadj : " << name(xadj_param().cdata()) << nl
311  << "velotab : " << name(velotab.cdata()) << nl
312  << "adjncySize : " << adjncy_param().size() << nl
313  << "adjncy : " << name(adjncy_param().cdata()) << nl
314  << endl;
315  }
316 
317  check
318  (
319  SCOTCH_dgraphBuild
320  (
321  &grafdat, // Graph to build
322  0, // Base for indexing (C-style)
323 
324  numCells, // vertlocnbr [== nCells]
325  numCells, // vertlocmax
326 
327  xadj_param.constCast().data(),
328  // vertloctab, start index per cell into
329  // adjncy
330  (xadj_param.constCast().data()+1),
331  // vendloctab, end index ,,
332 
333  velotab.data(), // veloloctab, vtx weights
334  nullptr, // vlblloctab
335 
336  adjncy.size(), // edgelocnbr, number of arcs
337  adjncy.size(), // edgelocsiz
338  adjncy_param.constCast().data(), // edgeloctab
339  nullptr, // edgegsttab
340  nullptr // edlotab, edge weights
341  ),
342  "SCOTCH_dgraphBuild"
343  );
344 
345  if (debug & 2)
346  {
347  Pout<< "SCOTCH_dgraphCheck" << endl;
348  }
349  check
350  (
351  SCOTCH_dgraphCheck(&grafdat),
352  "SCOTCH_dgraphCheck"
353  );
354 
355 
356  // Architecture
357  // ~~~~~~~~~~~~
358  // (fully connected network topology since using switch)
359 
360  if (debug & 2)
361  {
362  Pout<< "SCOTCH_archInit" << endl;
363  }
364  SCOTCH_Arch archdat;
365  check
366  (
367  SCOTCH_archInit(&archdat),
368  "SCOTCH_archInit"
369  );
370 
371  List<SCOTCH_Num> procWeights;
372  if
373  (
374  coeffsDict_.readIfPresent("processorWeights", procWeights)
375  && !procWeights.empty()
376  )
377  {
378  if (procWeights.size() != nDomains_)
379  {
380  FatalIOErrorInFunction(coeffsDict_)
381  << "processorWeights (" << procWeights.size()
382  << ") != number of domains (" << nDomains_ << ")" << nl
383  << exit(FatalIOError);
384  }
385 
386  DebugInfo
387  << "ptscotchDecomp : Using procesor weights "
388  << procWeights << endl;
389 
390  check
391  (
392  SCOTCH_archCmpltw(&archdat, nDomains_, procWeights.cdata()),
393  "SCOTCH_archCmpltw"
394  );
395  }
396  else
397  {
398  if (debug & 2)
399  {
400  Pout<< "SCOTCH_archCmplt" << endl;
401  }
402  check
403  (
404  SCOTCH_archCmplt(&archdat, nDomains_),
405  "SCOTCH_archCmplt"
406  );
407  }
408 
409 
410  // Hack:switch off fpu error trapping
411  #ifdef FE_NOMASK_ENV
412  int oldExcepts = fedisableexcept
413  (
414  FE_DIVBYZERO
415  | FE_INVALID
416  | FE_OVERFLOW
417  );
418  #endif
419 
420  if (debug & 2)
421  {
422  Pout<< "SCOTCH_dgraphMap" << endl;
423  }
424  check
425  (
426  SCOTCH_dgraphMap
427  (
428  &grafdat,
429  &archdat,
430  &stradat, // const SCOTCH_Strat *
431  decomp_param.ref().data() // parttab
432  ),
433  "SCOTCH_graphMap"
434  );
435 
436  #ifdef FE_NOMASK_ENV
437  feenableexcept(oldExcepts);
438  #endif
439 
440  //check
441  //(
442  // SCOTCH_dgraphPart
443  // (
444  // &grafdat,
445  // nDomains_, // partnbr
446  // &stradat, // const SCOTCH_Strat *
447  // decomp_param.ref().data() // parttab
448  // ),
449  // "SCOTCH_graphPart"
450  //);
451 
452  if (debug & 2)
453  {
454  Pout<< "SCOTCH_dgraphExit" << endl;
455  }
456 
457  SCOTCH_dgraphExit(&grafdat); // Release storage for graph
458  SCOTCH_stratExit(&stradat); // Release storage for strategy
459  SCOTCH_archExit(&archdat); // Release storage for network topology
460 
461  return 0;
462 }
463 
464 
465 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
466 
467 Foam::ptscotchDecomp::ptscotchDecomp
468 (
469  const dictionary& decompDict,
470  const word& regionName
471 )
472 :
473  decompositionMethod(decompDict, regionName),
474  coeffsDict_(findCoeffsDict("scotchCoeffs", selectionType::NULL_DICT))
475 {}
476 
477 
478 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
479 
480 Foam::labelList Foam::ptscotchDecomp::decompose
481 (
482  const polyMesh& mesh,
483  const pointField& points,
484  const scalarField& pointWeights
485 ) const
486 {
487  // Where to write graph
488  graphPath_ = getGraphPathBase(mesh);
489 
490  if (points.size() != mesh.nCells())
491  {
493  << "Can only use this decomposition method for entire mesh" << nl
494  << "and supply one coordinate (cellCentre) for every cell." << nl
495  << "The number of coordinates " << points.size() << nl
496  << "The number of cells in the mesh " << mesh.nCells() << nl
497  << exit(FatalError);
498  }
499 
500 
501  // Make Metis CSR (Compressed Storage Format) storage
502  // adjncy : contains neighbours (= edges in graph)
503  // xadj(celli) : start of information in adjncy for celli
504 
505  CompactListList<label> cellCells;
507  (
508  mesh,
509  identity(mesh.nCells()),
510  mesh.nCells(),
511  true,
512  cellCells
513  );
514 
515  // Decompose using default weights
516  labelList decomp;
517  decompose
518  (
519  cellCells.m(),
520  cellCells.offsets(),
521  pointWeights,
522  decomp
523  );
524 
525  return decomp;
526 }
527 
528 
529 Foam::labelList Foam::ptscotchDecomp::decompose
530 (
531  const polyMesh& mesh,
532  const labelList& agglom,
533  const pointField& agglomPoints,
534  const scalarField& pointWeights
535 ) const
536 {
537  // Where to write graph
538  graphPath_ = getGraphPathBase(mesh);
539 
540  if (agglom.size() != mesh.nCells())
541  {
543  << "Size of cell-to-coarse map " << agglom.size()
544  << " differs from number of cells in mesh " << mesh.nCells()
545  << exit(FatalError);
546  }
547 
548 
549  // Make Metis CSR (Compressed Storage Format) storage
550  // adjncy : contains neighbours (= edges in graph)
551  // xadj(celli) : start of information in adjncy for celli
552  CompactListList<label> cellCells;
553  calcCellCells
554  (
555  mesh,
556  agglom,
557  agglomPoints.size(),
558  true,
559  cellCells
560  );
561 
562  // Decompose using weights
563  labelList decomp;
564  decompose
565  (
566  cellCells.m(),
567  cellCells.offsets(),
568  pointWeights,
569  decomp
570  );
571 
572  // Rework back into decomposition for original mesh
573  labelList fineDistribution(agglom.size());
574 
575  forAll(fineDistribution, i)
576  {
577  fineDistribution[i] = decomp[agglom[i]];
578  }
579 
580  return fineDistribution;
581 }
582 
583 
584 Foam::labelList Foam::ptscotchDecomp::decompose
585 (
586  const labelListList& globalCellCells,
587  const pointField& cellCentres,
588  const scalarField& cWeights
589 ) const
590 {
591  // Where to write graph
592  graphPath_ = "ptscotch";
593 
594  if (cellCentres.size() != globalCellCells.size())
595  {
597  << "Inconsistent number of cells (" << globalCellCells.size()
598  << ") and number of cell centres (" << cellCentres.size()
599  << ")." << exit(FatalError);
600  }
601 
602 
603  // Make Metis CSR (Compressed Storage Format) storage
604  // adjncy : contains neighbours (= edges in graph)
605  // xadj(celli) : start of information in adjncy for celli
606 
607  CompactListList<label> cellCells(globalCellCells);
608 
609  // Decompose using weights
610  labelList decomp;
611  decompose
612  (
613  cellCells.m(),
614  cellCells.offsets(),
615  cWeights,
616  decomp
617  );
618 
619  return decomp;
620 }
621 
622 
623 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
Foam::decompositionMethod::nDomains_
label nDomains_
Number of domains for the decomposition.
Definition: decompositionMethod.H:95
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::HashTable::size
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::decompositionMethod::calcCellCells
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.
Definition: decompositionMethod.C:468
Foam::fileName
A class for handling file names.
Definition: fileName.H:73
Foam::returnReduce
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Definition: PstreamReduceOps.H:94
Foam::List::resize
void resize(const label len)
Adjust allocated size of list.
Definition: ListI.H:139
PrecisionAdaptor.H
Foam::FatalIOError
IOerror FatalIOError
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::gSum
Type gSum(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:594
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
OFstream.H
regionName
Foam::word regionName
Definition: createNamedDynamicFvMesh.H:1
Foam::primitiveMesh::nCells
label nCells() const noexcept
Number of mesh cells.
Definition: primitiveMeshI.H:96
Foam::check
static void check(const int retVal, const char *what)
Definition: ptscotchDecomp.C:80
Foam::OSstream::name
virtual const fileName & name() const
Get the name of the stream.
Definition: OSstream.H:107
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::FatalError
error FatalError
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:121
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Time.H
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
DebugInfo
#define DebugInfo
Report an information message using Foam::Info.
Definition: messageStream.H:382
Foam::UPstream::myProcNo
static int myProcNo(const label communicator=worldComm)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:463
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::Time::path
fileName path() const
Return path.
Definition: Time.H:358
Foam::List
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: BitOps.H:63
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::identity
labelList identity(const label len, label start=0)
Create identity map of the given length with (map[i] == i)
Definition: labelList.C:38
floatScalar.H
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Foam::fvMesh::time
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:280
ptscotchDecomp.H
FatalIOErrorInFunction
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:473
Foam::gMin
Type gMin(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:593
Foam::getGraphPathBase
static Foam::fileName getGraphPathBase(const polyMesh &mesh)
Definition: ptscotchDecomp.C:92
Foam::dictionary::getOrDefault
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:148
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::UPstream::nProcs
static label nProcs(const label communicator=worldComm)
Number of processes in parallel run, and 1 for serial run.
Definition: UPstream.H:445
Foam::dictionary::readIfPresent
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:405
Foam::fvMesh::name
const word & name() const
Return reference to name.
Definition: fvMesh.H:300