Merge pull request #4748 from wthrowe/fix_dense_output

Have time steppers overwrite in update_u again

src/Evolution/Actions/RunEventsAndDenseTriggers.hpp

@@ -207,13 +207,8 @@ struct RunEventsAndDenseTriggers {
         return {Parallel::AlgorithmExecution::Continue, std::nullopt};
       }
 
-      // This can only be true the first time through the loop,
-      // because triggers are not allowed to reschedule for the time
-      // they just triggered at.  This check is primarily to avoid
-      // special-case bookkeeping for the initial simulation time.
-      const bool already_at_correct_time =
-          db::get<::Tags::Time>(box) == next_trigger;
-      if (not already_at_correct_time) {
+      // Avoid invalidating compute items unless necessary.
+      if (db::get<::Tags::Time>(box) != next_trigger) {
         time_restorer.save();
         db::mutate<::Tags::Time>(
             make_not_null(&box),
@@ -228,53 +223,46 @@ struct RunEventsAndDenseTriggers {
       switch (triggered) {
         case TriggeringState::NotReady:
           return {Parallel::AlgorithmExecution::Retry, std::nullopt};
-        case TriggeringState::NeedsEvolvedVariables:
-          if (not already_at_correct_time) {
-            bool ready = true;
-            tmpl::for_each<Postprocessors>([&](auto postprocessor_v) {
-              using postprocessor = tmpl::type_from<decltype(postprocessor_v)>;
-              if (ready) {
-                if (not postprocessor::is_ready(make_not_null(&box),
-                                                make_not_null(&inboxes), cache,
-                                                array_index, component)) {
-                  ready = false;
-                }
+        case TriggeringState::NeedsEvolvedVariables: {
+          bool ready = true;
+          tmpl::for_each<Postprocessors>([&](auto postprocessor_v) {
+            using postprocessor = tmpl::type_from<decltype(postprocessor_v)>;
+            if (ready) {
+              if (not postprocessor::is_ready(make_not_null(&box),
+                                              make_not_null(&inboxes), cache,
+                                              array_index, component)) {
+                ready = false;
               }
-            });
-            if (not ready) {
-              return {Parallel::AlgorithmExecution::Retry, std::nullopt};
-            }
-
-            using history_tag = ::Tags::HistoryEvolvedVariables<variables_tag>;
-            bool dense_output_succeeded = false;
-            // Restore the evolved variables if they have been
-            // previously altered so the time stepper starts in the
-            // correct state, then save them if they have not been
-            // previously altered.  (So one of these two lines is a
-            // no-op.)
-            variables_restorer.restore();
-            variables_restorer.save();
-            db::mutate<variables_tag>(
-                make_not_null(&box),
-                [&dense_output_succeeded, &next_trigger](
-                    gsl::not_null<typename variables_tag::type*> vars,
-                    const TimeStepper& stepper,
-                    const typename history_tag::type& history) {
-                  dense_output_succeeded =
-                      stepper.dense_update_u(vars, history, next_trigger);
-                },
-                db::get<::Tags::TimeStepper<>>(box), db::get<history_tag>(box));
-            if (not dense_output_succeeded) {
-              // Need to take another time step
-              return {Parallel::AlgorithmExecution::Continue, std::nullopt};
             }
+          });
+          if (not ready) {
+            return {Parallel::AlgorithmExecution::Retry, std::nullopt};
+          }
 
-            postprocessor_restorer.save();
-            tmpl::for_each<Postprocessors>([&box](auto postprocessor_v) {
-              using postprocessor = tmpl::type_from<decltype(postprocessor_v)>;
-              db::mutate_apply<postprocessor>(make_not_null(&box));
-            });
+          using history_tag = ::Tags::HistoryEvolvedVariables<variables_tag>;
+          bool dense_output_succeeded = false;
+          variables_restorer.save();
+          db::mutate<variables_tag>(
+              make_not_null(&box),
+              [&dense_output_succeeded, &next_trigger](
+                  gsl::not_null<typename variables_tag::type*> vars,
+                  const TimeStepper& stepper,
+                  const typename history_tag::type& history) {
+                dense_output_succeeded =
+                    stepper.dense_update_u(vars, history, next_trigger);
+              },
+              db::get<::Tags::TimeStepper<>>(box), db::get<history_tag>(box));
+          if (not dense_output_succeeded) {
+            // Need to take another time step
+            return {Parallel::AlgorithmExecution::Continue, std::nullopt};
           }
+
+          postprocessor_restorer.save();
+          tmpl::for_each<Postprocessors>([&box](auto postprocessor_v) {
+            using postprocessor = tmpl::type_from<decltype(postprocessor_v)>;
+            db::mutate_apply<postprocessor>(make_not_null(&box));
+          });
+        }
           [[fallthrough]];
         default:
           break;

src/Time/TimeSteppers/AdamsBashforth.cpp

@@ -198,7 +198,9 @@ void clean_history(const MutableUntypedHistory<T>& history) {
   while (history.size() > history.integration_order()) {
     history.pop_front();
   }
-  history.discard_value(history.back().time_step_id);
+  if (history.size() > 1) {
+    history.discard_value(history[history.size() - 2].time_step_id);
+  }
 }
 }  // namespace
 
@@ -257,7 +259,6 @@ bool AdamsBashforth::update_u_impl(
     const gsl::not_null<T*> u, const gsl::not_null<T*> u_error,
     const MutableUntypedHistory<T>& history, const TimeDelta& time_step) const {
   clean_history(history);
-  *u_error = *u;
   update_u_common(u, history, time_step, history.integration_order());
   // the error estimate is only useful once the history has enough elements to
   // do more than one order of step
@@ -294,6 +295,7 @@ void AdamsBashforth::update_u_common(const gsl::not_null<T*> u,
       get_coefficients(history_time_iterator(history_start),
                        history_time_iterator(history.end()), time_step);
 
+  *u = *history.back().value;
   auto coefficient = coefficients.rbegin();
   for (auto history_entry = history_start;
        history_entry != history.end();
@@ -364,6 +366,11 @@ void AdamsBashforth::boundary_dense_output_impl(
     const gsl::not_null<T*> result,
     const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
     const double time) const {
+  if ((coupling.local_end() - 1)->value() == time) {
+    // Nothing to do.  The requested time is the start of the step,
+    // which is the input value of `result`.
+    return;
+  }
   return boundary_impl(result, coupling, ApproximateTime{time});
 }
 

src/Time/TimeSteppers/RungeKutta.cpp

@@ -110,9 +110,11 @@ void update_u_impl_with_tableau(const gsl::not_null<T*> u,
     if (history.size() > 1) {
       history.pop_front();
     }
-    history.discard_value(history.back().time_step_id);
+  } else if (history.substeps().size() == 1) {
+    history.discard_value(history.front().time_step_id);
   } else {
-    history.discard_value(history.substeps().back().time_step_id);
+    history.discard_value(
+        history.substeps()[history.substeps().size() - 2].time_step_id);
   }
   ASSERT(history.size() == 1, "Have more than one step after cleanup.");
 
@@ -125,12 +127,15 @@ void update_u_impl_with_tableau(const gsl::not_null<T*> u,
   if (substep == 0) {
     ASSERT(tableau.substep_coefficients[0].size() == 1,
            "First substep should use one derivative.");
-    *u += tableau.substep_coefficients[0][0] * dt * history.front().derivative;
+    *u = *history.front().value +
+         tableau.substep_coefficients[0][0] * dt * history.front().derivative;
   } else if (substep == number_of_substeps - 1) {
+    *u = *history.substeps().back().value;
     update_between_substeps(u, history, dt,
                             tableau.substep_coefficients[substep - 1],
                             tableau.result_coefficients);
   } else if (substep < number_of_substeps - 1) {
+    *u = *history.substeps().back().value;
     update_between_substeps(u, history, dt,
                             tableau.substep_coefficients[substep - 1],
                             tableau.substep_coefficients[substep]);
@@ -192,6 +197,7 @@ bool RungeKutta::dense_update_u_impl(const gsl::not_null<T*> u,
   if (time == step_end) {
     // Special case necessary for dense output at the initial time,
     // before taking a step.
+    *u = *history.back().value;
     return true;
   }
   const evolution_less<double> before{step_end > step_start};
@@ -206,6 +212,7 @@ bool RungeKutta::dense_update_u_impl(const gsl::not_null<T*> u,
 
   const auto& tableau = butcher_tableau();
 
+  *u = *history.back().value;
   // The Butcher dense output coefficients are given in terms of the
   // start of the step, but we are passed the value at the end of the
   // step, so we have to undo the step first.

src/Time/TimeSteppers/TimeStepper.hpp

@@ -72,7 +72,7 @@ class TimeStepper : public PUP::able {
 #undef TIME_STEPPER_DECLARE_VIRTUALS_IMPL
   /// \endcond
 
-  /// Add the change for the current substep to u.
+  /// Set \p u to the value at the end of the current substep.
   ///
   /// Derived classes must implement this as a function with signature
   ///
@@ -90,8 +90,8 @@ class TimeStepper : public PUP::able {
                             time_step);
   }
 
-  /// Add the change for the current substep to u; report the error measure when
-  /// available.
+  /// Set \p u to the value at the end of the current substep; report the error
+  /// measure when available.
   ///
   /// For a substep method, the error measure will only be available on full
   /// steps. For a multistep method, the error measure will only be available

tests/Unit/Evolution/Actions/Test_RunEventsAndDenseTriggers.cpp

@@ -110,11 +110,13 @@ class TestTrigger : public DenseTrigger {
       Parallel::GlobalCache<Metavariables>& /*cache*/,
       const ArrayIndex& /*array_index*/, const Component* /*component*/,
       const double time) const {
-    if (time == init_time_) {
+    if (time == trigger_time_) {
+      return is_triggered_;
+    } else {
+      // First initialization call.
+      CHECK(time == init_time_);
       return false;
     }
-    CHECK(time == trigger_time_);
-    return is_triggered_;
   }
 
   using next_check_time_argument_tags = tmpl::list<Tags::Time>;
@@ -123,11 +125,13 @@ class TestTrigger : public DenseTrigger {
       Parallel::GlobalCache<Metavariables>& /*cache*/,
       const ArrayIndex& /*array_index*/, const Component* /*component*/,
       const double time) const {
-    if (time == init_time_) {
+    if (time == trigger_time_) {
+      return next_trigger_;
+    } else {
+      // First initialization call.
+      CHECK(time == init_time_);
       return trigger_time_;
     }
-    CHECK(time == trigger_time_);
-    return next_trigger_;
   }
 
   // NOLINTNEXTLINE(google-runtime-references)
@@ -323,11 +327,12 @@ void test(const bool time_runs_forward) {
   const double done_time = (time_runs_forward ? 1.0 : -1.0) *
                            std::numeric_limits<double>::infinity();
   const VarsType initial_vars{1, 8.0};
+  const VarsType stored_vars{1, 123.0};
   const DtVarsType deriv_vars{1, 1.0};
   const VarsType center_vars = initial_vars + 0.5 * step_size * deriv_vars;
 
   const auto set_up_component =
-      [&deriv_vars, &exact_step_size, &initial_vars, &start_time,
+      [&deriv_vars, &exact_step_size, &initial_vars, &start_time, &stored_vars,
        &time_step_id](
           const gsl::not_null<MockRuntimeSystem*> runner,
           const std::vector<std::tuple<double, std::optional<bool>,
@@ -352,7 +357,7 @@ void test(const bool time_runs_forward) {
           ActionTesting::emplace_array_component<component>(
               runner, ActionTesting::NodeId{0}, ActionTesting::LocalCoreId{0},
               0, time_step_id, exact_step_size, start_time,
-              std::optional<double>{}, initial_vars, std::move(history),
+              std::optional<double>{}, stored_vars, std::move(history),
               evolution::EventsAndDenseTriggers(
                   std::move(events_and_dense_triggers)),
               typename evolution::dg::Tags::NeighborMesh<1>::type{},
@@ -427,11 +432,24 @@ void test(const bool time_runs_forward) {
         reschedule ? std::optional{done_time} : std::nullopt;
     set_up_component(&runner, {{step_center, true, next_check, false}});
     CHECK(run_if_ready(make_not_null(&runner)) == reschedule);
-    TestEvent::check_calls({{step_center, initial_vars}});
+    TestEvent::check_calls({{step_center, stored_vars}});
   };
   check_not_needed(true);
   check_not_needed(false);
 
+  // Variables not needed at initial time
+  const auto check_not_needed_initial = [&](const bool reschedule) {
+    MockRuntimeSystem runner{
+        {std::make_unique<TimeSteppers::AdamsBashforth>(1)}};
+    const auto next_check =
+        reschedule ? std::optional{done_time} : std::nullopt;
+    set_up_component(&runner, {{start_time, true, next_check, false}});
+    CHECK(run_if_ready(make_not_null(&runner)) == reschedule);
+    TestEvent::check_calls({{start_time, stored_vars}});
+  };
+  check_not_needed_initial(true);
+  check_not_needed_initial(false);
+
   // Variables needed
   const auto check_needed = [&](const bool reschedule) {
     MockRuntimeSystem runner{
@@ -444,6 +462,18 @@ void test(const bool time_runs_forward) {
   check_needed(true);
   check_needed(false);
 
+  // Variables needed at initial time
+  const auto check_needed_initial = [&](const bool reschedule) {
+    MockRuntimeSystem runner{
+        {std::make_unique<TimeSteppers::AdamsBashforth>(1)}};
+    const auto next_check =
+        reschedule ? std::optional{done_time} : std::nullopt;
+    set_up_component(&runner, {{start_time, true, next_check, true}});
+    TestCase::check_dense(&runner, reschedule, {{start_time, initial_vars}});
+  };
+  check_needed_initial(true);
+  check_needed_initial(false);
+
   // Missing dense output data
   const auto check_missing_dense_data = [&](const bool data_needed) {
     MockRuntimeSystem runner{
@@ -473,7 +503,7 @@ void test(const bool time_runs_forward) {
     } else {
       // If we don't need the data, it shouldn't matter whether it is missing.
       CHECK(run_if_ready(make_not_null(&runner)));
-      TestEvent::check_calls({{step_center, initial_vars}});
+      TestEvent::check_calls({{step_center, stored_vars}});
     }
   };
   check_missing_dense_data(true);

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

@@ -42,6 +42,7 @@ void take_step(
        ++substep) {
     CHECK(time_id.substep() == substep);
     history->insert(time_id, *y, rhs(*y, time_id.substep_time().value()));
+    *y = std::numeric_limits<double>::signaling_NaN();
     stepper.update_u(y, history, step_size);
     time_id = stepper.next_time_id(time_id, step_size);
   }
@@ -60,6 +61,7 @@ void take_step_and_check_error(
        ++substep) {
     CHECK(time_id.substep() == substep);
     history->insert(time_id, *y, rhs(*y, time_id.substep_time().value()));
+    *y = std::numeric_limits<double>::signaling_NaN();
     bool error_updated = stepper.update_u(y, y_error, history, step_size);
     CAPTURE(substep);
     REQUIRE((substep == stepper.number_of_substeps_for_error() - 1) ==
@@ -454,6 +456,7 @@ void check_dense_output(const TimeStepper& stepper,
       auto step = step_size;
       for (;;) {
         history.insert(time_id, y, y);
+        y = std::numeric_limits<double>::signaling_NaN();
         if (not before((time_id.step_time() + step).value(), time)) {
           if (stepper.dense_update_u(make_not_null(&y), history, time)) {
             return y;