Add LTS support to AdamsMoultonPc

src/Time/TimeSteppers/AdamsLts.cpp

@@ -316,15 +316,6 @@ LtsCoefficients lts_coefficients(const ConstBoundaryHistoryTimes& local_times,
                                  const AdamsScheme& local_scheme,
                                  const AdamsScheme& remote_scheme,
                                  const AdamsScheme& small_step_scheme) {
-  ASSERT(
-      small_step_scheme == local_scheme or small_step_scheme == remote_scheme,
-      "FIXME just deduce?");
-  ASSERT(small_step_scheme == local_scheme or
-             small_step_scheme.order > local_scheme.order,
-         "FIXME just deduce?");
-  ASSERT(small_step_scheme == remote_scheme or
-             small_step_scheme.order > remote_scheme.order,
-         "FIXME just deduce?");
   if (start_time == end_time) {
     return {};
   }

src/Time/TimeSteppers/AdamsMoultonPc.cpp

@@ -13,6 +13,7 @@
 #include "Time/History.hpp"
 #include "Time/SelfStart.hpp"
 #include "Time/TimeSteppers/AdamsCoefficients.hpp"
+#include "Time/TimeSteppers/AdamsLts.hpp"
 #include "Utilities/Algorithm.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
@@ -145,9 +146,53 @@ TimeStepId AdamsMoultonPc<Monotonic>::next_time_id_for_error(
   return next_time_id(current_id, time_step);
 }
 
+template <bool Monotonic>
+bool AdamsMoultonPc<Monotonic>::neighbor_data_required(
+    const TimeStepId& next_substep_id,
+    const TimeStepId& neighbor_data_id) const {
+  // Because of self-start, step times may not be monotonic across
+  // slabs, so check that first.
+  if (neighbor_data_id.slab_number() != next_substep_id.slab_number()) {
+    return neighbor_data_id.slab_number() < next_substep_id.slab_number();
+  }
+
+  const evolution_less<Time> before{neighbor_data_id.time_runs_forward()};
+
+  if constexpr (Monotonic) {
+    const auto next_time = adams_lts::exact_substep_time(next_substep_id);
+    const auto neighbor_time = adams_lts::exact_substep_time(neighbor_data_id);
+    return before(neighbor_time, next_time) or
+           (neighbor_time == next_time and neighbor_data_id.substep() == 1 and
+            next_substep_id.substep() == 0);
+  } else {
+    if (next_substep_id.substep() == 1) {
+      // predictor
+      return before(neighbor_data_id.step_time(),
+                    next_substep_id.step_time()) or
+             (neighbor_data_id.step_time() == next_substep_id.step_time() and
+              neighbor_data_id.substep() == 0);
+    } else {
+      // corrector
+      return before(neighbor_data_id.step_time(), next_substep_id.step_time());
+    }
+  }
+}
+
+template <bool Monotonic>
+bool AdamsMoultonPc<Monotonic>::neighbor_data_required(
+    const double dense_output_time, const TimeStepId& neighbor_data_id) const {
+  const evolution_less<double> before{neighbor_data_id.time_runs_forward()};
+  if constexpr (Monotonic) {
+    return not before(dense_output_time,
+                      adams_lts::exact_substep_time(neighbor_data_id).value());
+  } else {
+    return before(neighbor_data_id.step_time().value(), dense_output_time);
+  }
+}
+
 template <bool Monotonic>
 void AdamsMoultonPc<Monotonic>::pup(PUP::er& p) {
-  TimeStepper::pup(p);
+  LtsTimeStepper::pup(p);
   p | order_;
 }
 
@@ -245,6 +290,206 @@ bool AdamsMoultonPc<Monotonic>::can_change_step_size_impl(
          });
 }
 
