ode-toolbox

ODE integration tools
git clone https://git.0xfab.ch/ode-toolbox.git
Log | Files | Refs | README | LICENSE

LSRK3.h (1821B)


      1 /* File:   LSRK3.h */
      2 /* Date:   Fri Feb  5 18:58:48 2016 */
      3 /* Author: Fabian Wermelinger */
      4 /* Tag:    3rd order low-storage Runge-Kutta.  Requires a compatible kernel */
      5 /* Copyright 2016 ETH Zurich. All Rights Reserved. */
      6 #ifndef LSRK3_H_MJ6ACLH4
      7 #define LSRK3_H_MJ6ACLH4
      8 
      9 #include <ODETB/TimeStepper/explicit/Euler.h>
     10 
     11 template <typename Tinput, typename Trhs>
     12 class LSRK3 : public Euler<Tinput, Trhs>
     13 {
     14     const Real m_A[3];
     15     const Real m_B[3];
     16     using Euler<Tinput,Trhs>::m_settings;
     17     using Euler<Tinput,Trhs>::m_U;
     18     using Euler<Tinput,Trhs>::m_rhs;
     19     using Euler<Tinput,Trhs>::m_kernel;
     20 
     21 public:
     22     LSRK3(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) :
     23         Euler<Tinput,Trhs>(U,S,kern),
     24         // m_A{0.0, -2.915492524638791, -0.000000093517376}, // Gottlieb & Shu (Optimal LSRK3 for CFL = 0.32)
     25         // m_B{0.924574000000000, 0.287713063186749, 0.626538109512740}
     26         m_A{ 0.0, -17./32, -32./27},
     27         m_B{1./4,    8./9,    3./4}
     28     { }
     29     virtual ~LSRK3() {}
     30 
     31     virtual void step(void const *const)
     32     {
     33         // stage 1
     34         Real lsrk_data[2] = {static_cast<Real>(m_settings.dt), m_A[0]};
     35         m_kernel->compute(m_U, m_rhs, m_settings.t, lsrk_data);
     36         m_U += m_B[0]*m_rhs;
     37 
     38         // stage 2
     39         lsrk_data[1] = m_A[1];
     40         m_kernel->compute(m_U, m_rhs, m_settings.t, lsrk_data);
     41         m_U += m_B[1]*m_rhs;
     42 
     43         // stage 3
     44         lsrk_data[1] = m_A[2];
     45         m_kernel->compute(m_U, m_rhs, m_settings.t, lsrk_data);
     46         m_U += m_B[2]*m_rhs;
     47 
     48         m_settings.t += m_settings.dt;
     49         ++(m_settings.step);
     50 
     51         if (m_settings.step % m_settings.reportGranularity == 0)
     52             std::printf(
     53                 "Time = %e;\tStep = %zu\n", m_settings.t, m_settings.step);
     54     }
     55 };
     56 
     57 #endif /* LSRK3_H_MJ6ACLH4 */