multiLevelDecomp.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-2022 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 "multiLevelDecomp.H"
28 #include "IFstream.H"
29 #include "globalIndex.H"
30 #include "distributionMap.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(multiLevelDecomp, 0);
37 
39  (
40  decompositionMethod,
41  multiLevelDecomp,
42  decomposer
43  );
44 
46  (
47  decompositionMethod,
48  multiLevelDecomp,
49  distributor
50  );
51 }
52 
53 
54 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
55 
56 // Given a subset of cells determine the new global indices. The problem
57 // is in the cells from neighbouring processors which need to be renumbered.
58 void Foam::multiLevelDecomp::subsetGlobalCellCells
59 (
60  const label nDomains,
61  const label domainI,
62  const labelList& dist,
63 
64  const labelListList& cellCells,
65  const labelList& set,
66  labelListList& subCellCells,
67  labelList& cutConnections
68 ) const
69 {
70  // Determine new index for cells by inverting subset
71  labelList oldToNew(invert(cellCells.size(), set));
72 
73  globalIndex globalCells(cellCells.size());
74 
75  // Subset locally the elements for which I have data
76  subCellCells = UIndirectList<labelList>(cellCells, set);
77 
78  // Get new indices for neighbouring processors
79  List<Map<label>> compactMap;
80  distributionMap map(globalCells, subCellCells, compactMap);
81  map.distribute(oldToNew);
82  labelList allDist(dist);
83  map.distribute(allDist);
84 
85  // Now we have:
86  // oldToNew : the locally-compact numbering of all our cellCells. -1 if
87  // cellCell is not in set.
88  // allDist : destination domain for all our cellCells
89  // subCellCells : indexes into oldToNew and allDist
90 
91  // Globally compact numbering for cells in set.
92  globalIndex globalSubCells(set.size());
93 
94  // Now subCellCells contains indices into oldToNew which are the
95  // new locations of the neighbouring cells.
96 
97  cutConnections.setSize(nDomains);
98  cutConnections = 0;
99 
100  forAll(subCellCells, subCelli)
101  {
102  labelList& cCells = subCellCells[subCelli];
103 
104  // Keep the connections to valid mapped cells
105  label newI = 0;
106  forAll(cCells, i)
107  {
108  // Get locally-compact cell index of neighbouring cell
109  label nbrCelli = oldToNew[cCells[i]];
110  if (nbrCelli == -1)
111  {
112  cutConnections[allDist[cCells[i]]]++;
113  }
114  else
115  {
116  // Reconvert local cell index into global one
117 
118  // Get original neighbour
119  label celli = set[subCelli];
120  label oldNbrCelli = cellCells[celli][i];
121  // Get processor from original neighbour
122  label proci = globalCells.whichProcID(oldNbrCelli);
123  // Convert into global compact numbering
124  cCells[newI++] = globalSubCells.toGlobal(proci, nbrCelli);
125  }
126  }
127  cCells.setSize(newI);
128  }
129 }
130 
131 
132 void Foam::multiLevelDecomp::decompose
133 (
134  const labelListList& pointPoints,
135  const pointField& points,
136  const scalarField& pointWeights,
137  const labelList& pointMap, // map back to original points
138  const label levelI,
139 
140  labelField& finalDecomp
141 )
142 {
143  labelList dist
144  (
145  methods_[levelI].decompose
146  (
147  pointPoints,
148  points,
149  pointWeights
150  )
151  );
152 
153  forAll(pointMap, i)
154  {
155  label orig = pointMap[i];
156  finalDecomp[orig] += dist[i];
157  }
158 
159  if (levelI != methods_.size()-1)
160  {
161  // Recurse
162 
163  // Determine points per domain
164  label n = methods_[levelI].nDomains();
165  labelListList domainToPoints(invertOneToMany(n, dist));
166 
167  // 'Make space' for new levels of decomposition
168  finalDecomp *= methods_[levelI+1].nDomains();
169 
170  // Extract processor+local index from point-point addressing
171  if (debug && Pstream::master())
172  {
173  Pout<< "Decomposition at level " << levelI << " :" << endl;
174  }
175 
176  for (label domainI = 0; domainI < n; domainI++)
177  {
178  // Extract elements for current domain
179  const labelList domainPoints(findIndices(dist, domainI));
180 
181  // Subset point-wise data.
182  pointField subPoints(points, domainPoints);
183  scalarField subWeights(pointWeights, domainPoints);
184  labelList subPointMap(UIndirectList<label>(pointMap, domainPoints));
185  // Subset point-point addressing (adapt global numbering)
186  labelListList subPointPoints;
187  labelList nOutsideConnections;
188  subsetGlobalCellCells
189  (
190  n,
191  domainI,
192  dist,
193 
194  pointPoints,
195  domainPoints,
196 
197  subPointPoints,
198  nOutsideConnections
199  );
200 
201  label nPoints = returnReduce(domainPoints.size(), plusOp<label>());
202  Pstream::listCombineGather(nOutsideConnections, plusEqOp<label>());
203  Pstream::listCombineScatter(nOutsideConnections);
204  label nPatches = 0;
205  label nFaces = 0;
206  forAll(nOutsideConnections, i)
207  {
208  if (nOutsideConnections[i] > 0)
209  {
210  nPatches++;
211  nFaces += nOutsideConnections[i];
212  }
213  }
214 
215  string oldPrefix;
216  if (debug && Pstream::master())
217  {
218  Pout<< " Domain " << domainI << nl
219  << " Number of cells = " << nPoints << nl
220  << " Number of inter-domain patches = " << nPatches
221  << nl
222  << " Number of inter-domain faces = " << nFaces << nl
223  << endl;
224  oldPrefix = Pout.prefix();
225  Pout.prefix() = " " + oldPrefix;
226  }
227 
228  decompose
229  (
230  subPointPoints,
231  subPoints,
232  subWeights,
233  subPointMap,
234  levelI+1,
235 
236  finalDecomp
237  );
238  if (debug && Pstream::master())
239  {
240  Pout.prefix() = oldPrefix;
241  }
242  }
243 
244 
245  if (debug)
246  {
247  // Do straight decompose of two levels
248  label nNext = methods_[levelI+1].nDomains();
249  label nTotal = n*nNext;
250 
251  // Retrieve original level0 dictionary and modify number of domains
253  decompositionDict_.optionalSubDict(typeName + "Coeffs").begin();
254  dictionary myDict = iter().dict();
255  myDict.set("numberOfSubdomains", nTotal);
256 
257  if (debug && Pstream::master())
258  {
259  Pout<< "Reference decomposition with " << myDict << " :"
260  << endl;
261  }
262 
263  autoPtr<decompositionMethod> method0 =
265  (
266  myDict
267  );
268 
269  labelList dist
270  (
271  method0().decompose
272  (
273  pointPoints,
274  points,
275  pointWeights
276  )
277  );
278 
279  for (label blockI = 0; blockI < n; blockI++)
280  {
281  // Count the number in between blocks of nNext size
282 
283  label nPoints = 0;
284  labelList nOutsideConnections(n, 0);
285  forAll(pointPoints, pointi)
286  {
287  if ((dist[pointi] / nNext) == blockI)
288  {
289  nPoints++;
290 
291  const labelList& pPoints = pointPoints[pointi];
292 
293  forAll(pPoints, i)
294  {
295  label distBlockI = dist[pPoints[i]] / nNext;
296  if (distBlockI != blockI)
297  {
298  nOutsideConnections[distBlockI]++;
299  }
300  }
301  }
302  }
303 
304  reduce(nPoints, plusOp<label>());
306  (
307  nOutsideConnections,
308  plusEqOp<label>()
309  );
310  Pstream::listCombineScatter(nOutsideConnections);
311  label nPatches = 0;
312  label nFaces = 0;
313  forAll(nOutsideConnections, i)
314  {
315  if (nOutsideConnections[i] > 0)
316  {
317  nPatches++;
318  nFaces += nOutsideConnections[i];
319  }
320  }
321 
322  if (debug && Pstream::master())
323  {
324  Pout<< " Domain " << blockI << nl
325  << " Number of cells = " << nPoints << nl
326  << " Number of inter-domain patches = "
327  << nPatches << nl
328  << " Number of inter-domain faces = " << nFaces
329  << nl << endl;
330  }
331  }
332  }
333  }
334 }
335 
336 
337 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
338 
340 :
341  decompositionMethod(decompositionDict),
342  methodsDict_(decompositionDict_.optionalSubDict(typeName + "Coeffs"))
343 {
344  methods_.setSize(methodsDict_.size());
345  label i = 0;
346  forAllConstIter(dictionary, methodsDict_, iter)
347  {
348  methods_.set(i++, decompositionMethod::NewDecomposer(iter().dict()));
349  }
350 
351  label n = 1;
352  Info<< "decompositionMethod " << type() << " :" << endl;
353  forAll(methods_, i)
354  {
355  Info<< " level " << i << " decomposing with " << methods_[i].type()
356  << " into " << methods_[i].nDomains() << " subdomains." << endl;
357 
358  n *= methods_[i].nDomains();
359  }
360 
361  if (n != nDomains())
362  {
364  << "Top level decomposition specifies " << nDomains()
365  << " domains which is not equal to the product of"
366  << " all sub domains " << n
367  << exit(FatalError);
368  }
369 }
370 
371 
372 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
373 
374 Foam::labelList Foam::multiLevelDecomp::decompose
375 (
376  const polyMesh& mesh,
377  const pointField& cc,
378  const scalarField& cWeights
379 )
380 {
381  CompactListList<label> cellCells;
382  calcCellCells(mesh, identity(cc.size()), cc.size(), true, cellCells);
383 
384  labelField finalDecomp(cc.size(), 0);
385  labelList cellMap(identity(cc.size()));
386 
387  decompose
388  (
389  cellCells.list(),
390  cc,
391  cWeights,
392  cellMap, // map back to original cells
393  0,
394 
395  finalDecomp
396  );
397 
398  return move(finalDecomp);
399 }
400 
401 
402 Foam::labelList Foam::multiLevelDecomp::decompose
403 (
404  const labelListList& globalPointPoints,
405  const pointField& points,
406  const scalarField& pointWeights
407 )
408 {
409  labelField finalDecomp(points.size(), 0);
410  labelList pointMap(identity(points.size()));
411 
412  decompose
413  (
414  globalPointPoints,
415  points,
416  pointWeights,
417  pointMap, // map back to original points
418  0,
419 
420  finalDecomp
421  );
422 
423  return move(finalDecomp);
424 }
425 
426 
427 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
dictionary dict
Field< label > labelField
Specialisation of Field<T> for label.
Definition: labelField.H:49
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
#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.
multiLevelDecomp(const dictionary &decompositionDict)
Construct given the decomposition dictionary.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
labelListList invertOneToMany(const label len, const labelUList &)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:67
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
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
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurrences of given element. Linear search.
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
Macros for easy insertion into run-time selection tables.
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.
const string & prefix() const
Return the prefix of the stream.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< Container > list() const
Convert to List<Container>
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
Definition: ListOps.C:37
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static autoPtr< decompositionMethod > NewDecomposer(const dictionary &decompositionDict)
Return a reference to the selected decomposition method.
List< label > labelList
A List of labels.
Definition: labelList.H:56
Abstract base class for decomposition.
A packed storage unstructured matrix of objects of type <T> using an offset table for access...
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
static const char nl
Definition: Ostream.H:260
defineTypeNameAndDebug(combustionModel, 0)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
messageStream Info
label n
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
friend class const_iterator
Definition: UILList.H:78
Namespace for OpenFOAM.