SpECTRE Documentation Coverage Report
Current view: top level - Time/TimeSteppers - DormandPrince5.hpp Hit Total Coverage
Commit: b1342d46f40e2d46bbd11d0cef68fd973031a24b Lines: 2 32 6.2 %
Date: 2020-09-24 20:24:42
Legend: Lines: hit not hit

          Line data    Source code
       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines class DormandPrince5.
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <cstddef>
      10             : #include <cstdint>
      11             : #include <ostream>
      12             : #include <pup.h>
      13             : 
      14             : #include "ErrorHandling/Assert.hpp"
      15             : #include "ErrorHandling/Error.hpp"
      16             : #include "Options/Options.hpp"
      17             : #include "Parallel/CharmPupable.hpp"
      18             : #include "Time/Time.hpp"
      19             : #include "Time/TimeSteppers/TimeStepper.hpp"
      20             : #include "Utilities/ConstantExpressions.hpp"
      21             : #include "Utilities/Gsl.hpp"
      22             : #include "Utilities/TMPL.hpp"
      23             : 
      24             : /// \cond
      25             : struct TimeStepId;
      26             : namespace TimeSteppers {
      27             : template <typename Vars, typename DerivVars>
      28             : class History;
      29             : }  // namespace TimeSteppers
      30             : /// \endcond
      31             : 
      32             : namespace TimeSteppers {
      33             : 
      34             : /*!
      35             :  * \ingroup TimeSteppersGroup
      36             :  *
      37             :  * The standard 5th-order Dormand-Prince time stepping method, given e.g. in
      38             :  * Sec. 7.2 of \cite NumericalRecipes.
      39             :  *
      40             :  * \f{eqnarray}{
      41             :  * \frac{du}{dt} & = & \mathcal{L}(t,u).
      42             :  * \f}
      43             :  * Given a solution \f$u(t^n)=u^n\f$, this stepper computes
      44             :  * \f$u(t^{n+1})=u^{n+1}\f$ using the following equations:
      45             :  *
      46             :  * \f{align}{
      47             :  * k^{(1)} & = dt \mathcal{L}(t^n, u^n),\\
      48             :  * k^{(i)} & = dt \mathcal{L}(t^n + c_i dt,
      49             :  *                              u^n + \sum_{j=1}^{i-1} a_{ij} k^{(j)}),
      50             :  *                              \mbox{ } 2 \leq i \leq 6,\\
      51             :  * u^{n+1} & = u^n + \sum_{i=1}^{6} b_i k^{(i)}.
      52             :  * \f}
      53             :  *
      54             :  * Here the coefficients \f$a_{ij}\f$, \f$b_i\f$, and \f$c_i\f$ are given
      55             :  * in e.g. Sec. 7.2 of \cite NumericalRecipes. Note that \f$c_1 = 0\f$.
      56             :  *
      57             :  */
      58           1 : class DormandPrince5 : public TimeStepper::Inherit {
      59             :  public:
      60           0 :   using options = tmpl::list<>;
      61           0 :   static constexpr Options::String help = {
      62             :       "The standard Dormand-Prince 5th-order time stepper."};
      63             : 
      64           0 :   DormandPrince5() = default;
      65           0 :   DormandPrince5(const DormandPrince5&) noexcept = default;
      66           0 :   DormandPrince5& operator=(const DormandPrince5&) noexcept = default;
      67           0 :   DormandPrince5(DormandPrince5&&) noexcept = default;
      68           0 :   DormandPrince5& operator=(DormandPrince5&&) noexcept = default;
      69           0 :   ~DormandPrince5() noexcept override = default;
      70             : 
      71             :   template <typename Vars, typename DerivVars>
      72           0 :   void update_u(gsl::not_null<Vars*> u,
      73             :                 gsl::not_null<History<Vars, DerivVars>*> history,
      74             :                 const TimeDelta& time_step) const noexcept;
      75             : 
      76             :   template <typename Vars, typename DerivVars>
      77           0 :   void dense_update_u(gsl::not_null<Vars*> u,
      78             :                       const History<Vars, DerivVars>& history,
      79             :                       double time) const noexcept;
      80             : 
      81           0 :   uint64_t number_of_substeps() const noexcept override;
      82             : 
      83           0 :   size_t number_of_past_steps() const noexcept override;
      84             : 
      85           0 :   double stable_step() const noexcept override;
      86             : 
      87           0 :   TimeStepId next_time_id(const TimeStepId& current_id,
      88             :                           const TimeDelta& time_step) const noexcept override;
      89             : 
      90             :   template <typename Vars, typename DerivVars>
      91           0 :   bool can_change_step_size(
      92             :       const TimeStepId& time_id,
      93             :       const TimeSteppers::History<Vars, DerivVars>& /*history*/) const
      94             :       noexcept {
      95             :     return time_id.substep() == 0;
      96             :   }
      97             : 
      98           0 :   WRAPPED_PUPable_decl_template(DormandPrince5);  // NOLINT
      99             : 
     100           0 :   explicit DormandPrince5(CkMigrateMessage* /*unused*/) noexcept {}
     101             : 
     102             :   // clang-tidy: do not pass by non-const reference
     103           0 :   void pup(PUP::er& p) noexcept override {  // NOLINT
     104             :     TimeStepper::Inherit::pup(p);
     105             :   }
     106             : 
     107             :  private:
     108             :   // Coefficients from the Dormand-Prince 5 Butcher tableau (e.g. Sec. 7.2
     109             :   // of \cite NumericalRecipes).
     110           0 :   static constexpr double _a2 = 0.2;
     111           0 :   static constexpr std::array<double, 2> _a3{{3.0 / 40.0, 9.0 / 40.0}};
     112           0 :   static constexpr std::array<double, 3> _a4{
     113             :       {44.0 / 45.0, -56.0 / 15.0, 32.0 / 9.0}};
     114           0 :   static constexpr std::array<double, 4> _a5{
     115             :       {19372.0 / 6561.0, -25360.0 / 2187.0, 64448.0 / 6561.0, -212.0 / 729.0}};
     116           0 :   static constexpr std::array<double, 5> _a6{{9017.0 / 3168.0, -355.0 / 33.0,
     117             :                                               46732.0 / 5247.0, 49.0 / 176.0,
     118             :                                               -5103.0 / 18656.0}};
     119           0 :   static constexpr std::array<double, 6> _b{{35.0 / 384.0, 0.0, 500.0 / 1113.0,
     120             :                                              125.0 / 192.0, -2187.0 / 6784.0,
     121             :                                              11.0 / 84.0}};
     122           0 :   static const std::array<Time::rational_t, 5> _c;
     123             : 
     124             :   // Coefficients for dense output, taken from Sec. 7.2 of
     125             :   // \cite NumericalRecipes
     126           0 :   static constexpr std::array<double, 6> _d{
     127             :       {-12715105075.0 / 11282082432.0, 87487479700.0 / 32700410799.0,
     128             :        -10690763975.0 / 1880347072.0, 701980252875.0 / 199316789632.0,
     129             :        -1453857185.0 / 822651844.0, 69997945.0 / 29380423.0}};
     130             : };
     131             : 
     132           0 : inline bool constexpr operator==(const DormandPrince5& /*lhs*/,
     133             :                                  const DormandPrince5& /*rhs*/) noexcept {
     134             :   return true;
     135             : }
     136             : 
     137           0 : inline bool constexpr operator!=(const DormandPrince5& /*lhs*/,
     138             :                                  const DormandPrince5& /*rhs*/) noexcept {
     139             :   return false;
     140             : }
     141             : 
     142             : template <typename Vars, typename DerivVars>
     143           0 : void DormandPrince5::update_u(
     144             :     const gsl::not_null<Vars*> u,
     145             :     const gsl::not_null<History<Vars, DerivVars>*> history,
     146             :     const TimeDelta& time_step) const noexcept {
     147             :   const size_t substep = history->size() - 1;
     148             :   const auto& u0 = history->begin().value();
     149             :   const double dt = time_step.value();
     150             : 
     151             :   const auto increment_u = [&u, &history, &dt](const auto& coeffs) noexcept {
     152             :     for (size_t i = 0; i < coeffs.size(); ++i) {
     153             :       *u += (gsl::at(coeffs, i) * dt) *
     154             :             (history->begin() + static_cast<int>(i)).derivative();
     155             :     }
     156             :   };
     157             : 
     158             :   if (substep == 0) {
     159             :     *u += (_a2 * dt) * history->begin().derivative();
     160             :   } else if (substep < 6) {
     161             :     *u = u0;
     162             :     if (substep == 1) {
     163             :       increment_u(_a3);
     164             :     } else if (substep == 2) {
     165             :       increment_u(_a4);
     166             :     } else if (substep == 3) {
     167             :       increment_u(_a5);
     168             :     } else if (substep == 4) {
     169             :       increment_u(_a6);
     170             :     } else {
     171             :       increment_u(_b);
     172             :     }
     173             :   } else {
     174             :     ERROR("Substep in DP5 should be one of 0,1,2,3,4,5, not " << substep);
     175             :   }
     176             : 
     177             :   // Clean up old history
     178             :   if (history->size() == number_of_substeps()) {
     179             :     history->mark_unneeded(history->end());
     180             :   }
     181             : }
     182             : 
     183             : template <typename Vars, typename DerivVars>
     184           0 : void DormandPrince5::dense_update_u(const gsl::not_null<Vars*> u,
     185             :                                     const History<Vars, DerivVars>& history,
     186             :                                     const double time) const noexcept {
     187             :   ASSERT(history.size() == number_of_substeps(),
     188             :          "DP5 can only dense output on last substep ("
     189             :              << number_of_substeps() - 1 << "), not substep "
     190             :              << history.size() - 1);
     191             :   const double t0 = history[0].value();
     192             :   const double t_end = history[history.size() - 1].value();
     193             :   // The history does not contain the final step; specifically,
     194             :   // step_end = t0 + c[4] * dt, so dt = (t_end - t0) / c[4]. But since
     195             :   // c[4] = 1.0 for DP5, we don't need to divide by c[4] here.
     196             :   const double dt = t_end - t0;
     197             :   const double output_fraction = (time - t0) / dt;
     198             :   ASSERT(output_fraction >= 0.0, "Attempting dense output at time "
     199             :                                      << time << ", but already progressed past "
     200             :                                      << t0);
     201             :   ASSERT(output_fraction <= 1.0, "Requested time ("
     202             :                                      << time << ") not within step [" << t0
     203             :                                      << ", " << t0 + dt << "]");
     204             : 
     205             :   // Get the solution u1 at the final time
     206             :   const auto& u0 = history.begin().value();
     207             :   Vars u1 = u0;
     208             :   for (size_t i = 0; i < 6; ++i) {
     209             :     u1 += (gsl::at(_b, i) * dt) *
     210             :           (history.begin() + static_cast<int>(i)).derivative();
     211             :   }
     212             : 
     213             :   // We need the following: k1, k3, k4, k5, k6.
     214             :   // Here k1 = dt * l1, k3 = dt * l3, etc.
     215             :   const auto& l1 = history.begin().derivative();
     216             :   const auto& l3 = (history.begin() + 2).derivative();
     217             :   const auto& l4 = (history.begin() + 3).derivative();
     218             :   const auto& l5 = (history.begin() + 4).derivative();
     219             :   const auto& l6 = (history.begin() + 5).derivative();
     220             : 
     221             :   // Compute the updating coefficents, called rcontN in Numerical recipes,
     222             :   // that will be reused, so I don't have to compute them more than once.
     223             :   const Vars rcont2 = u1 - u0;
     224             :   const Vars rcont3 = dt * l1 - rcont2;
     225             : 
     226             :   // The formula for dense output is given in Numerical Recipes Sec. 17.2.3.
     227             :   // Note: L(t+dt, u1) after the step is unavailable in the history; so here,
     228             :   // approximate L(t+dt, u1) by l6.
     229             :   *u = u0 + output_fraction *
     230             :                 (rcont2 +
     231             :                  (1.0 - output_fraction) *
     232             :                      (rcont3 + output_fraction *
     233             :                                    ((rcont2 - dt * l6 - rcont3) +
     234             :                                     ((1.0 - output_fraction) * dt) *
     235             :                                         (_d[0] * l1 + _d[1] * l3 + _d[2] * l4 +
     236             :                                          _d[3] * l5 + (_d[4] + _d[5]) * l6))));
     237             : }
     238             : }  // namespace TimeSteppers

Generated by: LCOV version 1.14