SpECTRE Documentation Coverage Report
Current view: top level - Time/TimeSteppers - RungeKutta4.hpp Hit Total Coverage
Commit: d07bc7e978f0161d6eb259d687cd7878b5f9f6f0 Lines: 2 24 8.3 %
Date: 2020-09-24 02:51:00
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 RungeKutta4.
       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"  // IWYU pragma: keep
      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 4th-order Runge-Kutta method, given e.g. in
      38             :  * https://en.wikipedia.org/wiki/Runge-Kutta_methods
      39             :  * that solves equations of the form
      40             :  *
      41             :  * \f{eqnarray}{
      42             :  * \frac{du}{dt} & = & \mathcal{L}(t,u).
      43             :  * \f}
      44             :  * Given a solution \f$u(t^n)=u^n\f$, this stepper computes
      45             :  * \f$u(t^{n+1})=u^{n+1}\f$ using the following equations:
      46             :  *
      47             :  * \f{eqnarray}{
      48             :  * v^{(1)} & = & u^n + dt\cdot \mathcal{L}(t^n, u^n)/2,\\
      49             :  * v^{(2)} & = & u^n + dt\cdot \mathcal{L}(t^n + dt/2, v^{(1)})/2,\\
      50             :  * v^{(3)} & = & u^n + dt\cdot \mathcal{L}(t^n + dt/2, v^{(2)}),\\
      51             :  * v^{(4)} & = & u^n + dt\cdot \mathcal{L}(t^n + dt, v^{(3)}),\\
      52             :  * u^{n+1} & = & (2v^{(1)} + 4v^{(2)} + 2v^{(3)} + v^{(4)} - 3 u^n)/6.
      53             :  * \f}
      54             :  *
      55             :  * Note that in the implementation, the expression for \f$u^{n+1}\f$ is
      56             :  * computed simultaneously with \f$v^{(4)}\f$, so that there are
      57             :  * actually only four substeps per step.
      58             :  */
      59           1 : class RungeKutta4 : public TimeStepper::Inherit {
      60             :  public:
      61           0 :   using options = tmpl::list<>;
      62           0 :   static constexpr Options::String help = {
      63             :       "The standard fourth-order Runge-Kutta time-stepper."};
      64             : 
      65           0 :   RungeKutta4() = default;
      66           0 :   RungeKutta4(const RungeKutta4&) noexcept = default;
      67           0 :   RungeKutta4& operator=(const RungeKutta4&) noexcept = default;
      68           0 :   RungeKutta4(RungeKutta4&&) noexcept = default;
      69           0 :   RungeKutta4& operator=(RungeKutta4&&) noexcept = default;
      70           0 :   ~RungeKutta4() noexcept override = default;
      71             : 
      72             :   template <typename Vars, typename DerivVars>
      73           0 :   void update_u(gsl::not_null<Vars*> u,
      74             :                 gsl::not_null<History<Vars, DerivVars>*> history,
      75             :                 const TimeDelta& time_step) const noexcept;
      76             : 
      77             :   template <typename Vars, typename DerivVars>
      78           0 :   void dense_update_u(gsl::not_null<Vars*> u,
      79             :                       const History<Vars, DerivVars>& history,
      80             :                       double time) const noexcept;
      81             : 
      82           0 :   uint64_t number_of_substeps() const noexcept override;
      83             : 
      84           0 :   size_t number_of_past_steps() const noexcept override;
      85             : 
      86           0 :   double stable_step() const noexcept override;
      87             : 
      88           0 :   TimeStepId next_time_id(const TimeStepId& current_id,
      89             :                           const TimeDelta& time_step) const noexcept override;
      90             : 
      91             :   template <typename Vars, typename DerivVars>
      92           0 :   bool can_change_step_size(
      93             :       const TimeStepId& time_id,
      94             :       const TimeSteppers::History<Vars, DerivVars>& /*history*/) const
      95             :       noexcept {
      96             :     return time_id.substep() == 0;
      97             :   }
      98             : 
      99           0 :   WRAPPED_PUPable_decl_template(RungeKutta4);  // NOLINT
     100             : 
     101           0 :   explicit RungeKutta4(CkMigrateMessage* /*unused*/) noexcept {}
     102             : 
     103             :   // clang-tidy: do not pass by non-const reference
     104           0 :   void pup(PUP::er& p) noexcept override {  // NOLINT
     105             :     TimeStepper::Inherit::pup(p);
     106             :   }
     107             : };
     108             : 
     109           0 : inline bool constexpr operator==(const RungeKutta4& /*lhs*/,
     110             :                                  const RungeKutta4& /*rhs*/) noexcept {
     111             :   return true;
     112             : }
     113             : 
     114           0 : inline bool constexpr operator!=(const RungeKutta4& /*lhs*/,
     115             :                                  const RungeKutta4& /*rhs*/) noexcept {
     116             :   return false;
     117             : }
     118             : 
     119             : template <typename Vars, typename DerivVars>
     120           0 : void RungeKutta4::update_u(
     121             :     const gsl::not_null<Vars*> u,
     122             :     const gsl::not_null<History<Vars, DerivVars>*> history,
     123             :     const TimeDelta& time_step) const noexcept {
     124             :   const size_t substep = history->size() - 1;
     125             :   const auto& dt_vars = (history->end() - 1).derivative();
     126             :   const auto& u0 = history->begin().value();
     127             : 
     128             :   switch (substep) {
     129             :     case 0: {
     130             :       // from (17.1.3) of Numerical Recipes 3rd Edition
     131             :       // v^(1) = u^n + dt * \mathcal{L}(u^n,t^n)/2
     132             :       // On entry V = u^n, u0 = u^n, rhs0 = \mathcal{L}(u^n, t^n),
     133             :       // time = t^n
     134             :       *u += 0.5 * time_step.value() * dt_vars;
     135             :       // On exit v = v^(1), time = t^n + (1/2)*dt
     136             :       break;
     137             :     }
     138             :     case 1: {
     139             :       // from (17.1.3) of Numerical Recipes 3rd Edition
     140             :       // v^(2) = u^n + dt * \mathcal{L}(v^(1), t^n + (1/2)*dt)/2
     141             :       // On entry V = v^(1), u0 = u^n, rhs0 = \mathcal{L}(v^(1), t^n + dt/2),
     142             :       // time = t^n + dt
     143             :       *u = u0 + 0.5 * time_step.value() * dt_vars;
     144             :       // On exit v = v^(2), time = t^n + (1/2)*dt
     145             :       break;
     146             :     }
     147             :     case 2: {
     148             :       // from (17.1.3) of Numerical Recipes 3rd Edition
     149             :       // v^(3) = u^n + dt * \mathcal{L}(v^(2), t^n + (1/2)*dt))
     150             :       // On entry V = v^(2), u0 = u^n,
     151             :       // rhs0 = \mathcal{L}(v^(2), t^n + (1/2)*dt), time = t^n + (1/2)*dt
     152             :       *u = u0 + time_step.value() * dt_vars;
     153             :       // On exit v = v^(3), time = t^n + dt
     154             :       break;
     155             :     }
     156             :     case 3: {
     157             :       // from (17.1.3) of Numerical Recipes 3rd Edition
     158             :       // u^(n+1) = (2v^(1) + 4*v^(2) + 2*v^(3) + v^(4) - 3*u0)/6
     159             :       // On entry V = v^(3), u0 = u^n, rhs0 = \mathcal{L}(v^(3), t^n + dt),
     160             :       // time = t^n + dt
     161             :       // Note: v^(4) = u0 + dt * \mathcal{L}(t+dt, v^(3)); inserting this gives
     162             :       // u^(n+1) = (2v^(1) + 4*v^(2) + 2*v^(3)
     163             :       //         + dt*\mathcal{L}(t+dt,v^(3)) - 2*u0)/6
     164             :       constexpr double one_sixth = 1.0 / 6.0;
     165             :       *u = (2.0 * (history->begin() + 1).value() +
     166             :             4.0 * (history->begin() + 2).value() +
     167             :             2.0 * (history->begin() + 3).value() +
     168             :             (time_step.value() * dt_vars - 2.0 * u0)) *
     169             :            one_sixth;
     170             :       // On exit v = u^(n+1), time = t^n + dt
     171             :       break;
     172             :     }
     173             :     default:
     174             :       ERROR("Substep in RK4 should be one of 0,1,2,3, 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 RungeKutta4::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             :          "RK4 can only dense output on last substep ("
     189             :              << number_of_substeps() - 1 << "), not substep "
     190             :              << history.size() - 1);
     191             :   const double step_start = history[0].value();
     192             :   const double step_end = history[history.size() - 1].value();
     193             :   const double time_step = step_end - step_start;
     194             :   const double output_fraction = (time - step_start) / time_step;
     195             :   ASSERT(output_fraction >= 0.0, "Attempting dense output at time "
     196             :                                      << time << ", but already progressed past "
     197             :                                      << step_start);
     198             :   ASSERT(output_fraction <= 1.0, "Requested time ("
     199             :                                      << time << ") not within step ["
     200             :                                      << step_start << ", " << step_end << "]");
     201             : 
     202             :   // Numerical Recipes Eq. (17.2.15). This implements cubic interpolation
     203             :   // throughout the step. Because the history only is available through the
     204             :   // penultimate step, i) the value after the step is computed algebraically
     205             :   // from previous substeps, and ii) the derivative at the final step is
     206             :   // approximated as the derivative at the penultimate substep.
     207             :   constexpr double one_sixth = 1.0 / 6.0;
     208             :   const auto& u0 = history.begin().value();
     209             :   const auto& dt_vars = (history.end() - 1).derivative();
     210             :   const Vars u1 =
     211             :       (2.0 * (history.begin() + 1).value() +
     212             :        4.0 * (history.begin() + 2).value() +
     213             :        2.0 * (history.begin() + 3).value() + (time_step * dt_vars - 2.0 * u0)) *
     214             :       one_sixth;
     215             :   *u = (1.0 - output_fraction) * u0 + output_fraction * u1 +
     216             :        (output_fraction) * (output_fraction - 1.0) *
     217             :            ((1.0 - 2.0 * output_fraction) * (u1 - u0) +
     218             :             (output_fraction - 1.0) * time_step * history.begin().derivative() +
     219             :             output_fraction * time_step * dt_vars);
     220             : }
     221             : }  // namespace TimeSteppers

Generated by: LCOV version 1.14