+template <bool Monotonic>
+template <typename T>
+void AdamsMoultonPc<Monotonic>::add_boundary_delta_impl(
+    const gsl::not_null<T*> result,
+    const TimeSteppers::ConstBoundaryHistoryTimes& local_times,
+    const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
+    const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
+    const TimeDelta& time_step) const {
+  ASSERT(not local_times.empty(), "No local data provided.");
+  ASSERT(not remote_times.empty(), "No remote data provided.");
+  const auto current_order =
+      local_times.integration_order(local_times.size() - 1);
+  if constexpr (Monotonic) {
+    const adams_lts::AdamsScheme predictor_scheme{
+        adams_lts::SchemeType::Explicit, current_order - 1};
+    const adams_lts::AdamsScheme corrector_scheme{
+        adams_lts::SchemeType::Implicit, current_order};
+    const auto step_start = local_times.back().step_time();
+    const auto step_end = step_start + time_step;
+    const auto small_step_start =
+        std::max(local_times.back(), remote_times.back()).step_time();
+    const auto synchronization_time =
+        std::min(local_times.back(), remote_times.back()).step_time();
+    const auto is_synchronization_time = [&](const TimeStepId& id) {
+      return id.step_time() == synchronization_time;
+    };
+    ASSERT(alg::any_of(local_times, is_synchronization_time) and
+               alg::any_of(remote_times, is_synchronization_time),
+           "Only nested step patterns (N:1) are supported.");
+
+    if (local_times.number_of_substeps(local_times.size() - 1) == 1) {
+      // Predictor
+      auto lts_coefficients = adams_lts::lts_coefficients(
+          local_times, remote_times, small_step_start, step_end,
+          predictor_scheme, predictor_scheme, predictor_scheme);
+      if (not is_synchronization_time(remote_times.back())) {
+        lts_coefficients += adams_lts::lts_coefficients(
+            local_times, remote_times, synchronization_time, small_step_start,
+            predictor_scheme, corrector_scheme, corrector_scheme);
+      }
+      adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+    } else {
+      // Corrector
+      if (remote_times.number_of_substeps(remote_times.size() - 1) == 2) {
+        // Aligned corrector
+        ASSERT(adams_lts::exact_substep_time(
+                   remote_times[{remote_times.size() - 1, 1}]) == step_end,
+               "Have remote substep data, but it isn't aligned with the local "
+               "data.");
+        const auto lts_coefficients =
+            adams_lts::lts_coefficients(
+                local_times, remote_times, synchronization_time, step_end,
+                corrector_scheme, corrector_scheme, corrector_scheme) -
+            adams_lts::lts_coefficients(
+                local_times, remote_times, synchronization_time, step_start,
+                corrector_scheme, predictor_scheme, corrector_scheme);
+        adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+      } else {
+        // Unaligned corrector
+        ASSERT(step_start == small_step_start,
+               "Trying to take unaligned step, but remote side is smaller.");
+        const auto lts_coefficients = adams_lts::lts_coefficients(
+            local_times, remote_times, step_start, step_end, corrector_scheme,
+            predictor_scheme, corrector_scheme);
+        adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+      }
+    }
+  } else {
+    adams_lts::AdamsScheme scheme{adams_lts::SchemeType::Implicit,
+                                  current_order};
+    auto remote_scheme = scheme;
+
+    if (local_times.number_of_substeps(local_times.size() - 1) == 1) {
+      // Predictor
+      scheme = {adams_lts::SchemeType::Explicit, current_order - 1};
+      ASSERT(remote_times.back() <= local_times.back(),
+             "Unexpected remote values available.");
+      // If the sides are not aligned, we use the predictor data
+      // available from the neighbor.  If they are, that data has not
+      // been received.
+      remote_scheme = {remote_times.back() == local_times.back()
+                           ? adams_lts::SchemeType::Explicit
+                           : adams_lts::SchemeType::Implicit,
+                       current_order - 1};
+    }
+
+    const auto lts_coefficients = adams_lts::lts_coefficients(
+        local_times, remote_times, local_times.back().step_time(),
+        local_times.back().step_time() + time_step, scheme, remote_scheme,
+        scheme);
+    adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+  }
+}
+
+template <bool Monotonic>
+void AdamsMoultonPc<Monotonic>::clean_boundary_history_impl(
+    const TimeSteppers::MutableBoundaryHistoryTimes& local_times,
+    const TimeSteppers::MutableBoundaryHistoryTimes& remote_times) const {
+  if (local_times.empty() or
+      local_times.number_of_substeps(local_times.size() - 1) != 2) {
+    return;
+  }
+  ASSERT(not remote_times.empty(), "No remote data available.");
+
+  const bool synchronized =
+      remote_times.number_of_substeps(remote_times.size() - 1) == 2 and
+      adams_lts::exact_substep_time(local_times[{local_times.size() - 1, 1}]) ==
+          adams_lts::exact_substep_time(
+              remote_times[{remote_times.size() - 1, 1}]);
+
+  if (Monotonic and not synchronized) {
+    return;
+  }
+
+  const auto required_points =
+      local_times.integration_order(local_times.size() - 1) - 2;
+
+  while (local_times.size() > required_points) {
+    local_times.pop_front();
+  }
+  for (size_t i = 0; i < local_times.size(); ++i) {
+    local_times.clear_substeps(i);
+  }
+
+  // If the sides are not aligned, then we are in the middle of the
+  // remote step, so still need its data.
+  if (synchronized) {
+    while (remote_times.size() > required_points) {
+      remote_times.pop_front();
+    }
+    for (size_t i = 0; i < remote_times.size(); ++i) {
+      remote_times.clear_substeps(i);
+    }
+  }
+}
+
+template <bool Monotonic>
+template <typename T>
+void AdamsMoultonPc<Monotonic>::boundary_dense_output_impl(
+    const gsl::not_null<T*> result,
+    const TimeSteppers::ConstBoundaryHistoryTimes& local_times,
+    const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
+    const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
+    const double time) const {
+  if constexpr (Monotonic) {
+    ASSERT(local_times.number_of_substeps(local_times.size() - 1) == 1,
+           "Dense output must be done before predictor evaluation.");
+
+    const auto current_order =
+        local_times.integration_order(local_times.size() - 1);
+    const adams_lts::AdamsScheme predictor_scheme{
+        adams_lts::SchemeType::Explicit, current_order - 1};
+    const adams_lts::AdamsScheme corrector_scheme{
+        adams_lts::SchemeType::Implicit, current_order};
+    const auto small_step_start =
+        std::max(local_times.back(), remote_times.back()).step_time();
+    const auto synchronization_time =
+        std::min(local_times.back(), remote_times.back()).step_time();
+    const auto is_synchronization_time = [&](const TimeStepId& id) {
+      return id.step_time() == synchronization_time;
+    };
+    ASSERT(alg::any_of(local_times, is_synchronization_time) and
+               alg::any_of(remote_times, is_synchronization_time),
+           "Only nested step patterns (N:1) are supported.");
+
+    auto lts_coefficients = adams_lts::lts_coefficients(
+        local_times, remote_times, small_step_start, ApproximateTime{time},
+        predictor_scheme, predictor_scheme, predictor_scheme);
+    if (not is_synchronization_time(remote_times.back())) {
+      lts_coefficients += adams_lts::lts_coefficients(
+          local_times, remote_times, synchronization_time, small_step_start,
+          predictor_scheme, corrector_scheme, corrector_scheme);
+    }
+    adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+  } else {
+    if (local_times.back().step_time().value() == time) {
+      // Nothing to do.  The requested time is the start of the step,
+      // which is the input value of `result`.
+      return;
+    }
+    ASSERT(local_times.number_of_substeps(local_times.size() - 1) == 2,
+           "Dense output must be done after predictor evaluation.");
+
+    const auto current_order =
+        local_times.integration_order(local_times.size() - 1);
+    const adams_lts::AdamsScheme scheme{adams_lts::SchemeType::Implicit,
+                                        current_order};
+    const auto small_step_start =
+        std::max(local_times.back(), remote_times.back()).step_time();
+    const auto lts_coefficients =
+        adams_lts::lts_coefficients(local_times, remote_times,
+                                    local_times.back().step_time(),
+                                    small_step_start, scheme, scheme, scheme) +
+        adams_lts::lts_coefficients(local_times, remote_times, small_step_start,
+                                    ApproximateTime{time}, scheme, scheme,
+                                    scheme);
+    adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+  }
+}
+
 template <bool Monotonic>
 bool operator==(const AdamsMoultonPc<Monotonic>& lhs,
                 const AdamsMoultonPc<Monotonic>& rhs) {
@@ -259,6 +504,8 @@ bool operator!=(const AdamsMoultonPc<Monotonic>& lhs,
 
 TIME_STEPPER_DEFINE_OVERLOADS_TEMPLATED(AdamsMoultonPc<Monotonic>,
                                         bool Monotonic)
+LTS_TIME_STEPPER_DEFINE_OVERLOADS_TEMPLATED(AdamsMoultonPc<Monotonic>,
+                                            bool Monotonic)
 
 template <bool Monotonic>
 PUP::able::PUP_ID AdamsMoultonPc<Monotonic>::my_PUP_ID = 0;  // NOLINT

src/Time/TimeSteppers/AdamsMoultonPc.hpp

@@ -7,7 +7,7 @@
 #include <cstdint>
 
 #include "Options/String.hpp"
-#include "Time/TimeSteppers/TimeStepper.hpp"
+#include "Time/TimeSteppers/LtsTimeStepper.hpp"
 #include "Utilities/Serialization/CharmPupable.hpp"
 #include "Utilities/TMPL.hpp"
 
@@ -19,7 +19,11 @@ class er;
 }  // namespace PUP
 namespace TimeSteppers {
 template <typename T>
+class BoundaryHistoryEvaluator;
+class ConstBoundaryHistoryTimes;
+template <typename T>
 class ConstUntypedHistory;
+class MutableBoundaryHistoryTimes;
 template <typename T>
 class MutableUntypedHistory;
 }  // namespace TimeSteppers
@@ -81,7 +85,7 @@ namespace TimeSteppers {
  * </table>
  */
 template <bool Monotonic>
-class AdamsMoultonPc : public TimeStepper {
+class AdamsMoultonPc : public LtsTimeStepper {
  public:
   static std::string name() {
     return Monotonic ? "AdamsMoultonPcMonotonic" : "AdamsMoultonPc";
@@ -131,6 +135,14 @@ class AdamsMoultonPc : public TimeStepper {
   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(AdamsMoultonPc);  // NOLINT
 
   explicit AdamsMoultonPc(CkMigrateMessage* /*unused*/) {}
@@ -159,7 +171,29 @@ class AdamsMoultonPc : public TimeStepper {
   bool can_change_step_size_impl(const TimeStepId& time_id,
                                  const ConstUntypedHistory<T>& history) const;
 
+  template <typename T>
+  void add_boundary_delta_impl(
+      gsl::not_null<T*> result,
+      const TimeSteppers::ConstBoundaryHistoryTimes& local_times,
+      const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
+      const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
+      const TimeDelta& time_step) const;
+
+  void clean_boundary_history_impl(
+      const TimeSteppers::MutableBoundaryHistoryTimes& local_times,
+      const TimeSteppers::MutableBoundaryHistoryTimes& remote_times)
+      const override;
+
+  template <typename T>
+  void boundary_dense_output_impl(
+      gsl::not_null<T*> result,
+      const TimeSteppers::ConstBoundaryHistoryTimes& local_times,
+      const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
+      const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
+      double time) const;
+
   TIME_STEPPER_DECLARE_OVERLOADS
+  LTS_TIME_STEPPER_DECLARE_OVERLOADS
 
   size_t order_{};
 };

src/Time/TimeSteppers/Factory.hpp

@@ -27,7 +27,8 @@ using time_steppers =
                Rk5Owren, Rk5Tsitouras>;
 
 /// Typelist of available LtsTimeSteppers
-using lts_time_steppers = tmpl::list<AdamsBashforth>;
+using lts_time_steppers =
+    tmpl::list<AdamsBashforth, AdamsMoultonPc<false>, AdamsMoultonPc<true>>;
 
 /// Typelist of available ImexTimeSteppers
 using imex_time_steppers =