Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <boost/container/static_vector.hpp>
       7             : #include <cmath>
       8             : #include <cstddef>
       9             : #include <iterator>
      10             : #include <type_traits>
      11             : 
      12             : #include "Time/Time.hpp"
      13             : #include "Time/Utilities.hpp"
      14             : #include "Utilities/Algorithm.hpp"
      15             : 
      16             : /// \cond
      17             : struct ApproximateTime;
      18             : /// \endcond
      19             : 
      20             : /// Helpers for calculating Adams coefficients
      21           1 : namespace TimeSteppers::adams_coefficients {
      22           0 : constexpr size_t maximum_order = 8;
      23             : 
      24             : /// A vector holding one entry per order of integration.
      25             : template <typename T>
      26           1 : using OrderVector = boost::container::static_vector<T, maximum_order>;
      27             : 
      28             : /// The standard Adams-Bashforth coefficients for constant step size,
      29             : /// ordered from oldest to newest time, as one would find in a
      30             : /// reference table (except likely in the opposite order).
      31           1 : OrderVector<double> constant_adams_bashforth_coefficients(size_t order);
      32             : 
      33             : /// The standard Adams-Moulton coefficients for constant step size,
      34             : /// ordered from oldest to newest time, as one would find in a
      35             : /// reference table (except likely in the opposite order).
      36           1 : OrderVector<double> constant_adams_moulton_coefficients(size_t order);
      37             : 
      38             : /// \brief Generate coefficients for an Adams step.
      39             : ///
      40             : /// The coefficients are for a step using derivatives at \p
      41             : /// control_times, with the entries in the result vector corresponding
      42             : /// to the passed times in order.  The result includes the overall
      43             : /// factor of step size, so, for example, the coefficients for Euler's
      44             : /// method (`control_times = {0}, step_start=0, step_end=dt`) would be
      45             : /// `{dt}`, not `{1}`.
      46             : ///
      47             : /// No requirements are imposed on \p control_times, except that the
      48             : /// entries are all distinct.
      49             : ///
      50             : /// Only `T = double` is used by the time steppers, but `T = Rational`
      51             : /// can be used to generate coefficient tables.
      52             : template <typename T>
      53           1 : OrderVector<T> variable_coefficients(OrderVector<T> control_times,
      54             :                                      const T& step_start, const T& step_end);
      55             : 
      56             : /// \brief Get coefficients for a time step.
      57             : ///
      58             : /// Arguments are an iterator pair to past times (of type `Time`),
      59             : /// with the most recent last, and the start and end of the time step
      60             : /// to take, with the end a `Time` or `ApproximateTime`.  This
      61             : /// performs the same calculation as `variable_coefficients`, except
      62             : /// that it works with `Time`s and will detect and optimize the
      63             : /// constant-step-size case.
      64             : template <typename Iterator, typename TimeType>
      65           1 : OrderVector<double> coefficients(const Iterator& times_begin,
      66             :                                  const Iterator& times_end,
      67             :                                  const Time& step_start,
      68             :                                  const TimeType& step_end) {
      69             :   static_assert(std::is_same_v<TimeType, Time> or
      70             :                 std::is_same_v<TimeType, ApproximateTime>);
      71             :   if (times_begin == times_end) {
      72             :     return {};
      73             :   }
      74             :   const double step_size = (step_end - step_start).value();
      75             :   bool constant_step_size = true;
      76             :   // We shift the control times to be near zero, which gives smaller
      77             :   // errors from the variable_coefficients function.
      78             :   OrderVector<double> control_times{0.0};
      79             :   Time previous_time = *times_begin;
      80             :   for (auto t = std::next(times_begin); t != times_end; ++t) {
      81             :     const Time this_time = *t;
      82             :     const double this_step = (this_time - previous_time).value();
      83             :     control_times.push_back(control_times.back() + this_step);
      84             :     if (constant_step_size and
      85             :         std::abs(this_step - step_size) > slab_rounding_error(this_time)) {
      86             :       constant_step_size = false;
      87             :     }
      88             :     previous_time = this_time;
      89             :   }
      90             :   if (constant_step_size and step_start == previous_time) {
      91             :     auto result = constant_adams_bashforth_coefficients(control_times.size());
      92             :     alg::for_each(result, [&](double& coef) { coef *= step_size; });
      93             :     return result;
      94             :   } else if (constant_step_size and step_end == previous_time) {
      95             :     auto result = constant_adams_moulton_coefficients(control_times.size());
      96             :     alg::for_each(result, [&](double& coef) { coef *= step_size; });
      97             :     return result;
      98             :   }
      99             : 
     100             :   return variable_coefficients(
     101             :       control_times,
     102             :       control_times.back() + (step_start - previous_time).value(),
     103             :       control_times.back() + (step_end - previous_time).value());
     104             : }
     105             : }  // namespace TimeSteppers::adams_coefficients