Avoid assumptions about history in dense output

The generic dense output code previously made assumptions about how
the time stepper used substeps when initializing the variables for
dense output.  I have accidentally broken this multiple times while
experimenting with various things, and the problems are sometimes
subtle.

src/Evolution/Actions/RunEventsAndDenseTriggers.hpp

@@ -253,7 +253,7 @@ struct RunEventsAndDenseTriggers {
                   gsl::not_null<typename variables_tag::type*> vars,
                   const TimeStepper& stepper,
                   const typename history_tag::type& history) {
-                *vars = *history.complete_step_start().value;
+                *vars = *history.step_start(next_trigger).value;
                 dense_output_succeeded =
                     stepper.dense_update_u(vars, history, next_trigger);
               },

src/Time/History.hpp

@@ -14,6 +14,7 @@
 #include "DataStructures/DataBox/Prefixes.hpp"
 #include "DataStructures/MathWrapper.hpp"
 #include "DataStructures/StaticDeque.hpp"
+#include "Time/EvolutionOrdering.hpp"
 #include "Time/TimeStepId.hpp"
 #include "Utilities/ContainsAllocations.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
@@ -622,10 +623,10 @@ class History
   /// the time steppers can manage their history.
   const Vars& latest_value() const;
 
-  /// Get the record for the start of the latest complete step.  This
-  /// is the step with substeps stored if there is one, and otherwise
-  /// the most recent step.
-  const StepRecord<Vars>& complete_step_start() const;
+  /// Get the record for the start of the step containing the passed
+  /// time.  If the time matches a step time exactly, the record for
+  /// that step is returned.
+  const StepRecord<Vars>& step_start(double time) const;
 
   /// Check whether we are at the start of a step, i.e, the most
   /// recent entry in the history is not a substep.
@@ -792,16 +793,22 @@ const Vars& History<Vars>::latest_value() const {
 }
 
 template <typename Vars>
