Merge pull request #24840 from GiudGiud/PR_extreme_fix

Fix parallel logic for layered extremum
idaholab · Jun 28, 2023 · 467ef88 · 467ef88
2 parents 21e49dd + 3327884
commit 467ef88
Show file tree

Hide file tree

Showing 3 changed files with 41 additions and 15 deletions.
diff --git a/framework/include/userobjects/LayeredBase.h b/framework/include/userobjects/LayeredBase.h
@@ -135,22 +135,22 @@ class LayeredBase : private Restartable
   Real _direction_min;
   Real _direction_max;
 
-private:
   /// Value of the integral for each layer
   std::vector<Real> & _layer_values;
 
   /// Whether or not each layer has had any value summed into it
   std::vector<int> & _layer_has_value;
 
-  /// Subproblem for the child object
-  SubProblem & _layered_base_subproblem;
-
   /// Whether the values are cumulative over the layers
   bool _cumulative;
 
   /// Whether the cumulative values should be summed in the positive or negative direction
   const bool _positive_cumulative_direction;
 
+private:
+  /// Subproblem for the child object
+  SubProblem & _layered_base_subproblem;
+
   /// List of SubdomainIDs, if given
   std::vector<SubdomainID> _layer_bounding_blocks;
 

diff --git a/framework/src/userobjects/LayeredBase.C b/framework/src/userobjects/LayeredBase.C
@@ -89,9 +89,9 @@ LayeredBase::LayeredBase(const InputParameters & parameters)
     _using_displaced_mesh(_layered_base_params.get<bool>("use_displaced_mesh")),
     _layer_values(declareRestartableData<std::vector<Real>>("layer_values")),
     _layer_has_value(declareRestartableData<std::vector<int>>("layer_has_value")),
-    _layered_base_subproblem(*parameters.getCheckedPointerParam<SubProblem *>("_subproblem")),
     _cumulative(parameters.get<bool>("cumulative")),
     _positive_cumulative_direction(parameters.get<bool>("positive_cumulative_direction")),
+    _layered_base_subproblem(*parameters.getCheckedPointerParam<SubProblem *>("_subproblem")),
     _layer_bounding_blocks(),
     _has_direction_max_min(false)
 {
@@ -250,7 +250,7 @@ LayeredBase::integralValue(Point p) const
 
       if (higher_layer != -1)
       {
-        for (unsigned int i = 0; i < _average_radius; i++)
+        for (const auto i : make_range(_average_radius))
         {
           int current_layer = higher_layer + i;
 
@@ -267,7 +267,7 @@ LayeredBase::integralValue(Point p) const
 
       if (lower_layer != -1)
       {
-        for (unsigned int i = 0; i < _average_radius; i++)
+        for (const auto i : make_range(_average_radius))
         {
           int current_layer = lower_layer - i;
 
@@ -303,7 +303,7 @@ LayeredBase::initialize()
   if (_using_displaced_mesh)
     getBounds();
 
-  for (unsigned int i = 0; i < _layer_values.size(); i++)
+  for (const auto i : index_range(_layer_values))
   {
     _layer_values[i] = 0.0;
     _layer_has_value[i] = false;
@@ -321,7 +321,7 @@ LayeredBase::finalize()
     Real value = 0;
 
     if (_positive_cumulative_direction)
-      for (unsigned i = 0; i < _num_layers; i++)
+      for (const auto i : make_range(_num_layers))
       {
         value += getLayerValue(i);
         setLayerValue(i, value);
@@ -339,7 +339,7 @@ void
 LayeredBase::threadJoin(const UserObject & y)
 {
   const LayeredBase & lb = dynamic_cast<const LayeredBase &>(y);
-  for (unsigned int i = 0; i < _layer_values.size(); i++)
+  for (const auto i : index_range(_layer_values))
     if (lb.layerHasValue(i))
       setLayerValue(i, getLayerValue(i) + lb._layer_values[i]);
 }
@@ -393,12 +393,12 @@ LayeredBase::computeLayerCenters()
   {
     Real dx = (_direction_max - _direction_min) / _num_layers;
 
-    for (unsigned int i = 0; i < _num_layers; ++i)
+    for (const auto i : make_range(_num_layers))
       _layer_centers[i] = (i + 0.5) * dx;
   }
   else
   {
-    for (unsigned int i = 0; i < _num_layers; ++i)
+    for (const auto i : make_range(_num_layers))
       _layer_centers[i] = 0.5 * (_layer_bounds[i + 1] + _layer_bounds[i]);
   }
 }

diff --git a/framework/src/userobjects/LayeredExtremumMaterialProperty.C b/framework/src/userobjects/LayeredExtremumMaterialProperty.C
@@ -64,14 +64,40 @@ LayeredExtremumMaterialProperty::extreme_value(const Real a, const Real b) const
 void
 LayeredExtremumMaterialProperty::finalize()
 {
-  LayeredBase::finalize();
+  if (_type == MIN)
+    comm().min(_layer_values);
+  else
+    comm().max(_layer_values);
+  comm().max(_layer_has_value);
+
+  if (_cumulative)
+  {
+    Real value =
+        _type == MIN ? std::numeric_limits<Real>::max() : -std::numeric_limits<Real>::max();
+
+    if (_positive_cumulative_direction)
+      for (unsigned i = 0; i < _num_layers; i++)
+      {
+        value = extreme_value(value, getLayerValue(i));
+        setLayerValue(i, value);
+      }
+    else
+      for (int i = _num_layers - 1; i >= 0; i--)
+      {
+        value = extreme_value(value, getLayerValue(i));
+        setLayerValue(i, value);
+      }
+  }
 }
 
 void
 LayeredExtremumMaterialProperty::threadJoin(const UserObject & y)
 {
-  ElementExtremeMaterialProperty::threadJoin(y);
-  LayeredBase::threadJoin(y);
+  const LayeredExtremumMaterialProperty & lb =
+      static_cast<const LayeredExtremumMaterialProperty &>(y);
+  for (const auto i : make_range(_num_layers))
+    if (lb.layerHasValue(i))
+      setLayerValue(i, extreme_value(getLayerValue(i), lb.getLayerValue(i)));
 }
 
 const std::vector<Point>