simpleGeomDecomp.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-2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
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 "simpleGeomDecomp.H"
31#include "globalIndex.H"
32#include "SubField.H"
33
34// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35
36namespace Foam
37{
40 (
44 );
45}
46
47
48// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
49
50namespace Foam
51{
52
53// A list compare binary predicate for normal sort by vector component
55{
58
60 :
61 values(list),
62 sortCmpt(cmpt)
63 {}
64
66 {
67 sortCmpt = cmpt;
68 }
69
70 bool operator()(const label a, const label b) const
71 {
72 return values[a][sortCmpt] < values[b][sortCmpt];
73 }
74};
75
76} // End namespace Foam
77
78
79// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
80
81// assignToProcessorGroup : given nCells cells and nProcGroup processor
82// groups to share them, how do we share them out? Answer : each group
83// gets nCells/nProcGroup cells, and the first few get one
84// extra to make up the numbers. This should produce almost
85// perfect load balancing
86
87void Foam::simpleGeomDecomp::assignToProcessorGroup
88(
89 labelList& processorGroup,
90 const label nProcGroup
91)
92{
93 label jump = processorGroup.size()/nProcGroup;
94 label jumpb = jump + 1;
95 label fstProcessorGroup = processorGroup.size() - jump*nProcGroup;
96
97 label ind = 0;
98 label j = 0;
99
100 // assign cells to the first few processor groups (those with
101 // one extra cell each
102 for (j=0; j<fstProcessorGroup; j++)
103 {
104 for (label k=0; k<jumpb; k++)
105 {
106 processorGroup[ind++] = j;
107 }
108 }
109
110 // and now to the `normal' processor groups
111 for (; j<nProcGroup; j++)
112 {
113 for (label k=0; k<jump; k++)
114 {
115 processorGroup[ind++] = j;
116 }
117 }
118}
119
120
121void Foam::simpleGeomDecomp::assignToProcessorGroup
122(
123 labelList& processorGroup,
124 const label nProcGroup,
125 const labelList& indices,
126 const scalarField& weights,
127 const scalar summedWeights
128)
129{
130 // This routine gets the sorted points.
131 // Easiest to explain with an example.
132 // E.g. 400 points, sum of weights : 513.
133 // Now with number of divisions in this direction (nProcGroup) : 4
134 // gives the split at 513/4 = 128
135 // So summed weight from 0..128 goes into bin 0,
136 // ,, 128..256 goes into bin 1
137 // etc.
138 // Finally any remaining ones go into the last bin (3).
139
140 const scalar jump = summedWeights/nProcGroup;
141 const label nProcGroupM1 = nProcGroup - 1;
142 scalar sumWeights = 0;
143 label ind = 0;
144 label j = 0;
145
146 // assign cells to all except last group.
147 for (j=0; j<nProcGroupM1; j++)
148 {
149 const scalar limit = jump*scalar(j + 1);
150 while (sumWeights < limit)
151 {
152 sumWeights += weights[indices[ind]];
153 processorGroup[ind++] = j;
154 }
155 }
156 // Ensure last included.
157 while (ind < processorGroup.size())
158 {
159 processorGroup[ind++] = nProcGroupM1;
160 }
161}
162
163
164Foam::labelList Foam::simpleGeomDecomp::decomposeOneProc
165(
166 const pointField& points
167) const
168{
169 // construct a list for the final result
170 labelList finalDecomp(points.size());
171
172 labelList processorGroups(points.size());
173
174 labelList pointIndices(identity(points.size()));
175
176 const pointField rotatedPoints(adjustPoints(points));
177
178 vectorLessOp sorter(rotatedPoints);
179
180 // and one to take the processor group id's. For each direction.
181 // we assign the processors to groups of processors labelled
182 // 0..nX to give a banded structure on the mesh. Then we
183 // construct the actual processor number by treating this as
184 // the units part of the processor number.
185
186 sorter.setComponent(vector::X);
187 Foam::sort(pointIndices, sorter);
188
189 assignToProcessorGroup(processorGroups, n_.x());
190
191 forAll(points, i)
192 {
193 finalDecomp[pointIndices[i]] = processorGroups[i];
194 }
195
196
197 // now do the same thing in the Y direction. These processor group
198 // numbers add multiples of nX to the proc. number (columns)
199
200 sorter.setComponent(vector::Y);
201 Foam::sort(pointIndices, sorter);
202
203 assignToProcessorGroup(processorGroups, n_.y());
204
205 forAll(points, i)
206 {
207 finalDecomp[pointIndices[i]] += n_.x()*processorGroups[i];
208 }
209
210
211 // finally in the Z direction. Now we add multiples of nX*nY to give
212 // layers
213
214 sorter.setComponent(vector::Z);
215 Foam::sort(pointIndices, sorter);
216
217 assignToProcessorGroup(processorGroups, n_.z());
218
219 forAll(points, i)
220 {
221 finalDecomp[pointIndices[i]] += n_.x()*n_.y()*processorGroups[i];
222 }
223
224 return finalDecomp;
225}
226
227
228Foam::labelList Foam::simpleGeomDecomp::decomposeOneProc
229(
230 const pointField& points,
231 const scalarField& weights
232) const
233{
234 // construct a list for the final result
235 labelList finalDecomp(points.size());
236
237 labelList processorGroups(points.size());
238
239 labelList pointIndices(identity(points.size()));
240
241 const pointField rotatedPoints(adjustPoints(points));
242
243 vectorLessOp sorter(rotatedPoints);
244
245 // and one to take the processor group id's. For each direction.
246 // we assign the processors to groups of processors labelled
247 // 0..nX to give a banded structure on the mesh. Then we
248 // construct the actual processor number by treating this as
249 // the units part of the processor number.
250
251 sorter.setComponent(vector::X);
252 Foam::sort(pointIndices, sorter);
253
254 const scalar summedWeights = sum(weights);
255 assignToProcessorGroup
256 (
257 processorGroups,
258 n_.x(),
259 pointIndices,
260 weights,
261 summedWeights
262 );
263
264 forAll(points, i)
265 {
266 finalDecomp[pointIndices[i]] = processorGroups[i];
267 }
268
269
270 // now do the same thing in the Y direction. These processor group
271 // numbers add multiples of nX to the proc. number (columns)
272
273 sorter.setComponent(vector::Y);
274 Foam::sort(pointIndices, sorter);
275
276 assignToProcessorGroup
277 (
278 processorGroups,
279 n_.y(),
280 pointIndices,
281 weights,
282 summedWeights
283 );
284
285 forAll(points, i)
286 {
287 finalDecomp[pointIndices[i]] += n_.x()*processorGroups[i];
288 }
289
290
291 // finally in the Z direction. Now we add multiples of nX*nY to give
292 // layers
293
294 sorter.setComponent(vector::Z);
295 Foam::sort(pointIndices, sorter);
296
297 assignToProcessorGroup
298 (
299 processorGroups,
300 n_.z(),
301 pointIndices,
302 weights,
303 summedWeights
304 );
305
306 forAll(points, i)
307 {
308 finalDecomp[pointIndices[i]] += n_.x()*n_.y()*processorGroups[i];
309 }
310
311 return finalDecomp;
312}
313
314
315// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
316
318(
319 const dictionary& decompDict,
320 const word& regionName
321)
322:
323 geomDecomp(typeName, decompDict, regionName)
324{}
325
326
327// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
328
330(
331 const pointField& points
332) const
333{
334 if (!Pstream::parRun())
335 {
336 return decomposeOneProc(points);
337 }
338 else
339 {
340 const globalIndex globalNumbers(points.size());
341
342 pointField allPoints(globalNumbers.gather(points));
343
344 labelList allDecomp;
345 if (Pstream::master())
346 {
347 allDecomp = decomposeOneProc(allPoints);
348 allPoints.clear(); // Not needed anymore
349 }
350
351 return globalNumbers.scatter(allDecomp);
352 }
353}
354
355
357(
358 const pointField& points,
359 const scalarField& weights
360) const
361{
362 if (returnReduce((points.size() != weights.size()), orOp<bool>()))
363 {
364 // Ignore zero-sized weights ... and poorly sized ones too
365 return decompose(points);
366 }
367 else if (!Pstream::parRun())
368 {
369 return decomposeOneProc(points, weights);
370 }
371 else
372 {
373 const globalIndex globalNumbers(points.size());
374
375 pointField allPoints(globalNumbers.gather(points));
376 scalarField allWeights(globalNumbers.gather(weights));
377
378 labelList allDecomp;
379 if (Pstream::master())
380 {
381 allDecomp = decomposeOneProc(allPoints, allWeights);
382 allPoints.clear(); // Not needed anymore
383 allWeights.clear(); // ...
384 }
385
386 return globalNumbers.scatter(allDecomp);
387 }
388}
389
390
391// ************************************************************************* //
label k
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:94
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
Abstract base class for domain decomposition.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Base for geometrical domain decomposition methods.
Definition: geomDecomp.H:89
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
static void gather(const labelUList &offsets, const label comm, const ProcIDsContainer &procIDs, const UList< Type > &fld, List< Type > &allFld, const int tag=UPstream::msgType(), const UPstream::commsTypes=UPstream::commsTypes::nonBlocking)
Collect data in processor order on master (== procIDs[0]).
static void scatter(const labelUList &offsets, const label comm, const ProcIDsContainer &procIDs, const UList< Type > &allFld, UList< Type > &fld, const int tag=UPstream::msgType(), const UPstream::commsTypes=UPstream::commsTypes::nonBlocking)
Distribute data in processor order.
Simple geometric decomposition, selectable as simple.
splitCell * master() const
Definition: splitCell.H:113
bool decompose() const noexcept
Query the decompose flag (normally off)
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
Foam::word regionName(Foam::polyMesh::defaultRegion)
const pointField & points
Namespace for OpenFOAM.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
List< label > labelList
A List of labels.
Definition: List.H:66
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
complex limit(const complex &, const complex &)
Definition: complexI.H:263
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void sort(UList< T > &list)
Sort the list.
Definition: UList.C:342
uint8_t direction
Definition: direction.H:56
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
volScalarField & b
Definition: createFields.H:27
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
void setComponent(direction cmpt)
bool operator()(const label a, const label b) const
vectorLessOp(const UList< vector > &list, direction cmpt=vector::X)
const UList< vector > & values