-const StepRecord<Vars>& History<Vars>::complete_step_start() const {
-  ASSERT(not this->empty(), "History is empty");
-  if (substep_values_.empty()) {
-    return this->back();
-  } else {
-    const TimeStepId substep_id = substep_values_.front().time_step_id;
-    const TimeStepId step_id(substep_id.time_runs_forward(),
-                             substep_id.slab_number(), substep_id.step_time());
-    return (*this)[step_id];
+const StepRecord<Vars>& History<Vars>::step_start(const double time) const {
+  // Search starting at the end to handle self-start correctly (and
+  // because the result under usual use is one of the last two
+  // entries).
+  const auto first_step = this->begin();
+  auto step = this->end();
+  ASSERT(step != first_step, "History is empty");
+  --step;
+  const evolution_less<double> before{step->time_step_id.time_runs_forward()};
+  while (before(time, step->time_step_id.step_time().value())) {
+    ASSERT(step != first_step,
+           "Start of step at time " << time << " is before start of history "
+           << first_step->time_step_id);
+    --step;
   }
+  return *step;
 }
 
 template <typename Vars>

src/Time/TimeSteppers/Rk3HesthavenSsp.cpp

@@ -122,13 +122,16 @@ template <typename T>
 bool Rk3HesthavenSsp::dense_update_u_impl(const gsl::not_null<T*> u,
                                           const ConstUntypedHistory<T>& history,
                                           const double time) const {
+  const double step_start = history.front().time_step_id.step_time().value();
+  if (time == step_start) {
+    return true;
+  }
   if (not history.at_step_start()) {
     return false;
   }
-  const double step_start = history.front().time_step_id.step_time().value();
   const double step_end = history.back().time_step_id.step_time().value();
   const evolution_less<double> before{step_end > step_start};
-  if (history.size() == 1 or before(step_end, time)) {
+  if (history.size() == 1 or not before(time, step_end)) {
     return false;
   }
   const double step_size = step_end - step_start;

src/Time/TimeSteppers/RungeKutta.cpp

@@ -178,13 +178,16 @@ template <typename T>
 bool RungeKutta::dense_update_u_impl(const gsl::not_null<T*> u,
                                      const ConstUntypedHistory<T>& history,
                                      const double time) const {
+  const double step_start = history.front().time_step_id.step_time().value();
+  if (time == step_start) {
+    return true;
+  }
   if (not history.at_step_start()) {
     return false;
   }
-  const double step_start = history.front().time_step_id.step_time().value();
   const double step_end = history.back().time_step_id.step_time().value();
   const evolution_less<double> before{step_end > step_start};
-  if (history.size() == 1 or before(step_end, time)) {
+  if (history.size() == 1 or not before(time, step_end)) {
     return false;
   }
   const double step_size = step_end - step_start;

src/Time/TimeSteppers/TimeStepper.hpp

@@ -161,7 +161,7 @@ class TimeStepper : public PUP::able {
   ///
   /// The change from the partial step will be added to the initial
   /// value, so \p u should generally be initialized to
-  /// `*history.complete_step_start().value`.  (TimeStepper
+  /// `*history.step_start(time).value`.  (TimeStepper
   /// implementations are required to keep this value in the history.)
   ///
   /// Derived classes must implement this as a function with signature

tests/Unit/Helpers/Time/TimeSteppers/LtsHelpers.cpp

@@ -221,7 +221,7 @@ class Element {
 
   double dense_output(const double time) {
     REQUIRE(process_messages(time));
-    double dense_result = *volume_history_.complete_step_start().value;
+    double dense_result = *volume_history_.step_start(time).value;
     // Can skip the volume update because all derivatives are zero.
     stepper_->boundary_dense_output(make_not_null(&dense_result),
                                     boundary_history_, time, coupling_);

tests/Unit/Helpers/Time/TimeSteppers/TimeStepperTestUtils.cpp

@@ -381,9 +381,9 @@ void check_dense_output(
         history.insert(time_id, y, y);
         if (before(time, (time_id.step_time() + step).value())) {
           // Make sure the initial value is preserved.
-          y = 2.0 * *history.complete_step_start().value;
+          y = 2.0 * *history.step_start(time).value;
           if (stepper.dense_update_u(make_not_null(&y), history, time)) {
-            return y - *history.complete_step_start().value;
+            return y - *history.step_start(time).value;
           }
           REQUIRE(not before(time, time_id.step_time().value()));
           CHECK(not stepper.monotonic());

tests/Unit/Time/Test_History.cpp

@@ -117,7 +117,13 @@ void test_history() {
                          make_deriv(-10.0));
   // [(0, -1, -10)] []
   CHECK(const_history.size() == 1);
-  CHECK(&const_history.complete_step_start() == &const_history[0]);
+  CHECK(&const_history.step_start(0.0) == &const_history[0]);
+  CHECK(&const_history.step_start(0.1) == &const_history[0]);
+#ifdef SPECTRE_DEBUG
+  CHECK_THROWS_WITH(
+      const_history.step_start(-0.1),
+      Catch::Matchers::ContainsSubstring("before start of history"));
+#endif  // SPECTRE_DEBUG
   history.insert_initial(
       TimeStepId(true, -1, slab.start() - Slab(-1.0, 0.0).duration() / 4),
       History::no_value, make_deriv(-20.0));
@@ -142,7 +148,14 @@ void test_history() {
   CHECK(const_history[2].derivative == make_deriv(-10.0));
   CHECK(const_history.latest_value() == *const_history.back().value);
   CHECK(&const_history.latest_value() == &*const_history.back().value);
-  CHECK(&const_history.complete_step_start() == &const_history[2]);
+  CHECK(&const_history.step_start(0.0) == &const_history[2]);
+  CHECK(&const_history.step_start(0.1) == &const_history[2]);
+  CHECK(&const_history.step_start(-0.1) == &const_history[1]);
+#ifdef SPECTRE_DEBUG
+  CHECK_THROWS_WITH(
+      const_history.step_start(-1.0),
+      Catch::Matchers::ContainsSubstring("before start of history"));
+#endif  // SPECTRE_DEBUG
 
   history[0].derivative = make_deriv(-300.0);
   // [(-1/2, -3, -300), (-1/4, X, -20), (0, -1, -10)] []
@@ -297,7 +310,10 @@ void test_history() {
   CHECK(as_const(const_history.substeps()).size() == 1);
   CHECK(as_const(const_history.substeps())[0].derivative == make_deriv(40.0));
   CHECK(as_const(history.substeps()).size() == 1);
-  CHECK(&const_history.complete_step_start() == &const_history[1]);
+  CHECK(&const_history.step_start(0.5) == &const_history[1]);
+  CHECK(&const_history.step_start(0.7) == &const_history[1]);
+  CHECK(&const_history.step_start(0.3) == &const_history[0]);
+
   as_const(history.substeps())[0].derivative = make_deriv(400.0);
   // [(1/4, X, 10), (1/2, 2, 20)] [1/2: (1, 4, 400)]
   CHECK(as_const(const_history.substeps())[0].derivative == make_deriv(400.0));
@@ -395,7 +411,10 @@ void test_history() {
   CHECK(const_history.at_step_start());
   CHECK(static_cast<const ConstUntyped&>(const_history.untyped())
             .at_step_start());
-  CHECK(&const_history.complete_step_start() == &const_history[1]);
+  CHECK(&const_history.step_start(0.5) == &const_history[1]);
+  CHECK(&const_history.step_start(0.7) == &const_history[1]);
+  CHECK(&const_history.step_start(1.1) == &const_history[2]);
+  CHECK(&const_history.step_start(1.0) == &const_history[2]);
 
   // Test transform when the substeps are not associated with the last
   // step.  This causes errors in a naive implementation.