Merge pull request #5864 from wthrowe/TimeStepId_ordering

Generalize LTS communication handling

src/Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp

@@ -210,19 +210,18 @@ bool receive_boundary_data_local_time_stepping(
           Dim>>(*inboxes);
 
   const auto needed_time = [&box]() {
-    const auto& local_next_temporal_id =
-        db::get<::Tags::Next<::Tags::TimeStepId>>(*box);
+    const LtsTimeStepper& time_stepper =
+        db::get<::Tags::TimeStepper<LtsTimeStepper>>(*box);
     if constexpr (DenseOutput) {
       const auto& dense_output_time = db::get<::Tags::Time>(*box);
-      return
-          [&dense_output_time, &local_next_temporal_id](const TimeStepId& id) {
-            return evolution_less<double>{
-                local_next_temporal_id.time_runs_forward()}(
-                id.step_time().value(), dense_output_time);
-          };
+      return [&dense_output_time, &time_stepper](const TimeStepId& id) {
+        return time_stepper.neighbor_data_required(dense_output_time, id);
+      };
     } else {
-      return [&local_next_temporal_id](const TimeStepId& id) {
-        return id < local_next_temporal_id;
+      const auto& next_temporal_id =
+          db::get<::Tags::Next<::Tags::TimeStepId>>(*box);
+      return [&next_temporal_id, &time_stepper](const TimeStepId& id) {
+        return time_stepper.neighbor_data_required(next_temporal_id, id);
       };
     }
   }();
@@ -243,24 +242,30 @@ bool receive_boundary_data_local_time_stepping(
               boundary_data_history,
           const gsl::not_null<
               std::unordered_map<Key, TimeStepId, boost::hash<Key>>*>
-              mortar_next_time_step_id,
+              mortar_next_time_step_ids,
           const gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*>
               neighbor_mesh,
           const Element<Dim>& element) {
         // Remove neighbor meshes for neighbors that don't exist anymore
         domain::remove_nonexistent_neighbors(neighbor_mesh, element);
 
         // Move received boundary data into boundary history.
-        for (auto received_data = inbox.begin();
-             received_data != inbox.end() and needed_time(received_data->first);
-             received_data = inbox.erase(received_data)) {
-          const auto& receive_temporal_id = received_data->first;
-          // Loop over all mortars for which we received data at this time
-          for (auto received_mortar_data = received_data->second.begin();
-               received_mortar_data != received_data->second.end();
-               received_mortar_data =
-                   received_data->second.erase(received_mortar_data)) {
-            const auto& mortar_id = received_mortar_data->first;
+        for (auto& [mortar_id, mortar_next_time_step_id] :
+             *mortar_next_time_step_ids) {
+          if (mortar_id.id == ElementId<Dim>::external_boundary_id()) {
+            continue;
+          }
+          while (needed_time(mortar_next_time_step_id)) {
+            const auto time_entry = inbox.find(mortar_next_time_step_id);
+            if (time_entry == inbox.end()) {
+              return false;
+            }
+            const auto received_mortar_data =
+                time_entry->second.find(mortar_id);
+            if (received_mortar_data == time_entry->second.end()) {
+              return false;
+            }
+
             MortarData<Dim> neighbor_mortar_data{};
             // Insert:
             // - the current TimeStepId of the neighbor
@@ -269,53 +274,33 @@ bool receive_boundary_data_local_time_stepping(
             ASSERT(std::get<3>(received_mortar_data->second).has_value(),
                    "Did not receive boundary correction data from the "
                    "neighbor\nMortarId: "
-                       << mortar_id << "\nTimeStepId: " << receive_temporal_id);
-            ASSERT(
-                mortar_next_time_step_id->at(mortar_id) >= receive_temporal_id,
-                "Expected to receive mortar data on mortar "
-                    << mortar_id << " at time "
-                    << mortar_next_time_step_id->at(mortar_id)
-                    << " but actually received at time "
-                    << receive_temporal_id);
-            if (mortar_next_time_step_id->at(mortar_id) !=
-                receive_temporal_id) {
-              // We've received messages from our neighbor
-              // out-of-order.  They are always sent in-order, but
-              // messages are not guaranteed to be received in the
-              // order they were sent.
-              return false;
-            }
+                       << mortar_id
+                       << "\nTimeStepId: " << mortar_next_time_step_id);
             neighbor_mesh->insert_or_assign(
                 mortar_id, std::get<0>(received_mortar_data->second));
-            mortar_next_time_step_id->at(mortar_id) =
-                std::get<4>(received_mortar_data->second);
             neighbor_mortar_data.insert_neighbor_mortar_data(
-                receive_temporal_id, std::get<1>(received_mortar_data->second),
+                mortar_next_time_step_id,
+                std::get<1>(received_mortar_data->second),
                 std::move(*std::get<3>(received_mortar_data->second)));
             // We don't yet communicate the integration order, because
             // we don't have any variable-order methods.  The
             // fixed-order methods ignore the field.
             boundary_data_history->at(mortar_id).remote().insert(
-                receive_temporal_id, std::numeric_limits<size_t>::max(),
+                mortar_next_time_step_id, std::numeric_limits<size_t>::max(),
                 std::move(neighbor_mortar_data));
+            mortar_next_time_step_id =
+                std::get<4>(received_mortar_data->second);
+            time_entry->second.erase(received_mortar_data);
+            if (time_entry->second.empty()) {
+              inbox.erase(time_entry);
+            }
           }
         }
         return true;
       },
       box, db::get<::domain::Tags::Element<Dim>>(*box));
 
-  if (not have_all_intermediate_messages) {
-    return false;
-  }
-
-  return alg::all_of(
-      db::get<evolution::dg::Tags::MortarNextTemporalId<Dim>>(*box),
-      [&needed_time](
-          const std::pair<Key, TimeStepId>& mortar_id_and_next_temporal_id) {
-        return mortar_id_and_next_temporal_id.first.id ==
-                   ElementId<Dim>::external_boundary_id() or
-               not needed_time(mortar_id_and_next_temporal_id.second);
-      });
+  return have_all_intermediate_messages;
 }
 
 /// Apply corrections from boundary communication.

src/Time/TimeStepId.cpp

@@ -122,9 +122,16 @@ bool operator<(const TimeStepId& a, const TimeStepId& b) {
     return evolution_less<Time>{a.time_runs_forward()}(a.step_time(),
                                                        b.step_time());
   }
-  ASSERT(a.substep() == 0 or b.substep() == 0 or a.step_size() == b.step_size(),
-         "Meaning of ordering for LTS substeps is unclear.");
-  return a.substep() < b.substep();
+  if (a.substep() != b.substep()) {
+    return a.substep() < b.substep();
+  }
+  if (a.substep() == 0) {
+    // Objects are equal.
+    return false;
+  }
+  // Arbitrary, but need a total ordering of TimeStepId to use it as a
+  // std::map key.
+  return a.step_size() < b.step_size();
 }
 bool operator<=(const TimeStepId& a, const TimeStepId& b) { return not(b < a); }
 bool operator>(const TimeStepId& a, const TimeStepId& b) { return b < a; }

src/Time/TimeStepId.hpp

@@ -69,12 +69,19 @@ class TimeStepId {
   double substep_time_{};
 };
 
+/// These comparisons define a total order on all TimeStepIds with the
+/// same direction of time flow.  For any two objects representing
+/// states of the same element, the ordering will match the
+/// computational ordering of those states.  For others, (e.g.,
+/// overlapping steps) it is arbitrary.
+/// @{
 bool operator==(const TimeStepId& a, const TimeStepId& b);
 bool operator!=(const TimeStepId& a, const TimeStepId& b);
 bool operator<(const TimeStepId& a, const TimeStepId& b);
 bool operator<=(const TimeStepId& a, const TimeStepId& b);
 bool operator>(const TimeStepId& a, const TimeStepId& b);
 bool operator>=(const TimeStepId& a, const TimeStepId& b);
+/// @}
 
 std::ostream& operator<<(std::ostream& s, const TimeStepId& id);
 

src/Time/TimeSteppers/AdamsBashforth.cpp

@@ -77,6 +77,10 @@ size_t AdamsBashforth::order() const { return order_; }
 
 size_t AdamsBashforth::error_estimate_order() const { return order_ - 1; }
 
+uint64_t AdamsBashforth::number_of_substeps() const { return 1; }
+
+uint64_t AdamsBashforth::number_of_substeps_for_error() const { return 1; }
+
 size_t AdamsBashforth::number_of_past_steps() const { return order_ - 1; }
 
 double AdamsBashforth::stable_step() const {
@@ -103,6 +107,23 @@ TimeStepId AdamsBashforth::next_time_id(const TimeStepId& current_id,
   return current_id.next_step(time_step);
 }
 
+TimeStepId AdamsBashforth::next_time_id_for_error(
+    const TimeStepId& current_id, const TimeDelta& time_step) const {
+  return next_time_id(current_id, time_step);
+}
+
+bool AdamsBashforth::neighbor_data_required(
+    const TimeStepId& next_substep_id,
+    const TimeStepId& neighbor_data_id) const {
+  return neighbor_data_id < next_substep_id;
+}
+
+bool AdamsBashforth::neighbor_data_required(
+    const double dense_output_time, const TimeStepId& neighbor_data_id) const {
+  return evolution_less<double>{neighbor_data_id.time_runs_forward()}(
+      neighbor_data_id.substep_time(), dense_output_time);
+}
+
 void AdamsBashforth::pup(PUP::er& p) {
   LtsTimeStepper::pup(p);
   p | order_;

src/Time/TimeSteppers/AdamsBashforth.hpp

@@ -3,7 +3,8 @@
 
 #pragma once
 
-#include <vector>
+#include <cstddef>
+#include <cstdint>
 
 #include "Options/String.hpp"
 #include "Time/TimeStepId.hpp"
@@ -216,13 +217,28 @@ class AdamsBashforth : public LtsTimeStepper {
 
   size_t error_estimate_order() const override;
 
+  uint64_t number_of_substeps() const override;
+
+  uint64_t number_of_substeps_for_error() const override;
+
   size_t number_of_past_steps() const override;
 
   double stable_step() const override;
 
   TimeStepId next_time_id(const TimeStepId& current_id,
                           const TimeDelta& time_step) const override;
 
+  TimeStepId next_time_id_for_error(const TimeStepId& current_id,
+                                    const TimeDelta& time_step) const override;
+
+  bool neighbor_data_required(
+      const TimeStepId& next_substep_id,
+      const TimeStepId& neighbor_data_id) const override;
+
+  bool neighbor_data_required(
+      double dense_output_time,
+      const TimeStepId& neighbor_data_id) const override;
+
   WRAPPED_PUPable_decl_template(AdamsBashforth);  // NOLINT
 
   explicit AdamsBashforth(CkMigrateMessage* /*unused*/) {}

src/Time/TimeSteppers/LtsTimeStepper.hpp

@@ -3,7 +3,6 @@
 
 #pragma once
 
-#include <cstdint>
 #include <pup.h>
 #include <type_traits>
 
@@ -81,6 +80,22 @@ class LtsTimeStepper : public virtual TimeStepper {
 #undef LTS_TIME_STEPPER_DECLARE_VIRTUALS_IMPL
   /// \endcond
 
+  /// \brief Check whether a neighbor record is needed for boundary
+  /// output.
+  ///
+  /// In order to perform boundary output, all records from the
+  /// neighbor with `TimeStepId`s for which this method returns true
+  /// should have been added to the history.  Versions are provided
+  /// for a substep and for dense output.
+  /// @{
+  virtual bool neighbor_data_required(
+      const TimeStepId& next_substep_id,
+      const TimeStepId& neighbor_data_id) const = 0;
+
+  virtual bool neighbor_data_required(
+      double dense_output_time, const TimeStepId& neighbor_data_id) const = 0;
+  /// @}
+
   /// \brief Compute the change in a boundary quantity due to the
   /// coupling on the interface.
   ///
@@ -138,17 +153,6 @@ class LtsTimeStepper : public virtual TimeStepper {
                                          history.local(), history.remote(),
                                          history.evaluator(coupling), time);
   }
-
-  /// Substep LTS integrators are not supported, so this is always 1.
-  uint64_t number_of_substeps() const final { return 1; }
-
-  /// Substep LTS integrators are not supported, so this is always 1.
-  uint64_t number_of_substeps_for_error() const final { return 1; }
-
-  TimeStepId next_time_id_for_error(const TimeStepId& current_id,
-                                    const TimeDelta& time_step) const final {
-    return next_time_id(current_id, time_step);
-  }
 };
 
 /// \cond

tests/Unit/Time/TimeSteppers/Test_AdamsBashforth.cpp

@@ -342,10 +342,31 @@ void check_lts_vts() {
     }
   }
 }
+
+void test_neighbor_data_required() {
+  // Test is order-independent
+  const TimeSteppers::AdamsBashforth stepper(4);
+  const Slab slab(0.0, 1.0);
+  CHECK(not stepper.neighbor_data_required(TimeStepId(true, 0, slab.start()),
+                                           TimeStepId(true, 0, slab.start())));
+  CHECK(not stepper.neighbor_data_required(TimeStepId(true, 0, slab.start()),
+                                           TimeStepId(true, 0, slab.end())));
+  CHECK(stepper.neighbor_data_required(TimeStepId(true, 0, slab.end()),
+                                       TimeStepId(true, 0, slab.start())));
+
+  CHECK(not stepper.neighbor_data_required(TimeStepId(false, 0, slab.end()),
+                                           TimeStepId(false, 0, slab.end())));
+  CHECK(not stepper.neighbor_data_required(TimeStepId(false, 0, slab.end()),
+                                           TimeStepId(false, 0, slab.start())));
+  CHECK(stepper.neighbor_data_required(TimeStepId(false, 0, slab.start()),
+                                       TimeStepId(false, 0, slab.end())));
+}
 }  // namespace
 
 SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsBashforth.Boundary",
                   "[Unit][Time]") {
+  test_neighbor_data_required();
+
   // No local stepping
   for (size_t order = 1; order < 9; ++order) {
     INFO(order);