More robust communication (idaholab#21786)

framework/src/userobject/RadialAverage.C

@@ -39,7 +39,7 @@ RadialAverage::validParams()
   params.addRequiredParam<Real>("radius", "Cut-off radius for the averaging");
   params.addParam<Real>(
       "padding",
-      0.1,
+      0.0,
       "Padding for communication. This gets added to the radius when determining which QPs to send "
       "to other processors. It should be of the order of the maximum element size.");
 
@@ -199,6 +199,15 @@ RadialAverage::finalize()
     Parallel::wait(send_requests);
   }
 
+  {
+    std::ofstream dbg1(std::to_string(_my_pid) + "comm.dat");
+    for (const auto & qd : _qp_data)
+    {
+      const auto p = qd._q_point;
+      dbg1 << p(0) << ' ' << p(1) << ' ' << p(2) << std::endl;
+    }
+  }
+
   // build KD-Tree using data we just received
   const unsigned int max_leaf_size = 20; // slightly affects runtime
   auto point_list = PointListAdaptor<QPData>(_qp_data.begin(), _qp_data.end());
@@ -286,7 +295,8 @@ RadialAverage::updateCommunicationLists()
   for (const auto & bn : _boundary_nodes)
   {
     ret_matches.clear();
-    kd_tree->radiusSearch(&(bn(0)), _radius + _padding, ret_matches, search_params);
+    kd_tree->radiusSearch(
+        &(bn(0)), Utility::pow<2>(_radius + _padding), ret_matches, search_params);
     for (auto & match : ret_matches)
       _boundary_data_indices.insert(match.first);
   }
@@ -312,14 +322,11 @@ RadialAverage::updateCommunicationLists()
   std::pair<Point, Point> mypp = mybb;
   _communicator.allgather(mypp, pps);
 
-  // inflate all processor bounding boxes by radius
-  const Point rpoint(_radius + _padding, _radius + _padding, _radius + _padding);
+  // inflate all processor bounding boxes by radius + padding
+  const auto rpoint = Point(1, 1, 1) * (_radius + _padding);
   std::vector<BoundingBox> bbs;
   for (const auto & pp : pps)
-  {
-    std::cout << pp.first << ", " << pp.second << std::endl;
     bbs.emplace_back(pp.first - rpoint, pp.second + rpoint);
-  }
 
   // get candidate processors (overlapping bounding boxes)
   std::vector<processor_id_type> candidate_procs;
@@ -409,55 +416,4 @@ StandardType<RadialAverage::QPData>::StandardType(const StandardType<RadialAvera
 #endif
 }
 
-// StandardType<BoundingBox>::StandardType(const BoundingBox * example)
-// {
-//   // We need an example for MPI_Address to use
-//   static const BoundingBox p;
-//   if (!example)
-//     example = &p;
-
-// #ifdef LIBMESH_HAVE_MPI
-
-//   // Get the sub-data-types, and make sure they live long enough
-//   // to construct the derived type
-//   StandardType<Point> d1(&example->min());
-//   StandardType<Point> d2(&example->max());
-
-//   MPI_Datatype types[] = {
-//       (data_type)d1, (data_type)d2};
-//   int blocklengths[] = {1, 1};
-//   MPI_Aint displs[5], start;
-
-//   libmesh_call_mpi(MPI_Get_address(const_cast<BoundingBox *>(example), &start));
-//   libmesh_call_mpi(MPI_Get_address(const_cast<Point *>(&example->_q_point), &displs[0]));
-//   libmesh_call_mpi(MPI_Get_address(const_cast<dof_id_type *>(&example->_elem_id), &displs[1]));
-//   libmesh_call_mpi(MPI_Get_address(const_cast<short *>(&example->_qp), &displs[2]));
-//   libmesh_call_mpi(MPI_Get_address(const_cast<Real *>(&example->_volume), &displs[3]));
-//   libmesh_call_mpi(MPI_Get_address(const_cast<Real *>(&example->_value), &displs[4]));
-
-//   for (std::size_t i = 0; i < 5; ++i)
-//     displs[i] -= start;
-
-//   // create a prototype structure
-//   MPI_Datatype tmptype;
-//   libmesh_call_mpi(MPI_Type_create_struct(5, blocklengths, displs, types, &tmptype));
-//   libmesh_call_mpi(MPI_Type_commit(&tmptype));
-
-//   // resize the structure type to account for padding, if any
-//   libmesh_call_mpi(MPI_Type_create_resized(tmptype, 0, sizeof(RadialAverage::QPData),
-//   &_datatype)); libmesh_call_mpi(MPI_Type_free(&tmptype));
-
-//   this->commit();
-
-// #endif // LIBMESH_HAVE_MPI
-// }
-
-// StandardType<BoundingBox>::StandardType(const StandardType<BoundingBox> & t)
-//   : DataType(t._datatype)
-// {
-// #ifdef LIBMESH_HAVE_MPI
-//   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)
     {

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.2
+    radius = 0.3
   []
 []