Merge pull request #6565 from wthrowe/variable_order_coeffs

Add coefficient calculations for variable-order Adams LTS

Author: wthrowe

src/Time/TimeSteppers/AdamsBashforth.cpp

@@ -235,31 +235,41 @@ void AdamsBashforth::add_boundary_delta_impl(
     const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
     const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
     const TimeDelta& time_step) const {
-  const auto step_start = local_times.back().step_time();
-  const size_t integration_order =
-      local_times.integration_order(local_times.size() - 1);
+  for (size_t i = 0; i < local_times.size(); ++i) {
+    ASSERT(local_times.integration_order(i) == order_ or
+               ::SelfStart::is_self_starting(local_times[i]),
+           "Incorrect local order " << local_times.integration_order(i)
+           << " at time " << local_times[i]);
+  }
+  for (size_t i = 0; i < remote_times.size(); ++i) {
+    ASSERT(remote_times.integration_order(i) == order_ or
+               ::SelfStart::is_self_starting(remote_times[i]),
+           "Incorrect remote order " << remote_times.integration_order(i)
+           << " at time " << remote_times[i]);
+  }
 
-  const adams_lts::AdamsScheme scheme{adams_lts::SchemeType::Explicit,
-                                      integration_order};
+  const auto step_start = local_times.back().step_time();
   const auto lts_coefficients = adams_lts::lts_coefficients(
-      local_times, remote_times, step_start, step_start + time_step, scheme,
-      scheme, scheme);
+      local_times, remote_times, step_start, step_start + time_step,
+      adams_lts::SchemeType::Explicit, adams_lts::SchemeType::Explicit, 0, 0);
   adams_lts::apply_coefficients(result, lts_coefficients, coupling);
 }
 
 void AdamsBashforth::clean_boundary_history_impl(
     const TimeSteppers::MutableBoundaryHistoryTimes& local_times,
     const TimeSteppers::MutableBoundaryHistoryTimes& remote_times) const {
-  const size_t integration_order =
+  const size_t local_order =
       local_times.integration_order(local_times.size() - 1);
-
-  while (local_times.size() >= integration_order) {
+  while (local_times.size() >= local_order) {
     local_times.pop_front();
   }
+
+  const size_t remote_order =
+      remote_times.integration_order(remote_times.size() - 1);
   // We're guaranteed to have a new local value inserted before the
   // next use, but not a new remote value, so we need to keep one more
   // of these.
-  while (remote_times.size() > integration_order) {
+  while (remote_times.size() > remote_order) {
     remote_times.pop_front();
   }
 }
@@ -271,24 +281,35 @@ void AdamsBashforth::boundary_dense_output_impl(
     const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
     const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
     const double time) const {
+  for (size_t i = 0; i < local_times.size(); ++i) {
+    ASSERT(local_times.integration_order(i) == order_ or
+               ::SelfStart::is_self_starting(local_times[i]),
+           "Incorrect local order " << local_times.integration_order(i)
+           << " at time " << local_times[i]);
+  }
+  for (size_t i = 0; i < remote_times.size(); ++i) {
+    ASSERT(remote_times.integration_order(i) == order_ or
+               ::SelfStart::is_self_starting(remote_times[i]),
+           "Incorrect remote order " << remote_times.integration_order(i)
+           << " at time " << remote_times[i]);
+  }
+
   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;
   }
-  const auto current_order =
-      local_times.integration_order(local_times.size() - 1);
-  const adams_lts::AdamsScheme scheme{adams_lts::SchemeType::Explicit,
-                                      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, local_times.back().step_time(),
+          small_step_start, adams_lts::SchemeType::Explicit,
+          adams_lts::SchemeType::Explicit, 0, 0) +
       adams_lts::lts_coefficients(local_times, remote_times, small_step_start,
-                                  ApproximateTime{time}, scheme, scheme,
-                                  scheme);
+                                  ApproximateTime{time},
+                                  adams_lts::SchemeType::Explicit,
+                                  adams_lts::SchemeType::Explicit, 0, 0);
   adams_lts::apply_coefficients(result, lts_coefficients, coupling);
 }
 

src/Time/TimeSteppers/AdamsLts.cpp

@@ -155,13 +155,6 @@ void apply_coefficients(const gsl::not_null<T*> result,
   }
 }
 
