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