hierarchGeomDecomp.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-2019 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 "hierarchGeomDecomp.H"
28 #include "PstreamReduceOps.H"
29 #include "SortableList.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35  defineTypeNameAndDebug(hierarchGeomDecomp, 0);
36 
38  (
39  decompositionMethod,
40  hierarchGeomDecomp,
41  dictionary
42  );
43 }
44 
45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
46 
47 void Foam::hierarchGeomDecomp::setDecompOrder()
48 {
49  const word order(geomDecomDict_.lookup("order"));
50 
51  if (order.size() != 3)
52  {
54  (
55  decompositionDict_
56  ) << "number of characters in order (" << order << ") != 3"
57  << exit(FatalIOError);
58  }
59 
60  for (label i = 0; i < 3; ++i)
61  {
62  if (order[i] == 'x')
63  {
64  decompOrder_[i] = 0;
65  }
66  else if (order[i] == 'y')
67  {
68  decompOrder_[i] = 1;
69  }
70  else if (order[i] == 'z')
71  {
72  decompOrder_[i] = 2;
73  }
74  else
75  {
77  (
78  decompositionDict_
79  ) << "Illegal decomposition order " << order << endl
80  << "It should only contain x, y or z" << exit(FatalError);
81  }
82  }
83 }
84 
85 
86 Foam::label Foam::hierarchGeomDecomp::findLower
87 (
88  const List<scalar>& l,
89  const scalar t,
90  const label initLow,
91  const label initHigh
92 )
93 {
94  if (initHigh <= initLow)
95  {
96  return initLow;
97  }
98 
99  label low = initLow;
100  label high = initHigh;
101 
102  while ((high - low) > 1)
103  {
104  label mid = (low + high)/2;
105 
106  if (l[mid] < t)
107  {
108  low = mid;
109  }
110  else
111  {
112  high = mid;
113  }
114  }
115 
116  // high and low can still differ by one. Choose best.
117 
118  label tIndex = -1;
119 
120  if (l[high-1] < t)
121  {
122  tIndex = high;
123  }
124  else
125  {
126  tIndex = low;
127  }
128 
129  return tIndex;
130 }
131 
132 
133 // Create a mapping between the index and the weighted size.
134 // For convenience, sortedWeightedSize is one size bigger than current. This
135 // avoids extra tests.
136 void Foam::hierarchGeomDecomp::calculateSortedWeightedSizes
137 (
138  const labelList& current,
139  const labelList& indices,
140  const scalarField& weights,
141  const label globalCurrentSize,
142 
143  scalarField& sortedWeightedSizes
144 )
145 {
146  // Evaluate cumulative weights.
147  sortedWeightedSizes[0] = 0;
148  forAll(current, i)
149  {
150  label pointi = current[indices[i]];
151  sortedWeightedSizes[i + 1] = sortedWeightedSizes[i] + weights[pointi];
152  }
153  // Non-dimensionalise and multiply by size.
154  scalar globalCurrentLength = returnReduce
155  (
156  sortedWeightedSizes[current.size()],
157  sumOp<scalar>()
158  );
159  // Normalise weights by global sum of weights and multiply through
160  // by global size.
161  sortedWeightedSizes *= (globalCurrentSize/globalCurrentLength);
162 }
163 
164 
165 // Find position in values so between minIndex and this position there
166 // are wantedSize elements.
167 void Foam::hierarchGeomDecomp::findBinary
168 (
169  const label sizeTol,
170  const List<scalar>& values,
171  const label minIndex, // index of previous value
172  const scalar minValue, // value at minIndex
173  const scalar maxValue, // global max of values
174  const scalar wantedSize, // wanted size
175 
176  label& mid, // index where size of bin is
177  // wantedSize (to within sizeTol)
178  scalar& midValue // value at mid
179 )
180 {
181  label low = minIndex;
182  scalar lowValue = minValue;
183 
184  scalar highValue = maxValue;
185  // (one beyond) index of highValue
186  label high = values.size();
187 
188  // Safeguards to avoid infinite loop.
189  scalar midValuePrev = vGreat;
190 
191  while (true)
192  {
193  label size = returnReduce(mid-minIndex, sumOp<label>());
194 
195  if (debug)
196  {
197  Pout<< " low:" << low << " lowValue:" << lowValue
198  << " high:" << high << " highValue:" << highValue
199  << " mid:" << mid << " midValue:" << midValue << endl
200  << " globalSize:" << size << " wantedSize:" << wantedSize
201  << " sizeTol:" << sizeTol << endl;
202  }
203 
204  if (wantedSize < size - sizeTol)
205  {
206  high = mid;
207  highValue = midValue;
208  }
209  else if (wantedSize > size + sizeTol)
210  {
211  low = mid;
212  lowValue = midValue;
213  }
214  else
215  {
216  break;
217  }
218 
219  // Update mid, midValue
220  midValue = 0.5*(lowValue+highValue);
221  mid = findLower(values, midValue, low, high);
222 
223  // Safeguard if same as previous.
224  bool hasNotChanged = (mag(midValue-midValuePrev) < small);
225 
226  if (returnReduce(hasNotChanged, andOp<bool>()))
227  {
229  << "unable to find desired decomposition split, making do!"
230  << endl;
231  break;
232  }
233 
234  midValuePrev = midValue;
235  }
236 }
237 
238 
239 // Find position in values so between minIndex and this position there
240 // are wantedSize elements.
241 void Foam::hierarchGeomDecomp::findBinary
242 (
243  const label sizeTol,
244  const List<scalar>& sortedWeightedSizes,
245  const List<scalar>& values,
246  const label minIndex, // index of previous value
247  const scalar minValue, // value at minIndex
248  const scalar maxValue, // global max of values
249  const scalar wantedSize, // wanted size
250 
251  label& mid, // index where size of bin is
252  // wantedSize (to within sizeTol)
253  scalar& midValue // value at mid
254 )
255 {
256  label low = minIndex;
257  scalar lowValue = minValue;
258 
259  scalar highValue = maxValue;
260  // (one beyond) index of highValue
261  label high = values.size();
262 
263  // Safeguards to avoid infinite loop.
264  scalar midValuePrev = vGreat;
265 
266  while (true)
267  {
268  scalar weightedSize = returnReduce
269  (
270  sortedWeightedSizes[mid] - sortedWeightedSizes[minIndex],
271  sumOp<scalar>()
272  );
273 
274  if (debug)
275  {
276  Pout<< " low:" << low << " lowValue:" << lowValue
277  << " high:" << high << " highValue:" << highValue
278  << " mid:" << mid << " midValue:" << midValue << endl
279  << " globalSize:" << weightedSize
280  << " wantedSize:" << wantedSize
281  << " sizeTol:" << sizeTol << endl;
282  }
283 
284  if (wantedSize < weightedSize - sizeTol)
285  {
286  high = mid;
287  highValue = midValue;
288  }
289  else if (wantedSize > weightedSize + sizeTol)
290  {
291  low = mid;
292  lowValue = midValue;
293  }
294  else
295  {
296  break;
297  }
298 
299  // Update mid, midValue
300  midValue = 0.5*(lowValue+highValue);
301  mid = findLower(values, midValue, low, high);
302 
303  // Safeguard if same as previous.
304  bool hasNotChanged = (mag(midValue-midValuePrev) < small);
305 
306  if (returnReduce(hasNotChanged, andOp<bool>()))
307  {
309  << "unable to find desired decomposition split, making do!"
310  << endl;
311  break;
312  }
313 
314  midValuePrev = midValue;
315  }
316 }
317 
318 
319 // Sort points into bins according to one component. Recurses to next component.
320 void Foam::hierarchGeomDecomp::sortComponent
321 (
322  const label sizeTol,
323  const pointField& points,
324  const labelList& current, // slice of points to decompose
325  const direction componentIndex, // index in decompOrder_
326  const label mult, // multiplication factor for finalDecomp
327  labelList& finalDecomp
328 )
329 {
330  // Current component
331  label compI = decompOrder_[componentIndex];
332 
333  if (debug)
334  {
335  Pout<< "sortComponent : Sorting slice of size " << current.size()
336  << " in component " << compI << endl;
337  }
338 
339  // Storage for sorted component compI
340  SortableList<scalar> sortedCoord(current.size());
341 
342  forAll(current, i)
343  {
344  label pointi = current[i];
345 
346  sortedCoord[i] = points[pointi][compI];
347  }
348  sortedCoord.sort();
349 
350  label globalCurrentSize = returnReduce(current.size(), sumOp<label>());
351 
352  scalar minCoord = returnReduce
353  (
354  (
355  sortedCoord.size()
356  ? sortedCoord[0]
357  : great
358  ),
359  minOp<scalar>()
360  );
361 
362  scalar maxCoord = returnReduce
363  (
364  (
365  sortedCoord.size()
366  ? sortedCoord.last()
367  : -great
368  ),
369  maxOp<scalar>()
370  );
371 
372  if (debug)
373  {
374  Pout<< "sortComponent : minCoord:" << minCoord
375  << " maxCoord:" << maxCoord << endl;
376  }
377 
378  // starting index (in sortedCoord) of bin (= local)
379  label leftIndex = 0;
380  // starting value of bin (= global since coordinate)
381  scalar leftCoord = minCoord;
382 
383  // Sort bins of size n
384  for (label bin = 0; bin < n_[compI]; bin++)
385  {
386  // Now we need to determine the size of the bin (dx). This is
387  // determined by the 'pivot' values - everything to the left of this
388  // value goes in the current bin, everything to the right into the next
389  // bins.
390 
391  // Local number of elements
392  label localSize = -1; // offset from leftOffset
393 
394  // Value at right of bin (leftIndex+localSize-1)
395  scalar rightCoord = -great;
396 
397  if (bin == n_[compI]-1)
398  {
399  // Last bin. Copy all.
400  localSize = current.size()-leftIndex;
401  rightCoord = maxCoord; // note: not used anymore
402  }
403  else if (Pstream::nProcs() == 1)
404  {
405  // No need for binary searching of bin size
406  localSize = label(current.size()/n_[compI]);
407  rightCoord = sortedCoord[leftIndex+localSize];
408  }
409  else
410  {
411  // For the current bin (starting at leftCoord) we want a rightCoord
412  // such that the sum of all sizes are globalCurrentSize/n_[compI].
413  // We have to iterate to obtain this.
414 
415  label rightIndex = current.size();
416  rightCoord = maxCoord;
417 
418  // Calculate rightIndex/rightCoord to have wanted size
419  findBinary
420  (
421  sizeTol,
422  sortedCoord,
423  leftIndex,
424  leftCoord,
425  maxCoord,
426  globalCurrentSize/n_[compI], // wanted size
427 
428  rightIndex,
429  rightCoord
430  );
431  localSize = rightIndex - leftIndex;
432  }
433 
434  if (debug)
435  {
436  Pout<< "For component " << compI << ", bin " << bin
437  << " copying" << endl
438  << "from " << leftCoord << " at local index "
439  << leftIndex << endl
440  << "to " << rightCoord << " localSize:"
441  << localSize << endl
442  << endl;
443  }
444 
445 
446  // Copy localSize elements starting from leftIndex.
447  labelList slice(localSize);
448 
449  forAll(slice, i)
450  {
451  label pointi = current[sortedCoord.indices()[leftIndex+i]];
452 
453  // Mark point into correct bin
454  finalDecomp[pointi] += bin*mult;
455 
456  // And collect for next sorting action
457  slice[i] = pointi;
458  }
459 
460  // Sort slice in next component
461  if (componentIndex < 2)
462  {
463  string oldPrefix;
464  if (debug)
465  {
466  oldPrefix = Pout.prefix();
467  Pout.prefix() = " " + oldPrefix;
468  }
469 
470  sortComponent
471  (
472  sizeTol,
473  points,
474  slice,
475  componentIndex+1,
476  mult*n_[compI], // Multiplier to apply to decomposition.
477  finalDecomp
478  );
479 
480  if (debug)
481  {
482  Pout.prefix() = oldPrefix;
483  }
484  }
485 
486  // Step to next bin.
487  leftIndex += localSize;
488  leftCoord = rightCoord;
489  }
490 }
491 
492 
493 // Sort points into bins according to one component. Recurses to next component.
494 void Foam::hierarchGeomDecomp::sortComponent
495 (
496  const label sizeTol,
497  const scalarField& weights,
498  const pointField& points,
499  const labelList& current, // slice of points to decompose
500  const direction componentIndex, // index in decompOrder_
501  const label mult, // multiplication factor for finalDecomp
502  labelList& finalDecomp
503 )
504 {
505  // Current component
506  label compI = decompOrder_[componentIndex];
507 
508  if (debug)
509  {
510  Pout<< "sortComponent : Sorting slice of size " << current.size()
511  << " in component " << compI << endl;
512  }
513 
514  // Storage for sorted component compI
515  SortableList<scalar> sortedCoord(current.size());
516 
517  forAll(current, i)
518  {
519  label pointi = current[i];
520 
521  sortedCoord[i] = points[pointi][compI];
522  }
523  sortedCoord.sort();
524 
525  label globalCurrentSize = returnReduce(current.size(), sumOp<label>());
526 
527  // Now evaluate local cumulative weights, based on the sorting.
528  // Make one bigger than the nodes.
529  scalarField sortedWeightedSizes(current.size()+1, 0);
530  calculateSortedWeightedSizes
531  (
532  current,
533  sortedCoord.indices(),
534  weights,
535  globalCurrentSize,
536  sortedWeightedSizes
537  );
538 
539  scalar minCoord = returnReduce
540  (
541  (
542  sortedCoord.size()
543  ? sortedCoord[0]
544  : great
545  ),
546  minOp<scalar>()
547  );
548 
549  scalar maxCoord = returnReduce
550  (
551  (
552  sortedCoord.size()
553  ? sortedCoord.last()
554  : -great
555  ),
556  maxOp<scalar>()
557  );
558 
559  if (debug)
560  {
561  Pout<< "sortComponent : minCoord:" << minCoord
562  << " maxCoord:" << maxCoord << endl;
563  }
564 
565  // starting index (in sortedCoord) of bin (= local)
566  label leftIndex = 0;
567  // starting value of bin (= global since coordinate)
568  scalar leftCoord = minCoord;
569 
570  // Sort bins of size n
571  for (label bin = 0; bin < n_[compI]; bin++)
572  {
573  // Now we need to determine the size of the bin (dx). This is
574  // determined by the 'pivot' values - everything to the left of this
575  // value goes in the current bin, everything to the right into the next
576  // bins.
577 
578  // Local number of elements
579  label localSize = -1; // offset from leftOffset
580 
581  // Value at right of bin (leftIndex+localSize-1)
582  scalar rightCoord = -great;
583 
584  if (bin == n_[compI]-1)
585  {
586  // Last bin. Copy all.
587  localSize = current.size()-leftIndex;
588  rightCoord = maxCoord; // note: not used anymore
589  }
590  else
591  {
592  // For the current bin (starting at leftCoord) we want a rightCoord
593  // such that the sum of all weighted sizes are
594  // globalCurrentLength/n_[compI].
595  // We have to iterate to obtain this.
596 
597  label rightIndex = current.size();
598  rightCoord = maxCoord;
599 
600  // Calculate rightIndex/rightCoord to have wanted size
601  findBinary
602  (
603  sizeTol,
604  sortedWeightedSizes,
605  sortedCoord,
606  leftIndex,
607  leftCoord,
608  maxCoord,
609  globalCurrentSize/n_[compI], // wanted size
610 
611  rightIndex,
612  rightCoord
613  );
614  localSize = rightIndex - leftIndex;
615  }
616 
617  if (debug)
618  {
619  Pout<< "For component " << compI << ", bin " << bin
620  << " copying" << endl
621  << "from " << leftCoord << " at local index "
622  << leftIndex << endl
623  << "to " << rightCoord << " localSize:"
624  << localSize << endl
625  << endl;
626  }
627 
628 
629  // Copy localSize elements starting from leftIndex.
630  labelList slice(localSize);
631 
632  forAll(slice, i)
633  {
634  label pointi = current[sortedCoord.indices()[leftIndex+i]];
635 
636  // Mark point into correct bin
637  finalDecomp[pointi] += bin*mult;
638 
639  // And collect for next sorting action
640  slice[i] = pointi;
641  }
642 
643  // Sort slice in next component
644  if (componentIndex < 2)
645  {
646  string oldPrefix;
647  if (debug)
648  {
649  oldPrefix = Pout.prefix();
650  Pout.prefix() = " " + oldPrefix;
651  }
652 
653  sortComponent
654  (
655  sizeTol,
656  weights,
657  points,
658  slice,
659  componentIndex+1,
660  mult*n_[compI], // Multiplier to apply to decomposition.
661  finalDecomp
662  );
663 
664  if (debug)
665  {
666  Pout.prefix() = oldPrefix;
667  }
668  }
669 
670  // Step to next bin.
671  leftIndex += localSize;
672  leftCoord = rightCoord;
673  }
674 }
675 
676 
678 (
679  const pointField& points
680 )
681 {
682  // construct a list for the final result
683  labelList finalDecomp(points.size(), 0);
684 
685  // Start off with every point sorted onto itself.
686  labelList slice(points.size());
687  forAll(slice, i)
688  {
689  slice[i] = i;
690  }
691 
692  pointField rotatedPoints(rotDelta_ & points);
693 
694  // Calculate tolerance of cell distribution. For large cases finding
695  // distribution to the cell exact would cause too many iterations so allow
696  // some slack.
697  label allSize = points.size();
698  reduce(allSize, sumOp<label>());
699 
700  const label sizeTol = max(1, label(1e-3*allSize/nProcessors_));
701 
702  // Sort recursive
703  sortComponent
704  (
705  sizeTol,
706  rotatedPoints,
707  slice,
708  0, // Sort first component in decompOrder.
709  1, // Offset for different x bins.
710  finalDecomp
711  );
712 
713  return finalDecomp;
714 }
715 
716 
718 (
719  const pointField& points,
720  const scalarField& weights
721 )
722 {
723  // construct a list for the final result
724  labelList finalDecomp(points.size(), 0);
725 
726  // Start off with every point sorted onto itself.
727  labelList slice(points.size());
728  forAll(slice, i)
729  {
730  slice[i] = i;
731  }
732 
733  pointField rotatedPoints(rotDelta_ & points);
734 
735  // Calculate tolerance of cell distribution. For large cases finding
736  // distribution to the cell exact would cause too many iterations so allow
737  // some slack.
738  label allSize = points.size();
739  reduce(allSize, sumOp<label>());
740 
741  const label sizeTol = max(1, label(1e-3*allSize/nProcessors_));
742 
743  // Sort recursive
744  sortComponent
745  (
746  sizeTol,
747  weights,
748  rotatedPoints,
749  slice,
750  0, // Sort first component in decompOrder.
751  1, // Offset for different x bins.
752  finalDecomp
753  );
754 
755  return finalDecomp;
756 }
757 
758 
759 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
760 
762 (
763  const dictionary& decompositionDict
764 )
765 :
766  geomDecomp(decompositionDict, typeName),
767  decompOrder_()
768 {
769  setDecompOrder();
770 }
771 
772 
773 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Inter-processor communication reduction functions.
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
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
scalar minValue
uint8_t direction
Definition: direction.H:45
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
Macros for easy insertion into run-time selection tables.
const string & prefix() const
Return the prefix of the stream.
scalar maxValue
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Geometrical domain decomposition.
Definition: geomDecomp.H:47
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
hierarchGeomDecomp(const dictionary &decompositionDict)
Construct given the decomposition dictionary.
List< label > labelList
A List of labels.
Definition: labelList.H:56
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
virtual labelList decompose(const pointField &, const scalarField &weights)
Return for every coordinate the wanted processor number.
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
#define WarningInFunction
Report a warning using Foam::Warning.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:335
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
label findLower(const ListType &, typename ListType::const_reference, const label stary, const BinaryOp &bop)
Find last element < given value in sorted list and return index,.
dimensioned< scalar > mag(const dimensioned< Type > &)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
Namespace for OpenFOAM.
IOerror FatalIOError