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/SlabRoundingError.hpp"
      13              : #include "Time/Time.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