Limit initial time step to avoid control system deadlocks

src/ControlSystem/Actions/InitializeMeasurements.hpp

@@ -3,6 +3,7 @@
 
 #pragma once
 
+#include <algorithm>
 #include <cstddef>
 #include <limits>
 #include <memory>
@@ -23,11 +24,20 @@
 #include "Parallel/AlgorithmExecution.hpp"
 #include "Parallel/GlobalCache.hpp"
 #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
+#include "Time/ChooseLtsStepSize.hpp"
+#include "Time/Slab.hpp"
+#include "Time/Time.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/MakeVector.hpp"
 #include "Utilities/TMPL.hpp"
 
 /// \cond
+namespace Tags {
+struct TimeStep;
+}  // namespace Tags
+namespace domain::Tags {
+struct FunctionsOfTime;
+}  // namespace domain::Tags
 namespace evolution::Tags {
 struct EventsAndDenseTriggers;
 }  // namespace evolution::Tags
@@ -119,6 +129,39 @@ struct InitializeMeasurements {
               });
         },
         make_not_null(&box));
+
+    // Ensure that the initial time step is small enough that we don't
+    // need to perform any measurements to complete it.  This is only
+    // necessary for TimeSteppers that use the self-start procedure,
+    // because they cannot adjust their step size on the first step
+    // after self-start, but it isn't harmful to do it in other cases.
+    //
+    // Unlike in the steady-state step-limiting code, we don't do
+    // anything clever looking at measurement times or planning ahead
+    // for future steps.  Avoiding a single non-ideal step isn't worth
+    // the added complexity.
+    double earliest_expiration = std::numeric_limits<double>::infinity();
+    for (const auto& fot :
+         Parallel::get<domain::Tags::FunctionsOfTime>(cache)) {
+      earliest_expiration =
+          std::min(earliest_expiration, fot.second->time_bounds()[1]);
+    }
+    const auto& time_step = db::get<::Tags::TimeStep>(box);
+    const auto start_time = time_step.slab().start();
+    if ((start_time + time_step).value() > earliest_expiration) {
+      db::mutate<::Tags::TimeStep>(
+          [&](const gsl::not_null<TimeDelta*> step) {
+            if constexpr (Metavariables::local_time_stepping) {
+              *step = choose_lts_step_size(
+                  start_time,
+                  0.99 * (earliest_expiration - start_time.value()));
+            } else {
+              *step = Slab(start_time.value(), earliest_expiration).duration();
+            }
+          },
+          make_not_null(&box));
+    }
+
     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
   }
 };

src/ControlSystem/CMakeLists.txt

