Merge pull request #21944 from dschwen/nonlocal_21786

Fix radial averaging system

framework/include/userobject/RadialAverage.h

@@ -19,17 +19,13 @@
 #include <set>
 #include <map>
 #include <vector>
-#include <array>
 #include <memory>
-#include <tuple>
-
-using namespace TIMPI;
 
 class ThreadedRadialAverageLoop;
 
 /**
  * Gather and communicate a full list of all quadrature points and the values of
- * a selected variable at each point. Use a KD-Tree to get the weighted spatial
+ * a selected material property at each point. Use a KD-Tree to get the weighted spatial
  * average of a material property. This code borrows heavily from
  * RadialGreensConvolution in MAGPIE. This code does not include support for
  * periodic BCs, but RadialGreensConvolution shows how that can be supported.
@@ -74,7 +70,6 @@ class RadialAverage : public ElementUserObject
   const Result & getAverage() const { return _average; }
 
 protected:
-  void insertNotLocalPointNeighbors(dof_id_type);
   void updateCommunicationLists();
 
   /// material name to get averaged
@@ -85,6 +80,9 @@ class RadialAverage : public ElementUserObject
   /// cut-off radius
   const Real _radius;
 
+  /// communication padding distance
+  const Real _padding;
+
   /// gathered data
   std::vector<QPData> _qp_data;
 
@@ -100,17 +98,14 @@ class RadialAverage : public ElementUserObject
   /// spatial index (nanoflann guarantees this to be threadsafe under read-only operations)
   std::unique_ptr<KDTreeType> _kd_tree;
 
-  /// DOF map
-  const DofMap & _dof_map;
-
-  /// PointLocator for finding topological neighbors
-  std::unique_ptr<PointLocatorBase> _point_locator;
-
   /// The data structure used to find neighboring elements give a node ID
   std::vector<std::vector<const Elem *>> _nodes_to_elem_map;
 
-  // list of direct point neighbor elements of the current processor domain
-  std::set<const Elem *> _point_neighbors;
+  /// set of nodes on the boundary of the current processor domain
+  std::set<Point> _boundary_nodes;
+
+  /// set of all _qp_data indices that are within _radius of any _boundary_nodes
+  std::set<std::size_t> _boundary_data_indices;
 
   /// QPData indices to send to the various processors
   std::vector<std::set<std::size_t>> _communication_lists;
@@ -127,10 +122,6 @@ class RadialAverage : public ElementUserObject
   friend class ThreadedRadialAverageLoop;
 };
 
-template <>
-const Point &
-PointListAdaptor<RadialAverage::QPData>::getPoint(const RadialAverage::QPData & item) const;
-
 namespace TIMPI
 {
 
@@ -142,4 +133,5 @@ class StandardType<RadialAverage::QPData> : public DataType
   StandardType(const StandardType<RadialAverage::QPData> & t);
   ~StandardType() { this->free(); }
 };
+
 } // namespace TIMPI

framework/src/userobject/RadialAverage.C

@@ -9,12 +9,12 @@
 
 #include "RadialAverage.h"
 #include "ThreadedRadialAverageLoop.h"
-#include "NonlinearSystemBase.h"
 #include "FEProblemBase.h"
 
 #include "libmesh/nanoflann.hpp"
 #include "libmesh/parallel_algebra.h"
 #include "libmesh/mesh_tools.h"
+#include "libmesh/int_range.h"
 
 #include <list>
 #include <iterator>
@@ -37,6 +37,12 @@ RadialAverage::validParams()
   params.addClassDescription("Perform a radial equal weight average of a material property");
   params.addRequiredParam<std::string>("material_name", "Name of the material to average.");
   params.addRequiredParam<Real>("radius", "Cut-off radius for the averaging");
+  params.addRangeCheckedParam<Real>(
+      "padding",
+      0.0,
+      "padding >= 0",
+      "Padding for communication. This gets added to the radius when determining which QPs to send "
+      "to other processors. Increase this gradually if inconsistent parallel results occur.");
 
   // we run this object once at the beginning of the timestep by default
   params.set<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_BEGIN;
@@ -45,7 +51,7 @@ RadialAverage::validParams()
   params.addRelationshipManager("ElementPointNeighborLayers",
                                 Moose::RelationshipManagerType::GEOMETRIC |
                                     Moose::RelationshipManagerType::ALGEBRAIC);
-
+  params.addParamNamesToGroup("padding", "Advanced");
   return params;
 }
 
@@ -54,7 +60,7 @@ RadialAverage::RadialAverage(const InputParameters & parameters)
     _averaged_material_name(getParam<std::string>("material_name")),
     _averaged_material(getMaterialProperty<Real>(_averaged_material_name)),
     _radius(getParam<Real>("radius")),
-    _dof_map(_fe_problem.getNonlinearSystemBase().dofMap()),
+    _padding(getParam<Real>("padding")),
     _update_communication_lists(false),
     _my_pid(processor_id()),
     _perf_meshchanged(registerTimedSection("meshChanged", 3)),
@@ -82,7 +88,7 @@ RadialAverage::execute()
   auto id = _current_elem->id();
 
   // collect all QP data
-  for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
+  for (const auto qp : make_range(_qrule->n_points()))
     _qp_data.emplace_back(_q_point[qp], id, qp, _JxW[qp] * _coord[qp], _averaged_material[qp]);
 
   // make sure the result map entry for the current element is sized correctly
@@ -126,17 +132,7 @@ RadialAverage::finalize()
     Parallel::MessageTag send_tag = _communicator.get_unique_tag(4711);
     std::vector<QPData> receive;
 
-    const auto item_type = StandardType<QPData>(&(_qp_data[0]));
-
-#if 0
-    // output local qp locations
-    // _console << name() << ' ' << receive.size() << '\n' << name() << std::flush;
-    for (auto item : _qp_data)
-      _console << name() << ' ' << _my_pid << ' '
-               << item._q_point(0) << ' '
-               << item._q_point(1) << ' '
-               << item._q_point(2) << std::flush;
-#endif
+    const auto item_type = TIMPI::StandardType<QPData>(&(_qp_data[0]));
 
     // fill buffer and send structures
     for (auto i = beginIndex(non_zero_comm); i < non_zero_comm.size(); ++i)
@@ -147,17 +143,7 @@ RadialAverage::finalize()
       // fill send buffer for transfer to pid
       send[i].reserve(list.size());
       for (const auto & item : list)
-      {
         send[i].push_back(_qp_data[item]);
-#if 0
-        // output sent qp locations
-        _console << name() << ' '
-                 << _qp_data[item]._q_point(0) << ' '
-                 << _qp_data[item]._q_point(1) << ' '
-                 << _qp_data[item]._q_point(2) << ' '
-                 << pid << std::flush;
-#endif
-      }
 
       // issue non-blocking send
       _communicator.send(pid, send[i], send_requests[i], send_tag);
@@ -175,16 +161,6 @@ RadialAverage::finalize()
       receive.resize(message_size);
       _communicator.receive(source_pid, receive, send_tag);
 
-#if 0
-      // output received qp locations
-      // _console << name() << ' ' << receive.size() << '\n' << name() << std::flush;
-      for (auto item : receive)
-        _console << name() << ' ' << source_pid << ' '
-                 << item._q_point(0) << ' '
-                 << item._q_point(1) << ' '
-                 << item._q_point(2) << std::flush;
-#endif
-
       // append communicated data
       _qp_data.insert(_qp_data.end(), receive.begin(), receive.end());
     }
@@ -221,16 +197,6 @@ RadialAverage::threadJoin(const UserObject & y)
   _average.insert(uo._average.begin(), uo._average.end());
 }
 
-void
-RadialAverage::insertNotLocalPointNeighbors(dof_id_type node)
-{
-  mooseAssert(!_nodes_to_elem_map[node].empty(), "Node not found in _nodes_to_elem_map");
-
-  for (const auto * elem : _nodes_to_elem_map[node])
-    if (elem->processor_id() != _my_pid && elem->active())
-      _point_neighbors.insert(elem);
-}
-
 void
 RadialAverage::meshChanged()
 {
@@ -239,39 +205,31 @@ RadialAverage::meshChanged()
   // get underlying libMesh mesh
   auto & mesh = _mesh.getMesh();
 
-  // get a fresh point locator
-  _point_locator = _mesh.getPointLocator();
-
   // Build a new node to element map
   _nodes_to_elem_map.clear();
   MeshTools::build_nodes_to_elem_map(_mesh.getMesh(), _nodes_to_elem_map);
 
-  // clear point neighbor set
-  _point_neighbors.clear();
+  // clear procesor boundary nodes set
+  _boundary_nodes.clear();
 
-  // my processor id
-  const auto my_pid = processor_id();
-
-  // iterate over active local elements
+  //
+  // iterate over active local elements and store all processor boundary node locations
+  //
   const auto end = mesh.active_local_elements_end();
   for (auto it = mesh.active_local_elements_begin(); it != end; ++it)
-    // find a face that faces either a boundary (nullptr) or a different processor
-    for (unsigned int s = 0; s < (*it)->n_sides(); ++s)
+    // find faces at processor boundaries
+    for (const auto s : make_range((*it)->n_sides()))
     {
       const auto * neighbor = (*it)->neighbor_ptr(s);
-      if (neighbor)
-      {
-        if (neighbor->processor_id() != my_pid)
-        {
-          // add all face node touching elements directly
-          for (unsigned int n = 0; n < (*it)->n_nodes(); ++n)
-            if ((*it)->is_node_on_side(n, s))
-              insertNotLocalPointNeighbors((*it)->node_id(n));
-        }
-      }
+      if (neighbor && neighbor->processor_id() != _my_pid)
+        // add all nodes on the processor boundary
+        for (const auto n : make_range((*it)->n_nodes()))
+          if ((*it)->is_node_on_side(n, s))
+            _boundary_nodes.insert((*it)->node_ref(n));
+
+      // request communication list update
+      _update_communication_lists = true;
     }
-  // request communication list update
-  _update_communication_lists = true;
 }
 
 void
@@ -294,17 +252,42 @@ RadialAverage::updateCommunicationLists()
   std::vector<std::pair<std::size_t, Real>> ret_matches;
   nanoflann::SearchParams search_params;
 
-  // iterate over non periodic point neighbor elements
-  for (auto elem : _point_neighbors)
+  // iterate over all boundary nodes and collect all boundary-near data points
+  _boundary_data_indices.clear();
+  for (const auto & bn : _boundary_nodes)
   {
     ret_matches.clear();
-    Point centroid = elem->vertex_average();
-    const Real radius2 = _radius + elem->hmax() / 2.0;
-    kd_tree->radiusSearch(&(centroid(0)), radius2, ret_matches, search_params);
+    kd_tree->radiusSearch(
+        &(bn(0)), Utility::pow<2>(_radius + _padding), ret_matches, search_params);
     for (auto & match : ret_matches)
-      _communication_lists[elem->processor_id()].insert(match.first);
+      _boundary_data_indices.insert(match.first);
   }
 
+  // gather all processor bounding boxes (communicate as pairs)
+  std::vector<std::pair<Point, Point>> pps(n_processors());
+  const auto mybb = _mesh.getInflatedProcessorBoundingBox(0);
+  std::pair<Point, Point> mypp = mybb;
+  _communicator.allgather(mypp, pps);
+
+  // inflate all processor bounding boxes by radius (no padding)
+  const auto rpoint = Point(1, 1, 1) * _radius;
+  std::vector<BoundingBox> bbs;
+  for (const auto & pp : pps)
+    bbs.emplace_back(pp.first - rpoint, pp.second + rpoint);
+
+  // get candidate processors (overlapping bounding boxes)
+  std::vector<processor_id_type> candidate_procs;
+  for (const auto pid : index_range(bbs))
+    if (pid != _my_pid && bbs[pid].intersects(mypp))
+      candidate_procs.push_back(pid);
+
+  // go over all boundary data items and send them to the proc they overlap with
+  for (const auto i : _boundary_data_indices)
+    for (const auto pid : candidate_procs)
+      if (bbs[pid].contains_point(_qp_data[i]._q_point))
+        _communication_lists[pid].insert(i);
+
+  // done
   _update_communication_lists = false;
 }
 
@@ -356,7 +339,6 @@ StandardType<RadialAverage::QPData>::StandardType(const RadialAverage::QPData *
 
 #endif // LIBMESH_HAVE_MPI
 }
-} // namespace TIMPI
 
 StandardType<RadialAverage::QPData>::StandardType(const StandardType<RadialAverage::QPData> & t)
   : DataType(t._datatype)
@@ -365,3 +347,5 @@ StandardType<RadialAverage::QPData>::StandardType(const StandardType<RadialAvera
   libmesh_call_mpi(MPI_Type_dup(t._datatype, &_datatype));
 #endif
 }
+
+} // namespace TIMPI

framework/src/userobject/ThreadedRadialAverageLoop.C

@@ -50,7 +50,7 @@ ThreadedRadialAverageLoop::operator()(const QPDataRange & qpdata_range)
 
     ret_matches.clear();
     std::size_t n_result =
-        kd_tree->radiusSearch(&(local_qp._q_point(0)), radius, ret_matches, search_params);
+        kd_tree->radiusSearch(&(local_qp._q_point(0)), radius * radius, ret_matches, search_params);
     Real total_vol = 0.0;
     for (std::size_t j = 0; j < n_result; ++j)
     {

modules/tensor_mechanics/test/tests/scalar_material_damage/ad_nonlocal_scalar_damage.i

@@ -63,7 +63,7 @@
     material_name = local_damage_reg
     execute_on = "INITIAL timestep_end"
     block = 0
-    radius = 0.3
+    radius = 0.55
   []
 []
 

modules/tensor_mechanics/test/tests/scalar_material_damage/nonlocal_scalar_damage.i

@@ -22,7 +22,6 @@
   []
 []
 
-
 [BCs]
   [symmy]
     type = DirichletBC
@@ -63,7 +62,7 @@
     material_name = local_damage
     execute_on = "INITIAL timestep_end"
     block = 0
-    radius = 0.3
+    radius = 0.55
   []
 []
 

test/tests/userobjects/radial_average/tests

@@ -1,13 +1,13 @@
 [Tests]
-    issues = '#21786'
-    design = '/RadialAverage.md'
+  issues = '#21786'
+  design = '/RadialAverage.md'
 
-    [./radial_average]
-      type = 'Exodiff'
-      input = 'time_changing_test.i'
-      exodiff = 'time_changing_test_out.e'
-      recover = false
-      requirement = 'The system shall will compute a radial average of a material that changes over time.'
-      max_parallel = 2 #See comment on issue #21786
-    [../]
+  [radial_average]
+    type = 'Exodiff'
+    input = 'time_changing_test.i'
+    exodiff = 'time_changing_test_out.e'
+    recover = false
+    requirement = 'The system shall will compute a radial average of a material that changes over time.'
+    abs_zero = 1e-10
   []
+[]

test/tests/userobjects/radial_average/time_changing_test.i

@@ -51,7 +51,7 @@
     material_name = local
     execute_on = "INITIAL timestep_end"
     block = 0
-    radius = 0.25
+    radius = 0.3
   []
 []