multiLevelDecomp.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) 2017-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 "multiLevelDecomp.H"
31 #include "IFstream.H"
32 #include "globalIndex.H"
33 #include "mapDistribute.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39  defineTypeNameAndDebug(multiLevelDecomp, 0);
41  (
42  decompositionMethod,
43  multiLevelDecomp,
44  dictionary
45  );
46 }
47 
48 
49 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
50 
51 void Foam::multiLevelDecomp::createMethodsDict()
52 {
53  methodsDict_.clear();
54 
55  word defaultMethod;
56  labelList domains;
57 
58  label nTotal = 0;
59  label nLevels = 0;
60 
61  // Found (non-recursive, no patterns) "method" and "domains" ?
62  // Allow as quick short-cut entry
63  if
64  (
65  // non-recursive, no patterns
66  coeffsDict_.readIfPresent("method", defaultMethod, keyType::LITERAL)
67  // non-recursive, no patterns
68  && coeffsDict_.readIfPresent("domains", domains, keyType::LITERAL)
69  )
70  {
71  // Short-cut version specified by method, domains only
72 
73  nTotal = (domains.empty() ? 0 : 1);
74 
75  for (const label n : domains)
76  {
77  nTotal *= n;
78  ++nLevels;
79  }
80 
81  if (nTotal == 1)
82  {
83  // Emit Warning
84  nTotal = nDomains();
85  nLevels = 1;
86 
87  domains.setSize(1);
88  domains[0] = nTotal;
89  }
90  else if (nTotal > 0 && nTotal < nDomains() && !(nDomains() % nTotal))
91  {
92  // nTotal < nDomains, but with an integral factor,
93  // which we insert as level 0
94  ++nLevels;
95 
96  labelList old(std::move(domains));
97 
98  domains.setSize(old.size()+1);
99 
100  domains[0] = nDomains() / nTotal;
101  forAll(old, i)
102  {
103  domains[i+1] = old[i];
104  }
105  nTotal *= domains[0];
106 
107  Info<<" inferred level0 with " << domains[0]
108  << " domains" << nl << nl;
109  }
110 
111  if (!nLevels || nTotal != nDomains())
112  {
114  << "Top level decomposition specifies " << nDomains()
115  << " domains which is not equal to the product of"
116  << " all sub domains " << nTotal
117  << exit(FatalError);
118  }
119 
120  // Create editable methods dictionaries
121  nLevels = 0;
122 
123  // Common coeffs dictionary
124  const dictionary& subMethodCoeffsDict
125  (
127  (
128  coeffsDict_,
129  defaultMethod + "Coeffs",
130  selectionType::NULL_DICT
131  )
132  );
133 
134  for (const label n : domains)
135  {
136  const word levelName("level" + Foam::name(nLevels++));
137 
138  entry* dictptr = methodsDict_.set(levelName, dictionary());
139 
140  dictionary& dict = dictptr->dict();
141  dict.add("method", defaultMethod);
142  dict.add("numberOfSubdomains", n);
143 
144  // Inject coeffs dictionary too
145  if (subMethodCoeffsDict.size())
146  {
147  dict.add(subMethodCoeffsDict.dictName(), subMethodCoeffsDict);
148  }
149  }
150  }
151  else
152  {
153  // Specified by full dictionaries
154 
155  // Create editable methods dictionaries
156  // - Only consider sub-dictionaries with a "numberOfSubdomains" entry
157  // This automatically filters out any coeffs dictionaries
158 
159  for (const entry& dEntry : coeffsDict_)
160  {
161  word methodName;
162 
163  if
164  (
165  dEntry.isDict()
166  // non-recursive, no patterns
167  && dEntry.dict().found("numberOfSubdomains", keyType::LITERAL)
168  )
169  {
170  // No method specified? can use a default method?
171 
172  const bool addDefaultMethod
173  (
174  !(dEntry.dict().found("method", keyType::LITERAL))
175  && !defaultMethod.empty()
176  );
177 
178  entry* e = methodsDict_.add(dEntry);
179 
180  if (addDefaultMethod && e && e->isDict())
181  {
182  e->dict().add("method", defaultMethod);
183  }
184  }
185  }
186  }
187 }
188 
189 
190 void Foam::multiLevelDecomp::setMethods()
191 {
192  // Assuming methodsDict_ has be properly created, convert the method
193  // dictionaries to actual methods
194 
195  label nLevels = 0;
196 
197  methods_.clear();
198  methods_.setSize(methodsDict_.size());
199  for (const entry& dEntry : methodsDict_)
200  {
201  // Dictionary entries only
202  // - these method dictionaries are non-regional
203  if (dEntry.isDict())
204  {
205  methods_.set
206  (
207  nLevels++,
208  // non-verbose would be nicer
209  decompositionMethod::New(dEntry.dict())
210  );
211  }
212  }
213 
214  methods_.setSize(nLevels);
215 
216  // Verify that nTotal is correct based on what each method delivers
217 
218  Info<< nl
219  << "Decompose " << type() << " [" << nDomains() << "] in "
220  << nLevels << " levels:" << endl;
221 
222  label nTotal = 1;
223  forAll(methods_, i)
224  {
225  Info<< " level " << i << " : " << methods_[i].type()
226  << " [" << methods_[i].nDomains() << "]" << endl;
227 
228  nTotal *= methods_[i].nDomains();
229  }
230 
231  if (nTotal != nDomains())
232  {
234  << "Top level decomposition specifies " << nDomains()
235  << " domains which is not equal to the product of"
236  << " all sub domains " << nTotal
237  << exit(FatalError);
238  }
239 }
240 
241 
242 // Given a subset of cells determine the new global indices. The problem
243 // is in the cells from neighbouring processors which need to be renumbered.
244 void Foam::multiLevelDecomp::subsetGlobalCellCells
245 (
246  const label nDomains,
247  const label domainI,
248  const labelList& dist,
249 
250  const labelListList& cellCells,
251  const labelList& set,
252  labelListList& subCellCells,
253  labelList& cutConnections
254 ) const
255 {
256  // Determine new index for cells by inverting subset
257  labelList oldToNew(invert(cellCells.size(), set));
258 
259  globalIndex globalCells(cellCells.size());
260 
261  // Subset locally the elements for which I have data
262  subCellCells = UIndirectList<labelList>(cellCells, set);
263 
264  // Get new indices for neighbouring processors
265  List<Map<label>> compactMap;
266  mapDistribute map(globalCells, subCellCells, compactMap);
267  map.distribute(oldToNew);
268  labelList allDist(dist);
269  map.distribute(allDist);
270 
271  // Now we have:
272  // oldToNew : the locally-compact numbering of all our cellCells. -1 if
273  // cellCell is not in set.
274  // allDist : destination domain for all our cellCells
275  // subCellCells : indexes into oldToNew and allDist
276 
277  // Globally compact numbering for cells in set.
278  globalIndex globalSubCells(set.size());
279 
280  // Now subCellCells contains indices into oldToNew which are the
281  // new locations of the neighbouring cells.
282 
283  cutConnections.setSize(nDomains);
284  cutConnections = 0;
285 
286  forAll(subCellCells, subCelli)
287  {
288  labelList& cCells = subCellCells[subCelli];
289 
290  // Keep the connections to valid mapped cells
291  label newI = 0;
292  forAll(cCells, i)
293  {
294  // Get locally-compact cell index of neighbouring cell
295  const label nbrCelli = oldToNew[cCells[i]];
296  if (nbrCelli == -1)
297  {
298  cutConnections[allDist[cCells[i]]]++;
299  }
300  else
301  {
302  // Reconvert local cell index into global one
303 
304  // Get original neighbour
305  const label celli = set[subCelli];
306  const label oldNbrCelli = cellCells[celli][i];
307  // Get processor from original neighbour
308  const label proci = globalCells.whichProcID(oldNbrCelli);
309  // Convert into global compact numbering
310  cCells[newI++] = globalSubCells.toGlobal(proci, nbrCelli);
311  }
312  }
313  cCells.setSize(newI);
314  }
315 }
316 
317 
318 void Foam::multiLevelDecomp::decompose
319 (
320  const labelListList& pointPoints,
321  const pointField& points,
322  const scalarField& pointWeights,
323  const labelUList& pointMap, // map back to original points
324  const label currLevel,
325  const label leafOffset,
326 
327  labelList& finalDecomp
328 ) const
329 {
330  labelList dist
331  (
332  methods_[currLevel].decompose
333  (
334  pointPoints,
335  points,
336  pointWeights
337  )
338  );
339 
340  // The next recursion level
341  const label nextLevel = currLevel+1;
342 
343  // Number of domains at this current level
344  const label nCurrDomains = methods_[currLevel].nDomains();
345 
346  // Calculate the domain remapping.
347  // The decompose() method delivers a distribution of [0..nDomains-1]
348  // which we map to the final location according to the decomposition
349  // leaf we are on.
350 
351  labelList domainLookup(nCurrDomains);
352  {
353  label sizes = 1; // Cumulative number of domains
354  for (label i = 0; i <= currLevel; ++i)
355  {
356  sizes *= methods_[i].nDomains();
357  }
358 
359  // Distribution of domains at this level
360  sizes = this->nDomains() / sizes;
361 
362  forAll(domainLookup, i)
363  {
364  domainLookup[i] = i * sizes + leafOffset;
365  }
366  }
367 
368  if (debug)
369  {
370  Info<< "Distribute at level " << currLevel
371  << " to domains" << nl
372  << flatOutput(domainLookup) << endl;
373  }
374 
375  // Extract processor+local index from point-point addressing
376  forAll(pointMap, i)
377  {
378  const label orig = pointMap[i];
379  finalDecomp[orig] = domainLookup[dist[i]];
380  }
381 
382  if (nextLevel < methods_.size())
383  {
384  // Recurse
385 
386  // Determine points per domain
387  labelListList domainToPoints(invertOneToMany(nCurrDomains, dist));
388 
389  // Extract processor+local index from point-point addressing
390  if (debug && Pstream::master())
391  {
392  Pout<< "Decomposition at level " << currLevel << " :" << endl;
393  }
394 
395  for (label domainI = 0; domainI < nCurrDomains; ++domainI)
396  {
397  // Extract elements for current domain
398  const labelList domainPoints(findIndices(dist, domainI));
399 
400  // Subset point-wise data.
401  pointField subPoints(points, domainPoints);
402  scalarField subWeights(pointWeights, domainPoints);
403  labelList subPointMap(labelUIndList(pointMap, domainPoints));
404  // Subset point-point addressing (adapt global numbering)
405  labelListList subPointPoints;
406  labelList nOutsideConnections;
407  subsetGlobalCellCells
408  (
409  nCurrDomains,
410  domainI,
411  dist,
412 
413  pointPoints,
414  domainPoints,
415 
416  subPointPoints,
417  nOutsideConnections
418  );
419 
420  label nPoints = returnReduce(domainPoints.size(), plusOp<label>());
421  Pstream::listCombineGather(nOutsideConnections, plusEqOp<label>());
422  Pstream::listCombineScatter(nOutsideConnections);
423  label nPatches = 0;
424  label nFaces = 0;
425  for (const label nConnect : nOutsideConnections)
426  {
427  if (nConnect > 0)
428  {
429  ++nPatches;
430  nFaces += nConnect;
431  }
432  }
433 
434  string oldPrefix;
435  if (debug && Pstream::master())
436  {
437  Pout<< " Domain " << domainI << nl
438  << " Number of cells = " << nPoints << nl
439  << " Number of inter-domain patches = " << nPatches
440  << nl
441  << " Number of inter-domain faces = " << nFaces << nl
442  << endl;
443  oldPrefix = Pout.prefix();
444  Pout.prefix() = " " + oldPrefix;
445  }
446 
447  decompose
448  (
449  subPointPoints,
450  subPoints,
451  subWeights,
452  subPointMap,
453  nextLevel,
454  domainLookup[domainI], // The offset for this level and leaf
455 
456  finalDecomp
457  );
458  if (debug && Pstream::master())
459  {
460  Pout.prefix() = oldPrefix;
461  }
462  }
463 
464 
465  if (debug)
466  {
467  // Do straight decompose of two levels
468  const label nNext = methods_[nextLevel].nDomains();
469  const label nTotal = nCurrDomains * nNext;
470 
471  // Get original level0 dictionary and modify numberOfSubdomains
472  dictionary level0Dict;
473  for (const entry& dEntry : methodsDict_)
474  {
475  if (dEntry.isDict())
476  {
477  level0Dict = dEntry.dict();
478  break;
479  }
480  }
481  level0Dict.set("numberOfSubdomains", nTotal);
482 
483  if (debug && Pstream::master())
484  {
485  Pout<< "Reference decomposition with " << level0Dict << " :"
486  << endl;
487  }
488 
489  autoPtr<decompositionMethod> method0 = decompositionMethod::New
490  (
491  level0Dict
492  );
493  labelList dist
494  (
495  method0().decompose
496  (
497  pointPoints,
498  points,
499  pointWeights
500  )
501  );
502 
503  for (label blockI = 0; blockI < nCurrDomains; ++blockI)
504  {
505  // Count the number inbetween blocks of nNext size
506 
507  label nPoints = 0;
508  labelList nOutsideConnections(nCurrDomains, Zero);
509  forAll(pointPoints, pointi)
510  {
511  if ((dist[pointi] / nNext) == blockI)
512  {
513  nPoints++;
514 
515  const labelList& pPoints = pointPoints[pointi];
516 
517  forAll(pPoints, i)
518  {
519  const label distBlockI = dist[pPoints[i]] / nNext;
520  if (distBlockI != blockI)
521  {
522  nOutsideConnections[distBlockI]++;
523  }
524  }
525  }
526  }
527 
528  reduce(nPoints, plusOp<label>());
530  (
531  nOutsideConnections,
532  plusEqOp<label>()
533  );
534  Pstream::listCombineScatter(nOutsideConnections);
535  label nPatches = 0;
536  label nFaces = 0;
537  for (const label nConnect : nOutsideConnections)
538  {
539  if (nConnect > 0)
540  {
541  ++nPatches;
542  nFaces += nConnect;
543  }
544  }
545 
546  if (debug && Pstream::master())
547  {
548  Pout<< " Domain " << blockI << nl
549  << " Number of cells = " << nPoints << nl
550  << " Number of inter-domain patches = "
551  << nPatches << nl
552  << " Number of inter-domain faces = " << nFaces
553  << nl << endl;
554  }
555  }
556  }
557  }
558 }
559 
560 
561 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
562 
563 Foam::multiLevelDecomp::multiLevelDecomp
564 (
565  const dictionary& decompDict,
566  const word& regionName
567 )
568 :
569  decompositionMethod(decompDict, regionName),
570  coeffsDict_
571  (
572  findCoeffsDict
573  (
574  typeName + "Coeffs",
575  (selectionType::EXACT | selectionType::MANDATORY)
576  )
577  ),
578  methodsDict_(),
579  methods_()
580 {
581  createMethodsDict();
582  setMethods();
583 }
584 
585 
586 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
587 
589 {
590  for (const decompositionMethod& meth : methods_)
591  {
592  if (!meth.parallelAware())
593  {
594  return false;
595  }
596  }
597 
598  return true;
599 }
600 
601 
602 Foam::labelList Foam::multiLevelDecomp::decompose
603 (
604  const polyMesh& mesh,
605  const pointField& cc,
606  const scalarField& cWeights
607 ) const
608 {
609  CompactListList<label> cellCells;
610  calcCellCells(mesh, identity(cc.size()), cc.size(), true, cellCells);
611 
612  labelList finalDecomp(cc.size(), Zero);
613  labelList cellMap(identity(cc.size()));
614 
615  decompose
616  (
617  cellCells(),
618  cc,
619  cWeights,
620  cellMap, // map back to original cells
621  0,
622  0,
623 
624  finalDecomp
625  );
626 
627  return finalDecomp;
628 }
629 
630 
631 Foam::labelList Foam::multiLevelDecomp::decompose
632 (
633  const labelListList& globalPointPoints,
634  const pointField& points,
635  const scalarField& pointWeights
636 ) const
637 {
638  labelList finalDecomp(points.size(), Zero);
639  labelList pointMap(identity(points.size()));
640 
641  decompose
642  (
643  globalPointPoints,
644  points,
645  pointWeights,
646  pointMap, // map back to original points
647  0,
648  0,
649 
650  finalDecomp
651  );
652 
653  return finalDecomp;
654 }
655 
656 
657 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
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
multiLevelDecomp.H
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::BitOps::set
void set(List< bool > &bools, const labelRange &range)
Set the specified range 'on' in a boolList.
Definition: BitOps.C:37
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
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::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::CompactListList
A packed storage unstructured matrix of objects of type <T> using an offset table for access.
Definition: CompactListList.H:63
globalIndex.H
Foam::UPstream::master
static bool master(const label communicator=worldComm)
Am I the master process.
Definition: UPstream.H:457
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::dictionary::set
entry * set(entry *entryPtr)
Assign a new entry, overwriting any existing entry.
Definition: dictionary.C:780
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::invert
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
Definition: ListOps.C:36
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
Foam::multiLevelDecomp::parallelAware
virtual bool parallelAware() const
Is method parallel aware?
Definition: multiLevelDecomp.C:588
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
n
label n
Definition: TABSMDCalcMethod2.H:31
regionName
Foam::word regionName
Definition: createNamedDynamicFvMesh.H:1
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:51
nPatches
const label nPatches
Definition: printMeshSummary.H:30
Foam::Field< vector >
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
IFstream.H
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
Foam::decompositionMethod
Abstract base class for domain decomposition.
Definition: decompositionMethod.H:51
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::flatOutput
FlatOutput::OutputAdaptor< Container, Delimiters > flatOutput(const Container &obj, Delimiters delim)
Global flatOutput() function with specified output delimiters.
Definition: FlatOutput.H:216
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::Pstream::listCombineGather
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
Definition: combineGatherScatter.C:290
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
Foam::nl
constexpr char nl
Definition: Ostream.H:404
mapDistribute.H
Foam::List< label >
Foam::type
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:590
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
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
Foam::decompositionMethod::New
static autoPtr< decompositionMethod > New(const dictionary &decompDict, const word &regionName="")
Definition: decompositionMethod.C:344
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Foam::decompositionMethod::nDomains
label nDomains() const noexcept
Number of domains.
Definition: decompositionMethod.H:209
Foam::findIndices
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
Foam::dictionary::add
entry * add(entry *entryPtr, bool mergeEntry=false)
Add a new entry.
Definition: dictionary.C:640
Foam::keyType::LITERAL
String literal.
Definition: keyType.H:81
Foam::Pstream::listCombineScatter
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
Definition: combineGatherScatter.C:432
Foam::labelUList
UList< label > labelUList
A UList of labels.
Definition: UList.H:85
Foam::invertOneToMany
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:114
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
Foam::dictionary::clear
void clear()
Clear the dictionary.
Definition: dictionary.C:857
Foam::decompositionMethod::findCoeffsDict
static const dictionary & findCoeffsDict(const dictionary &dict, const word &coeffsName, int select=selectionType::DEFAULT)
Definition: decompositionMethod.C:244
Foam::dictionary::readIfPresent
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:405
Foam::prefixOSstream::prefix
const string & prefix() const noexcept
Return the stream prefix.
Definition: prefixOSstream.H:101
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: UIndirectList.H:58