@@ -57,6 +57,7 @@ target_link_libraries(
   Logging
   Parallel
   Serialization
+  Time
   Utilities
   INTERFACE
   Observer

src/Time/Actions/SelfStartActions.hpp

@@ -223,16 +223,13 @@ struct Initialize {
             ? 0
             : db::get<::Tags::TimeStepper<TimeStepper>>(box).order() - 1;
 
-    TimeDelta self_start_step = initial_step;
-
-    // Make sure the starting procedure fits in a slab.
-    if (abs(self_start_step * values_needed).fraction() >= 1) {
-      // Decrease the step so that the generated history will be
-      // entirely before the first step.  This ensures we will not
-      // generate any duplicate times in the history as we start the
-      // real evolution.
-      self_start_step /= (values_needed + 1);
-    }
+    // Decrease the step so that the generated history will be
+    // entirely before the first step.  This ensures we will not
+    // generate any duplicate times in the history as we start the
+    // real evolution and that the starting procedure does not require
+    // any more information (such as function-of-time values) than the
+    // first real step.
+    const TimeDelta self_start_step = initial_step / (values_needed + 1);
 
     StoreInitialValues<
         tmpl::push_back<detail::vars_to_save<System>, ::Tags::TimeStep,

src/Time/TimeSteppers/AdamsBashforth.cpp

@@ -173,18 +173,20 @@ void AdamsBashforth::update_u_common(const gsl::not_null<T*> u,
 template <typename T>
 bool AdamsBashforth::can_change_step_size_impl(
     const TimeStepId& time_id, const ConstUntypedHistory<T>& history) const {
-  // We need to forbid local time-stepping before initialization is
-  // complete.  The self-start procedure itself should never consider
-  // changing the step size, but we need to wait during the main
-  // evolution until the self-start history has been replaced with
-  // "real" values.
-  const evolution_less<Time> less{time_id.time_runs_forward()};
+  // We need to prevent the next step from occurring at the same time
+  // as one already in the history.  The self-start code ensures this
+  // can't happen during self-start, and it clearly can't happen
+  // during normal evolution where the steps are monotonic, but during
+  // the transition between them we have to worry about a step being
+  // placed on a self-start time.  The self-start algorithm guarantees
+  // the final state is safe for constant-time-step evolution, so we
+  // just force that until we've passed all the self-start times.
+  const evolution_less_equal<Time> less_equal{time_id.time_runs_forward()};
   return not ::SelfStart::is_self_starting(time_id) and
-         (history.size() == 0 or
-          (less(history.back().time_step_id.step_time(),
-                time_id.step_time()) and
-           std::is_sorted(history_time_iterator(history.begin()),
-                          history_time_iterator(history.end()), less)));
+         alg::all_of(history, [&](const auto& record) {
+           return less_equal(record.time_step_id.step_time(),
+                             time_id.step_time());
+         });
 }
 
 template <typename T>

src/Time/TimeSteppers/AdamsMoultonPc.cpp

@@ -199,10 +199,14 @@ bool AdamsMoultonPc::dense_update_u_impl(const gsl::not_null<T*> u,
 template <typename T>
 bool AdamsMoultonPc::can_change_step_size_impl(
     const TimeStepId& time_id, const ConstUntypedHistory<T>& history) const {
-  // We need to forbid local time-stepping before initialization is
-  // complete.  We need to wait during the main evolution so we don't
-  // increase the step size and step to a time still in use from
-  // self-start, as that would cause the coefficients to diverge.
+  // We need to prevent the next step from occurring at the same time
+  // as one already in the history.  The self-start code ensures this
+  // can't happen during self-start, and it clearly can't happen
+  // during normal evolution where the steps are monotonic, but during
+  // the transition between them we have to worry about a step being
+  // placed on a self-start time.  The self-start algorithm guarantees
+  // the final state is safe for constant-time-step evolution, so we
+  // just force that until we've passed all the self-start times.
   const evolution_less_equal<Time> less_equal{time_id.time_runs_forward()};
   return not ::SelfStart::is_self_starting(time_id) and
          alg::all_of(history, [&](const auto& record) {

tests/Unit/ControlSystem/Actions/Test_InitializeMeasurements.cpp

@@ -41,8 +41,11 @@
 #include "Parallel/Phase.hpp"
 #include "Parallel/PhaseDependentActionList.hpp"
 #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
+#include "Time/Slab.hpp"
 #include "Time/Tags/Time.hpp"
+#include "Time/Tags/TimeStep.hpp"
 #include "Time/Tags/TimeStepId.hpp"
+#include "Utilities/CartesianProduct.hpp"
 #include "Utilities/GetOutput.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/ProtocolHelpers.hpp"
@@ -138,7 +141,7 @@ struct Component {
   using simple_tags_from_options =
       tmpl::list<evolution::Tags::EventsAndDenseTriggers>;
 
-  using simple_tags = tmpl::list<Tags::TimeStepId, Tags::Time>;
+  using simple_tags = tmpl::list<Tags::TimeStepId, Tags::Time, Tags::TimeStep>;
   using compute_tags = tmpl::list<Parallel::Tags::FromGlobalCache<
       ::domain::Tags::FunctionsOfTimeInitialize>>;
   using const_global_cache_tags = tmpl::list<control_system::Tags::Verbosity>;
@@ -153,8 +156,10 @@ struct Component {
           control_system::Actions::InitializeMeasurements<control_systems>>>>;
 };
 
+template <bool LocalTimeStepping>
 struct Metavariables {
   using component_list = tmpl::list<Component<Metavariables>>;
+  static constexpr bool local_time_stepping = LocalTimeStepping;
 
   struct factory_creation
       : tt::ConformsTo<Options::protocols::FactoryCreation> {
@@ -166,24 +171,36 @@ struct Metavariables {
   };
 };
 
+template <bool LocalTimeStepping>
 void test_initialize_measurements(const bool ab_active, const bool c_active) {
+  CAPTURE(LocalTimeStepping);
   CAPTURE(ab_active);
   CAPTURE(c_active);
 
-  register_factory_classes_with_charm<Metavariables>();
+  register_factory_classes_with_charm<Metavariables<LocalTimeStepping>>();
   register_classes_with_charm<
       domain::FunctionsOfTime::PiecewisePolynomial<0>>();
 
-  using MockRuntimeSystem = ActionTesting::MockRuntimeSystem<Metavariables>;
-  using component = Component<Metavariables>;
+  using MockRuntimeSystem =
+      ActionTesting::MockRuntimeSystem<Metavariables<LocalTimeStepping>>;
+  using component = Component<Metavariables<LocalTimeStepping>>;
   const component* const component_p = nullptr;
 
   // Details shouldn't matter
   const double initial_time = 2.0;
   const size_t measurements_per_update = 6;
+  const double initial_timescale = 1.5;
+  const double fot_expiration = 4.0;
+
+  const double step_to_expiration_ratio = c_active ? 3.0 : 0.2;
+
+  const auto initial_slab = Slab::with_duration_from_start(
+      initial_time, step_to_expiration_ratio * (fot_expiration - initial_time));
+  const auto initial_time_step = initial_slab.duration();
+
   const auto timescale =
       std::make_unique<domain::FunctionsOfTime::PiecewisePolynomial<0>>(
-          0.0, std::array{DataVector{1.0}}, 2.0);
+          0.0, std::array{DataVector{initial_timescale}}, 2.0);
   const auto inactive_timescale =
       std::make_unique<domain::FunctionsOfTime::PiecewisePolynomial<0>>(
           0.0, std::array{DataVector{std::numeric_limits<double>::infinity()}},
@@ -195,19 +212,24 @@ void test_initialize_measurements(const bool ab_active, const bool c_active) {
                      (ab_active ? timescale : inactive_timescale)->get_clone());
   timescales.emplace("C",
                      (c_active ? timescale : inactive_timescale)->get_clone());
+
+  const auto fot =
+      std::make_unique<domain::FunctionsOfTime::PiecewisePolynomial<0>>(
+          0.0, std::array{DataVector{1.0}}, fot_expiration);
+
   std::unordered_map<std::string,
                      std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
       functions;
-  functions.emplace("A", timescale->get_clone());
-  functions.emplace("B", timescale->get_clone());
-  functions.emplace("C", timescale->get_clone());
+  functions.emplace("A", fot->get_clone());
+  functions.emplace("B", fot->get_clone());
+  functions.emplace("C", fot->get_clone());
 
   MockRuntimeSystem runner{{::Verbosity::Silent, measurements_per_update},
                            {std::move(functions), std::move(timescales)}};
   ActionTesting::emplace_array_component_and_initialize<component>(
       make_not_null(&runner), ActionTesting::NodeId{0},
       ActionTesting::LocalCoreId{0}, 0,
-      {Tags::TimeStepId::type{true, 0, {}}, initial_time},
+      {Tags::TimeStepId::type{true, 0, {}}, initial_time, initial_time_step},
       evolution::Tags::EventsAndDenseTriggers::type{});
 
   // InitializeRunEventsAndDenseTriggers
@@ -240,6 +262,16 @@ void test_initialize_measurements(const bool ab_active, const bool c_active) {
         (c_active ? std::nullopt
                   : std::optional(std::numeric_limits<double>::infinity())));
 
+  if (c_active) {
+    // 4 is because of the 3.0 in the step_to_expiration_ratio above.
+    const auto reduced_time_step =
+        LocalTimeStepping ? initial_time_step / 4
+                          : Slab(initial_time, fot_expiration).duration();
+    CHECK(db::get<::Tags::TimeStep>(box) == reduced_time_step);
+  } else {
+    CHECK(db::get<::Tags::TimeStep>(box) == initial_time_step);
+  }
+
   auto& events_and_dense_triggers =
       db::get_mutable_reference<evolution::Tags::EventsAndDenseTriggers>(
           make_not_null(&box));
@@ -266,9 +298,10 @@ void test_initialize_measurements(const bool ab_active, const bool c_active) {
 
 SPECTRE_TEST_CASE("Unit.ControlSystem.InitializeMeasurements",
                   "[ControlSystem][Unit]") {
-  test_initialize_measurements(true, true);
-  test_initialize_measurements(true, false);
-  test_initialize_measurements(false, true);
-  test_initialize_measurements(false, false);
+  for (const auto& [ab_active, c_active] :
+       cartesian_product(std::array{true, false}, std::array{true, false})) {
+    test_initialize_measurements<false>(ab_active, c_active);
+    test_initialize_measurements<true>(ab_active, c_active);
+  }
 }
 }  // namespace

tests/Unit/ControlSystem/CMakeLists.txt

@@ -44,5 +44,6 @@ target_link_libraries(
   Parallel
   PostNewtonianHelpers
   SphericalHarmonics
+  Time
   Utilities
 )

tests/Unit/Time/TimeSteppers/Test_AdamsBashforth.cpp

@@ -77,7 +77,7 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsBashforth", "[Unit][Time]") {
   };
   CHECK(can_change(start, mid, end));
   CHECK_FALSE(can_change(start, end, mid));
-  CHECK_FALSE(can_change(mid, start, end));
+  CHECK(can_change(mid, start, end));
   CHECK_FALSE(can_change(mid, end, start));
   CHECK_FALSE(can_change(end, start, mid));
   CHECK_FALSE(can_change(end, mid, start));
@@ -148,7 +148,7 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsBashforth.Backwards",
   CHECK_FALSE(can_change(start, mid, end));
   CHECK_FALSE(can_change(start, end, mid));
   CHECK_FALSE(can_change(mid, start, end));
-  CHECK_FALSE(can_change(mid, end, start));
+  CHECK(can_change(mid, end, start));
   CHECK_FALSE(can_change(end, start, mid));
   CHECK(can_change(end, mid, start));
 }