SpECTRE Documentation Coverage Report
Current view: top level - Time/TimeSteppers - AdamsBashforthN.cpp Hit Total Coverage
Commit: b1342d46f40e2d46bbd11d0cef68fd973031a24b Lines: 0 3 0.0 %
Date: 2020-09-24 20:24:42
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #include "Time/TimeSteppers/AdamsBashforthN.hpp"
       5             : 
       6             : #include <algorithm>
       7             : 
       8             : #include "Time/TimeStepId.hpp"
       9             : 
      10             : namespace TimeSteppers {
      11             : 
      12             : AdamsBashforthN::AdamsBashforthN(const size_t order) noexcept : order_(order) {
      13             :   if (order_ < 1 or order_ > maximum_order) {
      14             :     ERROR("The order for Adams-Bashforth Nth order must be 1 <= order <= "
      15             :           << maximum_order);
      16             :   }
      17             : }
      18             : 
      19             : size_t AdamsBashforthN::number_of_past_steps() const noexcept {
      20             :   return order_ - 1;
      21             : }
      22             : 
      23             : double AdamsBashforthN::stable_step() const noexcept {
      24             :   if (order_ == 1) {
      25             :     return 1.;
      26             :   }
      27             : 
      28             :   // This is the condition that the characteristic polynomial of the
      29             :   // recurrence relation defined by the method has the correct sign at
      30             :   // -1.  It is not clear whether this is actually sufficient.
      31             :   const auto& coefficients = constant_coefficients(order_);
      32             :   double invstep = 0.;
      33             :   double sign = 1.;
      34             :   for (const auto coef : coefficients) {
      35             :     invstep += sign * coef;
      36             :     sign = -sign;
      37             :   }
      38             :   return 1. / invstep;
      39             : }
      40             : 
      41             : TimeStepId AdamsBashforthN::next_time_id(
      42             :     const TimeStepId& current_id,
      43             :     const TimeDelta& time_step) const noexcept {
      44             :   ASSERT(current_id.substep() == 0, "Adams-Bashforth should not have substeps");
      45             :   return {current_id.time_runs_forward(), current_id.slab_number(),
      46             :           current_id.step_time() + time_step};
      47             : }
      48             : 
      49             : std::vector<double> AdamsBashforthN::get_coefficients_impl(
      50             :     const std::vector<double>& steps) noexcept {
      51             :   const size_t order = steps.size();
      52             :   ASSERT(order >= 1 and order <= maximum_order, "Bad order" << order);
      53             :   if (std::all_of(steps.begin(), steps.end(),
      54             :                   [=](const double s) { return s == 1.; })) {
      55             :     return constant_coefficients(order);
      56             :   }
      57             : 
      58             :   return variable_coefficients(steps);
      59             : }
      60             : 
      61             : std::vector<double> AdamsBashforthN::variable_coefficients(
      62             :     const std::vector<double>& steps) noexcept {
      63             :   const size_t order = steps.size();  // "k" in below equations
      64             : 
      65             :   // The `steps` vector contains the relative step sizes:
      66             :   //   steps = {dt_{n-k+1}/dt_n, ..., dt_n/dt_n}
      67             :   // Our goal is to calculate, for each j, the coefficient given by
      68             :   //   \int_0^1 dt ell_j(t; 1, (dt_n + dt_{n-1})/dt_n, ...,
      69             :   //                        (dt_n + ... + dt_{n-k+1})/dt_n)
      70             :   // (Where the ell_j are the Lagrange interpolating polynomials.)
      71             : 
      72             :   // Calculate coefficients of the numerators of the Lagrange interpolating
      73             :   // polynomial, in the standard form.
      74             :   std::vector<std::vector<double>> polynomials(order,
      75             :                                                std::vector<double>(order, 0.));
      76             :   for (auto& poly : polynomials) {
      77             :     poly[0] = 1.;
      78             :   }
      79             :   {
      80             :     double step_sum = 0.;
      81             :     for (size_t m = 0; m < order; ++m) {
      82             :       const double step = steps[order - m - 1];
      83             :       step_sum += step;
      84             :       for (size_t j = 0; j < order; ++j) {
      85             :         if (m == j) {
      86             :           continue;
      87             :         }
      88             :         auto& poly = polynomials[j];
      89             :         for (size_t i = m + (m > j ? 0 : 1); i > 0; --i) {
      90             :           poly[i] = poly[i - 1] - step_sum * poly[i];
      91             :         }
      92             :         poly[0] *= -step_sum;
      93             :       }
      94             :     }
      95             :   }
      96             : 
      97             :   // Calculate the denominators of the Lagrange interpolating polynomials.
      98             :   std::vector<double> denominators;
      99             :   denominators.reserve(order);
     100             :   for (size_t j = 0; j < order; ++j) {
     101             :     double denom = 1.;
     102             :     double step_sum = 0.;
     103             :     for (size_t m = 0; m < j; ++m) {
     104             :       const double step = steps[order - j + m - 1];
     105             :       step_sum += step;
     106             :       denom *= step_sum;
     107             :     }
     108             :     step_sum = 0.;
     109             :     for (size_t m = 0; m < order - j - 1; ++m) {
     110             :       const double step = steps[order - j - m - 2];
     111             :       step_sum += step;
     112             :       denom *= step_sum;
     113             :     }
     114             :     denominators.push_back(denom);
     115             :   }
     116             : 
     117             :   // At this point, the Lagrange interpolating polynomials are given by:
     118             :   //   ell_j(t; ...) = +/- sum_m t^m polynomials[j][m] / denominators[j]
     119             : 
     120             :   // Integrate, term by term.
     121             :   std::vector<double> result;
     122             :   result.reserve(order);
     123             :   double overall_sign = order % 2 == 0 ? -1. : 1.;
     124             :   for (size_t j = 0; j < order; ++j) {
     125             :     const auto& poly = polynomials[j];
     126             :     double integral = 0.;
     127             :     for (size_t i = 0; i < order; ++i) {
     128             :       integral += poly[i] / (i + 1);
     129             :     }
     130             :     result.push_back(overall_sign * integral / denominators[j]);
     131             :     overall_sign = -overall_sign;
     132             :   }
     133             :   return result;
     134             : }
     135             : 
     136             : std::vector<double> AdamsBashforthN::constant_coefficients(
     137             :     const size_t order) noexcept {
     138             :   switch (order) {
     139             :     case 1: return {1.};
     140             :     case 2: return {1.5, -0.5};
     141             :     case 3: return {23.0 / 12.0, -4.0 / 3.0, 5.0 / 12.0};
     142             :     case 4: return {55.0 / 24.0, -59.0 / 24.0, 37.0 / 24.0, -3.0 / 8.0};
     143             :     case 5: return {1901.0 / 720.0, -1387.0 / 360.0, 109.0 / 30.0,
     144             :           -637.0 / 360.0, 251.0 / 720.0};
     145             :     case 6: return {4277.0 / 1440.0, -2641.0 / 480.0, 4991.0 / 720.0,
     146             :           -3649.0 / 720.0, 959.0 / 480.0, -95.0 / 288.0};
     147             :     case 7: return {198721.0 / 60480.0, -18637.0 / 2520.0, 235183.0 / 20160.0,
     148             :           -10754.0 / 945.0, 135713.0 / 20160.0, -5603.0 / 2520.0,
     149             :           19087.0 / 60480.0};
     150             :     case 8: return {16083.0 / 4480.0, -1152169.0 / 120960.0, 242653.0 / 13440.0,
     151             :           -296053.0 / 13440.0, 2102243.0 / 120960.0, -115747.0 / 13440.0,
     152             :           32863.0 / 13440.0, -5257.0 / 17280.0};
     153             :     default:
     154             :       ERROR("Bad order: " << order);
     155             :   }
     156             : }
     157             : 
     158             : void AdamsBashforthN::pup(PUP::er& p) noexcept {
     159             :   LtsTimeStepper::Inherit::pup(p);
     160             :   p | order_;
     161             : }
     162             : 
     163           0 : bool operator==(const AdamsBashforthN& lhs,
     164             :                 const AdamsBashforthN& rhs) noexcept {
     165             :   return lhs.order_ == rhs.order_;
     166             : }
     167             : 
     168           0 : bool operator!=(const AdamsBashforthN& lhs,
     169             :                 const AdamsBashforthN& rhs) noexcept {
     170             :   return not(lhs == rhs);
     171             : }
     172             : }  // namespace TimeSteppers
     173             : 
     174             : /// \cond
     175             : PUP::able::PUP_ID TimeSteppers::AdamsBashforthN::my_PUP_ID =  // NOLINT
     176             :     0;
     177             : /// \endcond

Generated by: LCOV version 1.14