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