SpECTRE Documentation Coverage Report
Current view: top level - Time/TimeSteppers - RungeKutta3.hpp Hit Total Coverage
Commit: d07bc7e978f0161d6eb259d687cd7878b5f9f6f0 Lines: 2 23 8.7 %
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 RungeKutta3.
       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             : /// \ingroup TimeSteppersGroup
      35             : ///
      36             : /// A "strong stability-preserving" 3rd-order Runge-Kutta time-stepper.
      37             : /// Major reference:  J. Hesthaven & T. Warburton, Nodal Discontinuous
      38             : /// Galerkin Methods. section 5.7
      39           1 : class RungeKutta3 : public TimeStepper::Inherit {
      40             :  public:
      41           0 :   using options = tmpl::list<>;
      42           0 :   static constexpr Options::String help = {
      43             :       "A third-order strong stability-preserving Runge-Kutta time-stepper."};
      44             : 
      45           0 :   RungeKutta3() = default;
      46           0 :   RungeKutta3(const RungeKutta3&) noexcept = default;
      47           0 :   RungeKutta3& operator=(const RungeKutta3&) noexcept = default;
      48           0 :   RungeKutta3(RungeKutta3&&) noexcept = default;
      49           0 :   RungeKutta3& operator=(RungeKutta3&&) noexcept = default;
      50           0 :   ~RungeKutta3() noexcept override = default;
      51             : 
      52             :   template <typename Vars, typename DerivVars>
      53           0 :   void update_u(gsl::not_null<Vars*> u,
      54             :                 gsl::not_null<History<Vars, DerivVars>*> history,
      55             :                 const TimeDelta& time_step) const noexcept;
      56             : 
      57             :   template <typename Vars, typename DerivVars>
      58           0 :   void dense_update_u(gsl::not_null<Vars*> u,
      59             :                       const History<Vars, DerivVars>& history,
      60             :                       double time) const noexcept;
      61             : 
      62           0 :   uint64_t number_of_substeps() const noexcept override;
      63             : 
      64           0 :   size_t number_of_past_steps() const noexcept override;
      65             : 
      66           0 :   double stable_step() const noexcept override;
      67             : 
      68           0 :   TimeStepId next_time_id(const TimeStepId& current_id,
      69             :                           const TimeDelta& time_step) const noexcept override;
      70             : 
      71             :   template <typename Vars, typename DerivVars>
      72           0 :   bool can_change_step_size(
      73             :       const TimeStepId& time_id,
      74             :       const TimeSteppers::History<Vars, DerivVars>& /*history*/) const
      75             :       noexcept {
      76             :     return time_id.substep() == 0;
      77             :   }
      78             : 
      79           0 :   WRAPPED_PUPable_decl_template(RungeKutta3);  // NOLINT
      80             : 
      81           0 :   explicit RungeKutta3(CkMigrateMessage* /*unused*/) noexcept {}
      82             : 
      83             :   // clang-tidy: do not pass by non-const reference
      84           0 :   void pup(PUP::er& p) noexcept override {  // NOLINT
      85             :     TimeStepper::Inherit::pup(p);
      86             :   }
      87             : };
      88             : 
      89           0 : inline bool constexpr operator==(const RungeKutta3& /*lhs*/,
      90             :                                  const RungeKutta3& /*rhs*/) noexcept {
      91             :   return true;
      92             : }
      93             : 
      94           0 : inline bool constexpr operator!=(const RungeKutta3& /*lhs*/,
      95             :                                  const RungeKutta3& /*rhs*/) noexcept {
      96             :   return false;
      97             : }
      98             : 
      99             : template <typename Vars, typename DerivVars>
     100           0 : void RungeKutta3::update_u(
     101             :     const gsl::not_null<Vars*> u,
     102             :     const gsl::not_null<History<Vars, DerivVars>*> history,
     103             :     const TimeDelta& time_step) const noexcept {
     104             :   const size_t substep = history->size() - 1;
     105             :   const auto& vars = (history->end() - 1).value();
     106             :   const auto& dt_vars = (history->end() - 1).derivative();
     107             :   const auto& U0 = history->begin().value();
     108             : 
     109             :   switch (substep) {
     110             :     case 0: {
     111             :       // from (5.32) of Hesthaven
     112             :       // v^(1) = u^n + dt*RHS(u^n,t^n)
     113             :       // On entry V = u^n, U0 = u^n, rhs0 = RHS(u^n,t^n),
     114             :       // time = t^n
     115             :       *u += time_step.value() * dt_vars;
     116             :       // On exit v = v^(1), time = t^n + dt
     117             :       break;
     118             :     }
     119             :     case 1: {
     120             :       // from (5.32) of Hesthaven
     121             :       // v^(2) = (1/4)*( 3*u^n + v^(1) + dt*RHS(v^(1),t^n + dt) )
     122             :       // On entry V = v^(1), U0 = u^n, rhs0 = RHS(v^(1),t^n + dt),
     123             :       // time = t^n + dt
     124             :       *u += 0.25 * (3.0 * (U0 - vars) + time_step.value() * dt_vars);
     125             :       // On exit v = v^(2), time = t^n + (1/2)*dt
     126             :       break;
     127             :     }
     128             :     case 2: {
     129             :       // from (5.32) of Hesthaven
     130             :       // u^(n+1) = (1/3)*( u^n + 2*v^(2) + 2*dt*RHS(v^(2),t^n + (1/2)*dt) )
     131             :       // On entry V = v^(2), U0 = u^n, rhs0 = RHS(v^(2),t^n + (1/2)*dt),
     132             :       // time = t^n + (1/2)*dt
     133             :       *u += (1.0 / 3.0) * (U0 - vars + 2.0 * time_step.value() * dt_vars);
     134             :       // On exit v = u^(n+1), time = t^n + dt
     135             :       break;
     136             :     }
     137             :     default:
     138             :       ERROR("Bad substep value in RK3: " << substep);
     139             :   }
     140             : 
     141             :   // Clean up old history
     142             :   if (history->size() == number_of_substeps()) {
     143             :     history->mark_unneeded(history->end());
     144             :   }
     145             : }
     146             : 
     147             : template <typename Vars, typename DerivVars>
     148             : void RungeKutta3::dense_update_u(gsl::not_null<Vars*> u,
     149             :                                  const History<Vars, DerivVars>& history,
     150             :                                  const double time) const noexcept {
     151             :   ASSERT(history.size() == 3, "Can only dense output on last substep");
     152             :   const double step_start = history[0].value();
     153             :   const double step_end = history[1].value();
     154             :   const double time_step = step_end - step_start;
     155             :   const double output_fraction = (time - step_start) / time_step;
     156             :   ASSERT(output_fraction >= 0, "Attempting dense output at time " << time
     157             :          << ", but already progressed past " << step_start);
     158             :   ASSERT(output_fraction <= 1,
     159             :          "Requested time (" << time << " not within step [" << step_start
     160             :          << ", " << step_end << "]");
     161             : 
     162             :   // arXiv:1605.02429
     163             :   *u += (1.0 - output_fraction * (1.0 - output_fraction / 3.0)) *
     164             :             history.begin().value() +
     165             :         output_fraction * (1.0 - output_fraction) *
     166             :             (history.begin() + 1).value() +
     167             :         (2.0 / 3.0 * square(output_fraction) - 1.0) *
     168             :             (history.begin() + 2).value() +
     169             :         2.0 / 3.0 * square(output_fraction) * time_step *
     170             :             (history.begin() + 2).derivative();
     171             : }
     172             : }  // namespace TimeSteppers

Generated by: LCOV version 1.14