-bool operator==(const AdamsScheme& a, const AdamsScheme& b) {
-  return a.type == b.type and a.order == b.order;
-}
-bool operator!=(const AdamsScheme& a, const AdamsScheme& b) {
-  return not(a == b);
-}
-
 namespace {
 // Collect the ids used for interpolating during a step to `end_time`
 // from `times`.
@@ -171,29 +164,33 @@ namespace {
 template <typename TimeType>
 OrderVector<TimeStepId> find_relevant_ids(
     const ConstBoundaryHistoryTimes& times, const TimeType& end_time,
-    const AdamsScheme& scheme) {
+    const SchemeType scheme_type, const int order_offset) {
   OrderVector<TimeStepId> ids{};
   using difference_type = std::iterator_traits<
       ConstBoundaryHistoryTimes::const_iterator>::difference_type;
   const evolution_less<> less{times.front().time_runs_forward()};
   // Can't do a binary search because times are not sorted during self-start.
   auto used_range_end = times.end();
-  while (used_range_end != times.begin()) {
+  for (;;) {
+    ASSERT(used_range_end != times.begin(),
+           "All times in history were after " << end_time);
     if (less((used_range_end - 1)->step_time(), end_time)) {
       break;
     }
     --used_range_end;
   }
+  const auto last_step =
+      static_cast<size_t>(used_range_end - times.begin() - 1);
+  const auto order =
+      static_cast<difference_type>(times.integration_order(last_step)) +
+      order_offset;
   const auto number_of_past_steps =
-      static_cast<difference_type>(scheme.order) -
-      (scheme.type == SchemeType::Implicit ? 1 : 0);
+      order - (scheme_type == SchemeType::Implicit ? 1 : 0);
   ASSERT(used_range_end - times.begin() >= number_of_past_steps,
          "Insufficient past data.");
   std::copy(used_range_end - number_of_past_steps, used_range_end,
             std::back_inserter(ids));
-  if (scheme.type == SchemeType::Implicit) {
-    const auto last_step =
-        static_cast<size_t>(used_range_end - times.begin() - 1);
+  if (scheme_type == SchemeType::Implicit) {
     ASSERT(times.number_of_substeps(last_step) == 2,
            "Must have substep data for implicit stepping.");
     ids.push_back(times[{last_step, 1}]);
@@ -202,52 +199,42 @@ OrderVector<TimeStepId> find_relevant_ids(
 }
 
 // Choose the relevant times from `local` and `remote` for defining
-// small steps for `small_step_scheme`, using the specified schemes
-// for interpolation on the local and remote sides.
+// the small steps, using the specified schemes for interpolation on
+// the local and remote sides.
 //
-// This is the most recent values from the union of the `local` and
-// `remote` times with any values that should only be used for
-// interpolation removed.
+// This is the `num_steps` most recent values from the union of the
+// `local` and `remote` times, excluding any values that should only
+// be used for interpolation.
 OrderVector<Time> merge_to_small_steps(const OrderVector<Time>& local,
                                        const OrderVector<Time>& remote,
                                        const evolution_less<Time>& less,
-                                       const AdamsScheme& local_scheme,
-                                       const AdamsScheme& remote_scheme,
-                                       const AdamsScheme& small_step_scheme) {
-  OrderVector<Time> small_steps(small_step_scheme.order);
+                                       const SchemeType local_scheme,
+                                       const SchemeType remote_scheme,
+                                       const size_t num_steps) {
+  OrderVector<Time> small_steps(num_steps);
   auto local_it = local.rbegin();
   auto remote_it = remote.rbegin();
 
-  ASSERT(not(local_scheme.type == SchemeType::Implicit and
-             remote_scheme.type == SchemeType::Explicit and
+  ASSERT(not(local_scheme == SchemeType::Implicit and
+             remote_scheme == SchemeType::Explicit and
              less(*local_it, *remote_it)),
          "Explicit time " << *remote_it << " after implicit " << *local_it);
-  ASSERT(not(remote_scheme.type == SchemeType::Implicit and
-             local_scheme.type == SchemeType::Explicit and
+  ASSERT(not(remote_scheme == SchemeType::Implicit and
+             local_scheme == SchemeType::Explicit and
              less(*remote_it, *local_it)),
          "Explicit time " << *local_it << " after implicit " << *remote_it);
 
-  if (small_step_scheme.type == SchemeType::Explicit) {
-    // Don't use implicit interpolation points for an explicit step
-    if (local_scheme.type == SchemeType::Implicit) {
-      ++local_it;
-    }
-    if (remote_scheme.type == SchemeType::Implicit) {
+  if (local_scheme == SchemeType::Implicit and
+      remote_scheme == SchemeType::Implicit) {
+    // If both the interpolation schemes are implicit, we will get
+    // two times after the small step we are working on, one from
+    // each.  One of them (if they are different) belongs to a later
+    // small step, and we should ignore it.  If they are the same,
+    // ignoring one of them is harmless.
+    if (less(*local_it, *remote_it)) {
       ++remote_it;
-    }
-  } else {
-    if (local_scheme.type == SchemeType::Implicit and
-        remote_scheme.type == SchemeType::Implicit) {
-      // If both the interpolation schemes are implicit, we will get
-      // two times after the small step we are working on, one from
-      // each.  One of them (if they are different) belongs to a later
-      // small step, and we should ignore it.  If they are the same,
-      // ignoring one of them is harmless.
-      if (less(*local_it, *remote_it)) {
-        ++remote_it;
-      } else {
-        ++local_it;
-      }
+    } else {
+      ++local_it;
     }
   }
 
@@ -327,44 +314,61 @@ LtsCoefficients lts_coefficients(const ConstBoundaryHistoryTimes& local_times,
                                  const ConstBoundaryHistoryTimes& remote_times,
                                  const Time& start_time,
                                  const TimeType& end_time,
-                                 const AdamsScheme& local_scheme,
-                                 const AdamsScheme& remote_scheme,
-                                 const AdamsScheme& small_step_scheme) {
+                                 const SchemeType local_scheme,
+                                 const SchemeType remote_scheme,
+                                 const int local_order_offset,
+                                 const int remote_order_offset) {
+  ASSERT(local_order_offset == 0 or local_order_offset == -1,
+         "Must be 0 or -1, not " << local_order_offset);
+  ASSERT(remote_order_offset == 0 or remote_order_offset == -1,
+         "Must be 0 or -1, not " << remote_order_offset);
   if (start_time == end_time) {
     return {};
   }
+  const int integration_order_offset =
+      std::max(local_order_offset, remote_order_offset);
   const evolution_less<Time> time_less{local_times.front().time_runs_forward()};
 
   LtsCoefficients step_coefficients{};
 
   TimeType small_step_end = end_time;
   for (;;) {
-    const OrderVector<TimeStepId> local_ids =
-        find_relevant_ids(local_times, small_step_end, local_scheme);
-    const OrderVector<TimeStepId> remote_ids =
-        find_relevant_ids(remote_times, small_step_end, remote_scheme);
+    const OrderVector<TimeStepId> local_ids = find_relevant_ids(
+        local_times, small_step_end, local_scheme, local_order_offset);
+    const OrderVector<TimeStepId> remote_ids = find_relevant_ids(
+        remote_times, small_step_end, remote_scheme, remote_order_offset);
 
     // Check is the there is actually local time-stepping happening at
     // this boundary.  Only check for the latest small step, before we
     // have generated any coefficients.
-    if (step_coefficients.empty() and small_step_scheme == local_scheme and
-        small_step_scheme == remote_scheme and local_ids == remote_ids) {
+    if (step_coefficients.empty() and local_scheme == remote_scheme and
+        local_ids == remote_ids) {
       return lts_coefficients_for_gts(local_ids, start_time, end_time);
     }
 
-    OrderVector<Time> local_control_times(local_scheme.order);
+    OrderVector<Time> local_control_times(local_ids.size());
     alg::transform(local_ids, local_control_times.begin(), exact_substep_time);
-    OrderVector<Time> remote_control_times(remote_scheme.order);
+    OrderVector<Time> remote_control_times(remote_ids.size());
     alg::transform(remote_ids, remote_control_times.begin(),
                    exact_substep_time);
 
+    const size_t local_order =
+        local_ids.size() - static_cast<size_t>(local_order_offset);
+    const size_t remote_order =
+        remote_ids.size() - static_cast<size_t>(remote_order_offset);
+    const size_t integration_order =
+        std::max(local_order, remote_order) +
+        static_cast<size_t>(integration_order_offset);
+
     const OrderVector<Time> small_step_times = merge_to_small_steps(
         local_control_times, remote_control_times, time_less, local_scheme,
-        remote_scheme, small_step_scheme);
+        remote_scheme, integration_order);
     const Time current_small_step =
         small_step_times[small_step_times.size() -
-                         (small_step_scheme.type == SchemeType::Implicit ? 2
-                                                                         : 1)];
+                         (local_scheme == SchemeType::Implicit or
+                                  remote_scheme == SchemeType::Implicit
+                              ? 2
+                              : 1)];
     ASSERT(not time_less(current_small_step, start_time),
            "Reached time " << current_small_step
            << " without hitting start time " << start_time
@@ -459,8 +463,9 @@ GENERATE_INSTANTIATIONS(INSTANTIATE, (MATH_WRAPPER_TYPES))
   template LtsCoefficients lts_coefficients(                                 \
       const ConstBoundaryHistoryTimes& local_times,                          \
       const ConstBoundaryHistoryTimes& remote_times, const Time& start_time, \
-      const TIME_TYPE(data) & end_time, const AdamsScheme& local_scheme,     \
-      const AdamsScheme& remote_scheme, const AdamsScheme& small_step_scheme);
+      const TIME_TYPE(data) & end_time, SchemeType local_scheme,             \
+      SchemeType remote_scheme, int local_order_offset,                      \
+      int remote_order_offset);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (Time, ApproximateTime))
 #undef INSTANTIATE

src/Time/TimeSteppers/AdamsLts.hpp

@@ -69,24 +69,22 @@ void apply_coefficients(gsl::not_null<T*> result,
                         const LtsCoefficients& coefficients,
                         const BoundaryHistoryEvaluator<T>& coupling);
 
-/// Type of coefficients used for an AdamsScheme
+/// Type of coefficients used for an Adams scheme
 enum class SchemeType { Explicit, Implicit };
 
-struct AdamsScheme {
-  SchemeType type;
-  size_t order;
-};
-
-bool operator==(const AdamsScheme& a, const AdamsScheme& b);
-bool operator!=(const AdamsScheme& a, const AdamsScheme& b);
-
 /*!
  * Calculate the nonzero terms in an Adams LTS boundary contribution.
  *
  * The coefficients are generated for a step from \p start_time to \p
- * end_time, integrating using the \p small_step_scheme.
- * Interpolation from the elements is performed using polynomials
- * appropriate for the given \p local_scheme and \p remote_scheme.
+ * end_time.  Interpolation from the elements is performed using
+ * polynomials appropriate for the given \p local_scheme and \p
+ * remote_scheme.  The orders for the interpolations are taken from
+ * the \p local_times and \p remote_times, adjusted by \p
+ * local_order_offset and \p remote_order_offset, which can be either
+ * 0 or -1.  The small-step integration is performed with an implicit
+ * method if either \p local_scheme or \p remote_scheme is implicit,
+ * and the order is the larger of the integration orders on the two
+ * sides, reduced by one if both sides have offset -1.
  *
  * This function returns the sum of the coefficients for the small
  * steps, i.e., the steps between all the times either side updates
@@ -104,10 +102,11 @@ bool operator!=(const AdamsScheme& a, const AdamsScheme& b);
  * $n$, $y^L_{n,1}$ and $y^R_{n,1}$ are the predictor values for those
  * steps, and other index values count backwards in the history of
  * that side.  The range of $i$ and choice of the control times
- * $t^L_\cdots$ are determined by \p local_scheme, $j$ and
- * $t^R_\cdots$ are determined by \p remote_scheme, and the Adams
- * coefficients $\tilde{\alpha}_{nq}$ correspond to \p
- * small_step_scheme.
+ * $t^L_\cdots$ are determined by \p local_scheme and the local order,
+ * $j$ and $t^R_\cdots$ are determined by \p remote_scheme and the
+ * remote order, and the Adams coefficients $\tilde{\alpha}_{nq}$
+ * correspond to the larger of the two orders, and are implicit if
+ * either set of interpolation coefficients is.
  *
  * When called for dense output, the arguments must represent a
  * single small step, i.e., the step cannot cross any of the control
@@ -118,11 +117,9 @@ bool operator!=(const AdamsScheme& a, const AdamsScheme& b);
  * control times or `ApproximateTime` for dense output.
  */
 template <typename TimeType>
-LtsCoefficients lts_coefficients(const ConstBoundaryHistoryTimes& local_times,
-                                 const ConstBoundaryHistoryTimes& remote_times,
-                                 const Time& start_time,
-                                 const TimeType& end_time,
-                                 const AdamsScheme& local_scheme,
-                                 const AdamsScheme& remote_scheme,
-                                 const AdamsScheme& small_step_scheme);
+LtsCoefficients lts_coefficients(
+    const ConstBoundaryHistoryTimes& local_times,
+    const ConstBoundaryHistoryTimes& remote_times, const Time& start_time,
+    const TimeType& end_time, SchemeType local_scheme, SchemeType remote_scheme,
+    int local_order_offset, int remote_order_offset);
 }  // namespace TimeSteppers::adams_lts