Add option for monotonic AdamsMoulton dense output

src/ControlSystem/Actions/LimitTimeStep.hpp

@@ -96,10 +96,9 @@ struct LimitTimeStep {
     }
 
     const auto& time_stepper = db::get<::Tags::TimeStepper<TimeStepper>>(box);
-    if (time_stepper.number_of_substeps() == 1) {
-      // If there are no substeps, there is no reason to limit the
-      // step size so substeps can be evaluated.  Single-substep FSAL
-      // methods could be problematic, but we don't have any of those.
+    if (time_stepper.monotonic()) {
+      // Monotonic steppers order operations in the same manner at the
+      // control system, so they cannot introduce deadlocks.
       return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     }
 

src/Executables/TimeStepperSummary/TimeStepperSummary.cpp

@@ -26,6 +26,7 @@ class ImexTimeStepper;
 class LtsTimeStepper;
 namespace TimeSteppers {
 class AdamsBashforth;
+template <bool Monotonic>
 class AdamsMoultonPc;
 }  // namespace TimeSteppers
 
@@ -35,7 +36,9 @@ extern "C" void CkRegisterMainModule(void) {}
 
 namespace {
 using time_steppers_taking_order =
-    tmpl::list<TimeSteppers::AdamsBashforth, TimeSteppers::AdamsMoultonPc>;
+    tmpl::list<TimeSteppers::AdamsBashforth,
+               TimeSteppers::AdamsMoultonPc<false>,
+               TimeSteppers::AdamsMoultonPc<true>>;
 
 const char* const program_help =
     "Print various properties about SpECTRE's time steppers.  Abbreviations\n"

src/Time/TimeSteppers/AdamsBashforth.cpp

@@ -101,6 +101,8 @@ double AdamsBashforth::stable_step() const {
   return 1. / invstep;
 }
 
+bool AdamsBashforth::monotonic() const { return true; }
+
 TimeStepId AdamsBashforth::next_time_id(const TimeStepId& current_id,
                                         const TimeDelta& time_step) const {
   ASSERT(current_id.substep() == 0, "Adams-Bashforth should not have substeps");

src/Time/TimeSteppers/AdamsBashforth.hpp

@@ -225,6 +225,8 @@ class AdamsBashforth : public LtsTimeStepper {
 
   double stable_step() const override;
 
+  bool monotonic() const override;
+
   TimeStepId next_time_id(const TimeStepId& current_id,
                           const TimeDelta& time_step) const override;
 

src/Time/TimeSteppers/AdamsMoultonPc.cpp

@@ -16,18 +16,19 @@
 #include "Utilities/Algorithm.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/Gsl.hpp"
 
 namespace TimeSteppers {
 
 // Don't include AdamsCoefficients.hpp in the header just to get one
 // constant.
 static_assert(adams_coefficients::maximum_order ==
-              AdamsMoultonPc::maximum_order);
+              AdamsMoultonPc<false>::maximum_order);
 
 namespace {
 template <typename T>
-void clean_history(const MutableUntypedHistory<T>& history) {
+void clean_history_nonmonotonic(const MutableUntypedHistory<T>& history) {
   ASSERT(history.integration_order() > 1, "Cannot run below second order.");
   const auto required_points = history.integration_order() - 1;
   ASSERT(history.size() >= required_points,
@@ -47,6 +48,44 @@ void clean_history(const MutableUntypedHistory<T>& history) {
   }
 }
 
+template <typename T>
+void clean_history_monotonic(const MutableUntypedHistory<T>& history) {
+  ASSERT(history.integration_order() > 1, "Cannot run below second order.");
+  // Unlike in nonmonotonic mode, monotonic cleaning happens after the
+  // update.  This behavior is reasonable because it is normally kept
+  // for dense output, which has already happened.
+  //
+  // It is also a bit of a hack, since it is necessary to trick the
+  // dense output code into using the correct starting values.  During
+  // dense output, the variables are initialized with
+  // History::complete_step_start(), but monotonic dense output
+  // doesn't require a complete step, so this is wrong.  Clearing the
+  // substeps tricks that code into considering the current step to be
+  // complete.
+  //
+  // Unfortunately, this also breaks self-start, which sometimes jumps
+  // out of the middle of the algorithm to avoid triggering the
+  // clean-up code to elongate the history.  This is dealt with by
+  // cleaning one fewer entries during self-start.
+  if (not history.at_step_start()) {
+    const auto required_points =
+        ::SelfStart::is_self_starting(history.back().time_step_id)
+            ? history.integration_order() - 1
+            : history.integration_order() - 2;
+    ASSERT(history.size() >= required_points,
+           "Insufficient data to take an order-" << history.integration_order()
+           << " step.  Have " << history.size() << " times, need "
+           << required_points);
+    history.clear_substeps();
+    while (history.size() > required_points) {
+      history.pop_front();
+    }
+    if (not history.empty()) {
+      history.discard_value(history.back().time_step_id);
+    }
+  }
+}
+
 template <typename T, typename TimeType>
 void update_u_common(const gsl::not_null<T*> u,
                      const ConstUntypedHistory<T>& history,
@@ -88,24 +127,39 @@ void update_u_common(const gsl::not_null<T*> u,
 }
 }  // namespace
 
-AdamsMoultonPc::AdamsMoultonPc(const size_t order) : order_(order) {
+template <bool Monotonic>
+AdamsMoultonPc<Monotonic>::AdamsMoultonPc(const size_t order) : order_(order) {
   ASSERT(order >= minimum_order and order <= maximum_order,
          "Invalid order: " << order);
 }
 
-size_t AdamsMoultonPc::order() const { return order_; }
+template <bool Monotonic>
+size_t AdamsMoultonPc<Monotonic>::order() const {
+  return order_;
+}
 
-size_t AdamsMoultonPc::error_estimate_order() const { return order_ - 1; }
+template <bool Monotonic>
+size_t AdamsMoultonPc<Monotonic>::error_estimate_order() const {
+  return order_ - 1;
+}
 
-uint64_t AdamsMoultonPc::number_of_substeps() const { return 2; }
+template <bool Monotonic>
+uint64_t AdamsMoultonPc<Monotonic>::number_of_substeps() const {
+  return 2;
+}
 
-uint64_t AdamsMoultonPc::number_of_substeps_for_error() const {
+template <bool Monotonic>
+uint64_t AdamsMoultonPc<Monotonic>::number_of_substeps_for_error() const {
   return number_of_substeps();
 }
 
-size_t AdamsMoultonPc::number_of_past_steps() const { return order_ - 2; }
+template <bool Monotonic>
+size_t AdamsMoultonPc<Monotonic>::number_of_past_steps() const {
+  return order_ - 2;
+}
 
-double AdamsMoultonPc::stable_step() const {
+template <bool Monotonic>
+double AdamsMoultonPc<Monotonic>::stable_step() const {
   switch (order_) {
     case 2:
       return 1.0;
@@ -126,8 +180,14 @@ double AdamsMoultonPc::stable_step() const {
   }
 }
 
-TimeStepId AdamsMoultonPc::next_time_id(const TimeStepId& current_id,
-                                        const TimeDelta& time_step) const {
+template <bool Monotonic>
+bool AdamsMoultonPc<Monotonic>::monotonic() const {
+  return Monotonic;
+}
+
+template <bool Monotonic>
+TimeStepId AdamsMoultonPc<Monotonic>::next_time_id(
+    const TimeStepId& current_id, const TimeDelta& time_step) const {
   switch (current_id.substep()) {
     case 0:
       return current_id.next_substep(time_step, 1.0);
@@ -138,62 +198,87 @@ TimeStepId AdamsMoultonPc::next_time_id(const TimeStepId& current_id,
   }
 }
 
-TimeStepId AdamsMoultonPc::next_time_id_for_error(
+template <bool Monotonic>
+TimeStepId AdamsMoultonPc<Monotonic>::next_time_id_for_error(
     const TimeStepId& current_id, const TimeDelta& time_step) const {
   return next_time_id(current_id, time_step);
 }
 
-void AdamsMoultonPc::pup(PUP::er& p) {
+template <bool Monotonic>
+void AdamsMoultonPc<Monotonic>::pup(PUP::er& p) {
   TimeStepper::pup(p);
   p | order_;
 }
 
+template <bool Monotonic>
 template <typename T>
-void AdamsMoultonPc::update_u_impl(const gsl::not_null<T*> u,
-                                   const MutableUntypedHistory<T>& history,
-                                   const TimeDelta& time_step) const {
-  clean_history(history);
+void AdamsMoultonPc<Monotonic>::update_u_impl(
+    const gsl::not_null<T*> u, const MutableUntypedHistory<T>& history,
+    const TimeDelta& time_step) const {
+  if constexpr (not Monotonic) {
+    clean_history_nonmonotonic(history);
+  }
   const Time next_time = history.back().time_step_id.step_time() + time_step;
   *u = *history.back().value;
   update_u_common(u, history, next_time, history.integration_order(),
                   not history.at_step_start());
+  if constexpr (Monotonic) {
+    clean_history_monotonic(history);
+  }
 }
 
+template <bool Monotonic>
 template <typename T>
-bool AdamsMoultonPc::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);
+bool AdamsMoultonPc<Monotonic>::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 {
+  if constexpr (not Monotonic) {
+    clean_history_nonmonotonic(history);
+  }
   const bool predictor = history.at_step_start();
   const Time next_time = history.back().time_step_id.step_time() + time_step;
   *u = *history.back().value;
   update_u_common(u, history, next_time, history.integration_order(),
                   not predictor);
   if (predictor) {
+    // No monotonic cleaning necessary on predictor phase.
     return false;
   }
   *u_error = *history.back().value;
   update_u_common(u_error, history, next_time, history.integration_order() - 1,
                   true);
   *u_error = *u - *u_error;
+  if constexpr (Monotonic) {
+    clean_history_monotonic(history);
+  }
   return true;
 }
 
+template <bool Monotonic>
 template <typename T>
-bool AdamsMoultonPc::dense_update_u_impl(const gsl::not_null<T*> u,
-                                         const ConstUntypedHistory<T>& history,
-                                         const double time) const {
-  if (history.at_step_start()) {
-    return false;
+bool AdamsMoultonPc<Monotonic>::dense_update_u_impl(
+    const gsl::not_null<T*> u, const ConstUntypedHistory<T>& history,
+    const double time) const {
+  if constexpr (Monotonic) {
+    if (not history.at_step_start()) {
+      return false;
+    }
+    update_u_common(u, history, ApproximateTime{time},
+                    history.integration_order(), false);
+    return true;
+  } else {
+    if (history.at_step_start()) {
+      return false;
+    }
+    update_u_common(u, history, ApproximateTime{time},
+                    history.integration_order(), true);
+    return true;
   }
-  update_u_common(u, history, ApproximateTime{time},
-                  history.integration_order(), true);
-  return true;
 }
 
+template <bool Monotonic>
 template <typename T>
-bool AdamsMoultonPc::can_change_step_size_impl(
+bool AdamsMoultonPc<Monotonic>::can_change_step_size_impl(
     const TimeStepId& time_id, const ConstUntypedHistory<T>& history) const {
   // 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
@@ -212,15 +297,34 @@ bool AdamsMoultonPc::can_change_step_size_impl(
          });
 }
 
-bool operator==(const AdamsMoultonPc& lhs, const AdamsMoultonPc& rhs) {
+template <bool Monotonic>
+bool operator==(const AdamsMoultonPc<Monotonic>& lhs,
+                const AdamsMoultonPc<Monotonic>& rhs) {
   return lhs.order() == rhs.order();
 }
 
-bool operator!=(const AdamsMoultonPc& lhs, const AdamsMoultonPc& rhs) {
+template <bool Monotonic>
+bool operator!=(const AdamsMoultonPc<Monotonic>& lhs,
+                const AdamsMoultonPc<Monotonic>& rhs) {
   return not(lhs == rhs);
 }
 
-TIME_STEPPER_DEFINE_OVERLOADS(AdamsMoultonPc)
-}  // namespace TimeSteppers
+TIME_STEPPER_DEFINE_OVERLOADS_TEMPLATED(AdamsMoultonPc<Monotonic>,
+                                        bool Monotonic)
 
-PUP::able::PUP_ID TimeSteppers::AdamsMoultonPc::my_PUP_ID = 0;  // NOLINT
+template <bool Monotonic>
+PUP::able::PUP_ID AdamsMoultonPc<Monotonic>::my_PUP_ID = 0;  // NOLINT
+
+#define MONOTONIC(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATE(_, data)                                            \
+  template class AdamsMoultonPc<MONOTONIC(data)>;                       \
+  template bool operator==(const AdamsMoultonPc<MONOTONIC(data)>& lhs,  \
+                           const AdamsMoultonPc<MONOTONIC(data)>& rhs); \
+  template bool operator!=(const AdamsMoultonPc<MONOTONIC(data)>& lhs,  \
+                           const AdamsMoultonPc<MONOTONIC(data)>& rhs);
+
+GENERATE_INSTANTIATIONS(INSTANTIATE, (false, true))
+#undef INSTANTIATE
+#undef MONOTONIC
+}  // namespace TimeSteppers

src/Time/TimeSteppers/AdamsMoultonPc.hpp

@@ -36,6 +36,12 @@ namespace TimeSteppers {
  * An \f$N\$th order Adams-Moulton predictor-corrector method using an
  * \$(N - 1)\$th order Adams-Bashforth predictor.
  *
+ * If \p Monotonic is true, dense output is performed using the
+ * predictor stage, otherwise the corrector is used.  The corrector
+ * results are more accurate (but still formally the same order), but
+ * require a RHS evaluation at the end of the step before dense output
+ * can be performed.
+ *
  * The stable step size factors for different orders are (to
  * approximately 4-5 digits):
  *
@@ -74,8 +80,13 @@ namespace TimeSteppers {
  *  </tr>
  * </table>
  */
+template <bool Monotonic>
 class AdamsMoultonPc : public TimeStepper {
  public:
+  static std::string name() {
+    return Monotonic ? "AdamsMoultonPcMonotonic" : "AdamsMoultonPc";
+  }
+
   static constexpr size_t minimum_order = 2;
   static constexpr size_t maximum_order = 8;
 
@@ -86,8 +97,11 @@ class AdamsMoultonPc : public TimeStepper {
     static type upper_bound() { return maximum_order; }
   };
   using options = tmpl::list<Order>;
-  static constexpr Options::String help = {
-      "An Adams-Moulton predictor-corrector time-stepper."};
+  static constexpr Options::String help =
+      Monotonic
+          ? "An Adams-Moulton predictor-corrector time-stepper with monotonic "
+            "dense output."
+          : "An Adams-Moulton predictor-corrector time-stepper.";
 
   AdamsMoultonPc() = default;
   explicit AdamsMoultonPc(size_t order);
@@ -109,6 +123,8 @@ class AdamsMoultonPc : public TimeStepper {
 
   double stable_step() const override;
 
+  bool monotonic() const override;
+
   TimeStepId next_time_id(const TimeStepId& current_id,
                           const TimeDelta& time_step) const override;
 
@@ -146,6 +162,10 @@ class AdamsMoultonPc : public TimeStepper {
   size_t order_{};
 };
 
-bool operator==(const AdamsMoultonPc& lhs, const AdamsMoultonPc& rhs);
-bool operator!=(const AdamsMoultonPc& lhs, const AdamsMoultonPc& rhs);
+template <bool Monotonic>
+bool operator==(const AdamsMoultonPc<Monotonic>& lhs,
+                const AdamsMoultonPc<Monotonic>& rhs);
+template <bool Monotonic>
+bool operator!=(const AdamsMoultonPc<Monotonic>& lhs,
+                const AdamsMoultonPc<Monotonic>& rhs);
 }  // namespace TimeSteppers