simpleGeomDecomp.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 "simpleGeomDecomp.H"
28 #include "SortableList.H"
29 #include "globalIndex.H"
30 #include "SubField.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(simpleGeomDecomp, 0);
37 
39  (
40  decompositionMethod,
41  simpleGeomDecomp,
42  decomposer
43  );
44 
46  (
47  decompositionMethod,
48  simpleGeomDecomp,
49  distributor
50  );
51 }
52 
53 
54 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
55 
56 // assignToProcessorGroup : given nCells cells and nProcGroup processor
57 // groups to share them, how do we share them out? Answer : each group
58 // gets nCells/nProcGroup cells, and the first few get one
59 // extra to make up the numbers. This should produce almost
60 // perfect load balancing
61 
62 void Foam::simpleGeomDecomp::assignToProcessorGroup
63 (
64  labelList& processorGroup,
65  const label nProcGroup
66 ) const
67 {
68  label jump = processorGroup.size()/nProcGroup;
69  label jumpb = jump + 1;
70  label fstProcessorGroup = processorGroup.size() - jump*nProcGroup;
71 
72  label ind = 0;
73  label j = 0;
74 
75  // assign cells to the first few processor groups (those with
76  // one extra cell each
77  for (j=0; j<fstProcessorGroup; j++)
78  {
79  for (label k=0; k<jumpb; k++)
80  {
81  processorGroup[ind++] = j;
82  }
83  }
84 
85  // and now to the `normal' processor groups
86  for (; j<nProcGroup; j++)
87  {
88  for (label k=0; k<jump; k++)
89  {
90  processorGroup[ind++] = j;
91  }
92  }
93 }
94 
95 
96 void Foam::simpleGeomDecomp::assignToProcessorGroup
97 (
98  labelList& processorGroup,
99  const label nProcGroup,
100  const labelList& indices,
101  const scalarField& weights,
102  const scalar summedWeights
103 ) const
104 {
105  // This routine gets the sorted points.
106  // Easiest to explain with an example.
107  // E.g. 400 points, sum of weights : 513.
108  // Now with number of divisions in this direction (nProcGroup) : 4
109  // gives the split at 513/4 = 128
110  // So summed weight from 0..128 goes into bin 0,
111  // ,, 128..256 goes into bin 1
112  // etc.
113  // Finally any remaining ones go into the last bin (3).
114 
115  const scalar jump = summedWeights/nProcGroup;
116  const label nProcGroupM1 = nProcGroup - 1;
117  scalar sumWeights = 0;
118  label ind = 0;
119  label j = 0;
120 
121  // assign cells to all except last group.
122  for (j=0; j<nProcGroupM1; j++)
123  {
124  const scalar limit = jump*scalar(j + 1);
125  while (sumWeights < limit)
126  {
127  sumWeights += weights[indices[ind]];
128  processorGroup[ind++] = j;
129  }
130  }
131  // Ensure last included.
132  while (ind < processorGroup.size())
133  {
134  processorGroup[ind++] = nProcGroupM1;
135  }
136 }
137 
138 
139 Foam::labelList Foam::simpleGeomDecomp::decomposeOneProc
140 (
141  const pointField& points
142 ) const
143 {
144  // construct a list for the final result
145  labelList finalDecomp(points.size());
146 
147  labelList processorGroups(points.size());
148 
149  labelList pointIndices(points.size());
150  forAll(pointIndices, i)
151  {
152  pointIndices[i] = i;
153  }
154 
155  const pointField rotatedPoints(rotDelta_ & points);
156 
157  // and one to take the processor group id's. For each direction.
158  // we assign the processors to groups of processors labelled
159  // 0..nX to give a banded structure on the mesh. Then we
160  // construct the actual processor number by treating this as
161  // the units part of the processor number.
162  sort
163  (
164  pointIndices,
165  UList<scalar>::less(rotatedPoints.component(vector::X))
166  );
167 
168  assignToProcessorGroup(processorGroups, n_.x());
169 
170  forAll(points, i)
171  {
172  finalDecomp[pointIndices[i]] = processorGroups[i];
173  }
174 
175 
176  // now do the same thing in the Y direction. These processor group
177  // numbers add multiples of nX to the proc. number (columns)
178  sort
179  (
180  pointIndices,
181  UList<scalar>::less(rotatedPoints.component(vector::Y))
182  );
183 
184  assignToProcessorGroup(processorGroups, n_.y());
185 
186  forAll(points, i)
187  {
188  finalDecomp[pointIndices[i]] += n_.x()*processorGroups[i];
189  }
190 
191 
192  // finally in the Z direction. Now we add multiples of nX*nY to give
193  // layers
194  sort
195  (
196  pointIndices,
197  UList<scalar>::less(rotatedPoints.component(vector::Z))
198  );
199 
200  assignToProcessorGroup(processorGroups, n_.z());
201 
202  forAll(points, i)
203  {
204  finalDecomp[pointIndices[i]] += n_.x()*n_.y()*processorGroups[i];
205  }
206 
207  return finalDecomp;
208 }
209 
210 
211 Foam::labelList Foam::simpleGeomDecomp::decomposeOneProc
212 (
213  const pointField& points,
214  const scalarField& weights
215 ) const
216 {
217  // construct a list for the final result
218  labelList finalDecomp(points.size());
219 
220  labelList processorGroups(points.size());
221 
222  labelList pointIndices(points.size());
223  forAll(pointIndices, i)
224  {
225  pointIndices[i] = i;
226  }
227 
228  const pointField rotatedPoints(rotDelta_ & points);
229 
230  // and one to take the processor group id's. For each direction.
231  // we assign the processors to groups of processors labelled
232  // 0..nX to give a banded structure on the mesh. Then we
233  // construct the actual processor number by treating this as
234  // the units part of the processor number.
235  sort
236  (
237  pointIndices,
238  UList<scalar>::less(rotatedPoints.component(vector::X))
239  );
240 
241  const scalar summedWeights = sum(weights);
242  assignToProcessorGroup
243  (
244  processorGroups,
245  n_.x(),
246  pointIndices,
247  weights,
248  summedWeights
249  );
250 
251  forAll(points, i)
252  {
253  finalDecomp[pointIndices[i]] = processorGroups[i];
254  }
255 
256 
257  // now do the same thing in the Y direction. These processor group
258  // numbers add multiples of nX to the proc. number (columns)
259  sort
260  (
261  pointIndices,
262  UList<scalar>::less(rotatedPoints.component(vector::Y))
263  );
264 
265  assignToProcessorGroup
266  (
267  processorGroups,
268  n_.y(),
269  pointIndices,
270  weights,
271  summedWeights
272  );
273 
274  forAll(points, i)
275  {
276  finalDecomp[pointIndices[i]] += n_.x()*processorGroups[i];
277  }
278 
279 
280  // finally in the Z direction. Now we add multiples of nX*nY to give
281  // layers
282  sort
283  (
284  pointIndices,
285  UList<scalar>::less(rotatedPoints.component(vector::Z))
286  );
287 
288  assignToProcessorGroup
289  (
290  processorGroups,
291  n_.z(),
292  pointIndices,
293  weights,
294  summedWeights
295  );
296 
297  forAll(points, i)
298  {
299  finalDecomp[pointIndices[i]] += n_.x()*n_.y()*processorGroups[i];
300  }
301 
302  return finalDecomp;
303 }
304 
305 
306 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
307 
309 :
310  geomDecomp(decompositionDict, typeName)
311 {}
312 
313 
314 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
315 
317 (
318  const pointField& points
319 )
320 {
321  if (!Pstream::parRun())
322  {
323  return decomposeOneProc(points);
324  }
325  else
326  {
327  globalIndex globalNumbers(points.size());
328 
329  // Collect all points on master
330  if (Pstream::master())
331  {
332  pointField allPoints(globalNumbers.size());
333 
334  label nTotalPoints = 0;
335  // Master first
336  SubField<point>(allPoints, points.size()) = points;
337  nTotalPoints += points.size();
338 
339  // Add slaves
340  for (int slave=1; slave<Pstream::nProcs(); slave++)
341  {
342  IPstream fromSlave(Pstream::commsTypes::scheduled, slave);
343  pointField nbrPoints(fromSlave);
345  (
346  allPoints,
347  nbrPoints.size(),
348  nTotalPoints
349  ) = nbrPoints;
350  nTotalPoints += nbrPoints.size();
351  }
352 
353  // Decompose
354  labelList finalDecomp(decomposeOneProc(allPoints));
355 
356  // Send back
357  for (int slave=1; slave<Pstream::nProcs(); slave++)
358  {
360  toSlave << SubField<label>
361  (
362  finalDecomp,
363  globalNumbers.localSize(slave),
364  globalNumbers.offset(slave)
365  );
366  }
367  // Get my own part
368  finalDecomp.setSize(points.size());
369 
370  return finalDecomp;
371  }
372  else
373  {
374  // Send my points
375  {
376  OPstream toMaster
377  (
380  );
381  toMaster<< points;
382  }
383 
384  // Receive back decomposition
385  IPstream fromMaster
386  (
389  );
390  labelList finalDecomp(fromMaster);
391 
392  return finalDecomp;
393  }
394  }
395 }
396 
397 
399 (
400  const pointField& points,
401  const scalarField& weights
402 )
403 {
404  if (!Pstream::parRun())
405  {
406  return decomposeOneProc(points, weights);
407  }
408  else
409  {
410  globalIndex globalNumbers(points.size());
411 
412  // Collect all points on master
413  if (Pstream::master())
414  {
415  pointField allPoints(globalNumbers.size());
416  scalarField allWeights(allPoints.size());
417 
418  label nTotalPoints = 0;
419  // Master first
420  SubField<point>(allPoints, points.size()) = points;
421  SubField<scalar>(allWeights, points.size()) = weights;
422  nTotalPoints += points.size();
423 
424  // Add slaves
425  for (int slave=1; slave<Pstream::nProcs(); slave++)
426  {
427  IPstream fromSlave(Pstream::commsTypes::scheduled, slave);
428  pointField nbrPoints(fromSlave);
429  scalarField nbrWeights(fromSlave);
431  (
432  allPoints,
433  nbrPoints.size(),
434  nTotalPoints
435  ) = nbrPoints;
437  (
438  allWeights,
439  nbrWeights.size(),
440  nTotalPoints
441  ) = nbrWeights;
442  nTotalPoints += nbrPoints.size();
443  }
444 
445  // Decompose
446  labelList finalDecomp(decomposeOneProc(allPoints, allWeights));
447 
448  // Send back
449  for (int slave=1; slave<Pstream::nProcs(); slave++)
450  {
452  toSlave << SubField<label>
453  (
454  finalDecomp,
455  globalNumbers.localSize(slave),
456  globalNumbers.offset(slave)
457  );
458  }
459  // Get my own part
460  finalDecomp.setSize(points.size());
461 
462  return finalDecomp;
463  }
464  else
465  {
466  // Send my points
467  {
468  OPstream toMaster
469  (
472  );
473  toMaster<< points << weights;
474  }
475 
476  // Receive back decomposition
477  IPstream fromMaster
478  (
481  );
482  labelList finalDecomp(fromMaster);
483 
484  return finalDecomp;
485  }
486  }
487 }
488 
489 
490 // ************************************************************************* //
#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
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
label k
Boltzmann constant.
Pre-declare related SubField type.
Definition: Field.H:60
static bool less(const vector &x, const vector &y)
To compare normals.
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Input inter-processor communications stream.
Definition: IPstream.H:50
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Geometrical domain decomposition.
Definition: geomDecomp.H:47
const pointField & points
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void sort(UList< T > &)
Definition: UList.C:115
List< label > labelList
A List of labels.
Definition: labelList.H:56
virtual labelList decompose(const pointField &)
Like decompose but with uniform weights on the points.
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
Output inter-processor communications stream.
Definition: OPstream.H:50
void setSize(const label)
Reset size of List.
Definition: List.C:281
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
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
simpleGeomDecomp(const dictionary &decompositionDict)
Construct given the decomposition dictionary.
Namespace for OpenFOAM.