Reinterpolate to all elements when interpolator receives new points

This is believed to be what was causing a lot of deadlocks in BBH
simulations.

src/ParallelAlgorithms/ApparentHorizonFinder/Callbacks/FindApparentHorizon.hpp

@@ -267,9 +267,10 @@ struct FindApparentHorizon
         auto& interpolation_target = Parallel::get_parallel_component<
             intrp::InterpolationTarget<Metavariables, InterpolationTargetTag>>(
             *cache);
+        // The iteration of these new coords is the fast flow iteration
         Parallel::simple_action<
             Actions::SendPointsToInterpolator<InterpolationTargetTag>>(
-            interpolation_target, temporal_ids.front());
+            interpolation_target, temporal_ids.front(), info.iteration);
         // We return false because we don't want this iteration to clean
         // up the volume data, since we are using it for the next iteration
         // (i.e. the simple_action that we just called).

src/ParallelAlgorithms/Interpolation/Actions/InterpolatorReceivePoints.hpp

@@ -45,6 +45,14 @@ namespace Actions {
 /// After receiving the points, interpolates volume data onto them
 /// if it already has all the volume data.
 ///
+/// The `iteration` parameter is used to order receives of
+/// `block_logical_coords`. Because of the asynchronous nature of communication,
+/// it is possible that a more recent set of points arrives before an older set.
+/// It is assumed that if a more recent set arrives, then the old set is no
+/// longer needed. This `iteration` parameter tags each communication as "more
+/// recent" or "older" so if we receive an older set of points after a more
+/// recent set, we don't overwrite the more recent set.
+///
 /// Uses:
 /// - Databox:
 ///   - `Tags::NumberOfElements`
@@ -69,9 +77,10 @@ struct ReceivePoints {
       std::vector<std::optional<
           IdPair<domain::BlockId,
                  tnsr::I<double, VolumeDim, typename ::Frame::BlockLogical>>>>&&
-          block_logical_coords) {
+          block_logical_coords,
+      const size_t iteration = 0_st) {
     db::mutate<intrp::Tags::InterpolatedVarsHolders<Metavariables>>(
-        [&temporal_id, &block_logical_coords](
+        [&temporal_id, &block_logical_coords, &iteration](
             const gsl::not_null<typename intrp::Tags::InterpolatedVarsHolders<
                 Metavariables>::type*>
                 vars_holders) {
@@ -80,12 +89,34 @@ struct ReceivePoints {
                                          Metavariables>>(*vars_holders)
                   .infos;
 
-          // Add the target interpolation points at this temporal_id.
-          vars_infos.emplace(std::make_pair(
-              temporal_id,
-              intrp::Vars::Info<VolumeDim, typename InterpolationTargetTag::
-                                               vars_to_interpolate_to_target>{
-                  std::move(block_logical_coords)}));
+          // Add the new target interpolation points at this temporal_id. There
+          // are two conditions that allow us to overwrite the current target
+          // points. Either
+          //
+          //  1. There are no current target points at the temporal_id, OR
+          //  2. There are target points already at this temporal_id, but the
+          //     iteration of the new target points is greater than the
+          //     iteration of the current target points.
+          //
+          // If we already have target points and the iteration of the new
+          // points is less than or equal to the iteration of the current target
+          // points, then we ignore the new points. The new points are outdated
+          // and we definitely didn't have any of the new target points in our
+          // element by the fact that we have already received the next
+          // iteration of points.
+          //
+          // Whenever we overwrite the target points, we also empty the
+          // `interpolation_is_done_for_these_elements` (by virtue of a default
+          // constructed `intrp::Vars::Info`) so that we always check every
+          // element for this new set of target points.
+          if (vars_infos.count(temporal_id) == 0 or
+              vars_infos.at(temporal_id).iteration < iteration) {
+            vars_infos.insert_or_assign(
+                temporal_id,
+                intrp::Vars::Info<VolumeDim, typename InterpolationTargetTag::
+                                                 vars_to_interpolate_to_target>{
+                    std::move(block_logical_coords), iteration});
+          }
         },
         make_not_null(&box));
 

src/ParallelAlgorithms/Interpolation/Actions/SendPointsToInterpolator.hpp

@@ -3,6 +3,7 @@
 
 #pragma once
 
+#include <cstddef>
 #include <utility>
 
 #include "DataStructures/DataBox/DataBox.hpp"
@@ -18,6 +19,11 @@ namespace Actions {
 /// \brief Sets up points on an `InterpolationTarget` at a new `temporal_id`
 /// and sends these points to an `Interpolator`.
 ///
+/// The `iteration` parameter tags each set of points so the `Interpolator`
+/// knows which are newer points and which are older points.
+///
+/// \see `intrp::Actions::ReceivePoints`
+///
 /// Uses:
 /// - DataBox:
 ///   - `domain::Tags::Domain<3>`
@@ -38,15 +44,16 @@ struct SendPointsToInterpolator {
   static void apply(db::DataBox<DbTags>& box,
                     Parallel::GlobalCache<Metavariables>& cache,
                     const ArrayIndex& /*array_index*/,
-                    const TemporalId& temporal_id) {
+                    const TemporalId& temporal_id,
+                    const size_t iteration = 0_st) {
     auto coords = InterpolationTarget_detail::block_logical_coords<
         InterpolationTargetTag>(box, cache, temporal_id);
     InterpolationTarget_detail::set_up_interpolation<InterpolationTargetTag>(
         make_not_null(&box), temporal_id, coords);
     auto& receiver_proxy =
         Parallel::get_parallel_component<Interpolator<Metavariables>>(cache);
     Parallel::simple_action<Actions::ReceivePoints<InterpolationTargetTag>>(
-        receiver_proxy, temporal_id, std::move(coords));
+        receiver_proxy, temporal_id, std::move(coords), iteration);
   }
 };
 

src/ParallelAlgorithms/Interpolation/InterpolatedVars.hpp

@@ -45,6 +45,10 @@ struct Info {
       IdPair<domain::BlockId,
              tnsr::I<double, VolumeDim, typename ::Frame::BlockLogical>>>>
       block_coord_holders;
+  /// If a target needs to send points in a specific order, it should also send
+  /// along which iteration the `block_coord_holders` are for. That way they can
+  /// be properly ordered in the Interpolator.
+  size_t iteration{0_st};
   /// `vars` holds the interpolated `Variables` on some subset of the
   /// points in `block_coord_holders`.  The grid points inside vars
   /// are indexed according to `global_offsets` below.  The size of

tests/Unit/Helpers/ParallelAlgorithms/Interpolation/InterpolationTargetTestHelpers.hpp

@@ -101,9 +101,10 @@ struct MockReceivePoints {
       std::vector<std::optional<
           IdPair<domain::BlockId,
                  tnsr::I<double, VolumeDim, typename Frame::BlockLogical>>>>&&
-          block_coord_holders) {
+          block_coord_holders,
+      const size_t iteration = 0_st) {
     db::mutate<intrp::Tags::InterpolatedVarsHolders<Metavariables>>(
-        [&temporal_id, &block_coord_holders](
+        [&temporal_id, &block_coord_holders, &iteration](
             const gsl::not_null<typename intrp::Tags::InterpolatedVarsHolders<
                 Metavariables>::type*>
                 vars_holders) {
@@ -117,7 +118,7 @@ struct MockReceivePoints {
               temporal_id,
               intrp::Vars::Info<VolumeDim, typename InterpolationTargetTag::
                                                vars_to_interpolate_to_target>{
-                  std::move(block_coord_holders)}));
+                  std::move(block_coord_holders), iteration}));
         },
         make_not_null(&box));
   }
@@ -131,12 +132,11 @@ struct mock_interpolator {
   using phase_dependent_action_list = tmpl::list<
       Parallel::PhaseActions<
           Parallel::Phase::Initialization,
-          tmpl::list<
-              intrp::Actions::InitializeInterpolator<
-                  intrp::Tags::VolumeVarsInfo<
-                      Metavariables, typename Metavariables::
-                                         InterpolationTargetA::temporal_id>,
-                  intrp::Tags::InterpolatedVarsHolders<Metavariables>>>>,
+          tmpl::list<intrp::Actions::InitializeInterpolator<
+              intrp::Tags::VolumeVarsInfo<
+                  Metavariables,
+                  typename Metavariables::InterpolationTargetA::temporal_id>,
+              intrp::Tags::InterpolatedVarsHolders<Metavariables>>>>,
       Parallel::PhaseActions<Parallel::Phase::Testing, tmpl::list<>>>;
 
   using component_being_mocked = intrp::Interpolator<Metavariables>;
@@ -235,9 +235,10 @@ void test_interpolation_target(
   const auto& info = vars_infos.at(temporal_id);
   const auto& block_coord_holders = info.block_coord_holders;
 
-  // Check number of points
+  // Check number of points and iteration
   const size_t number_of_points = expected_block_coord_holders.size();
   CHECK(block_coord_holders.size() == number_of_points);
+  CHECK(info.iteration == 0_st);
 
   for (size_t i = 0; i < number_of_points; ++i) {
     CHECK(block_coord_holders[i].value().id ==
@@ -255,10 +256,12 @@ void test_interpolation_target(
   ActionTesting::invoke_queued_simple_action<interp_component>(
       make_not_null(&runner), 0);
 
-  // Should be two entries in the vars_infos
+  // Should be two entries in the vars_infos. Second should also have iteration
+  // 0
   CHECK(vars_infos.size() == 2);
-  const auto& new_block_coord_holders =
-      vars_infos.at(new_temporal_id).block_coord_holders;
+  const auto& new_info = vars_infos.at(new_temporal_id);
+  const auto& new_block_coord_holders = new_info.block_coord_holders;
+  CHECK(new_info.iteration == 0_st);
   for (size_t i = 0; i < number_of_points; ++i) {
     CHECK(new_block_coord_holders[i].value().id ==
           expected_block_coord_holders[i].value().id);