RK4.h (1804B)
1 /* File: RK4.h */ 2 /* Date: Fri Feb 5 19:14:44 2016 */ 3 /* Author: Fabian Wermelinger */ 4 /* Tag: Classical 4th order Runge-Kutta */ 5 /* Copyright 2016 ETH Zurich. All Rights Reserved. */ 6 #ifndef RK4_H_LROUVHDE 7 #define RK4_H_LROUVHDE 8 9 #include <ODETB/TimeStepper/explicit/Euler.h> 10 11 template <typename Tinput, typename Trhs> 12 class RK4 : public Euler<Tinput, Trhs> 13 { 14 using Euler<Tinput, Trhs>::m_settings; 15 using Euler<Tinput, Trhs>::m_U; 16 using Euler<Tinput, Trhs>::m_rhs; 17 using Euler<Tinput, Trhs>::m_kernel; 18 19 public: 20 RK4(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) : Euler<Tinput,Trhs>(U,S,kern) { } 21 virtual ~RK4() {} 22 23 virtual void step(void const* const data=nullptr) 24 { 25 const Real oneHalf = static_cast<Real>(0.5); 26 const Real oneThird = static_cast<Real>(1./3.); 27 Tinput tmpU = m_U; 28 29 // stage 1 30 m_kernel->compute(tmpU, m_rhs, m_settings.t, data); 31 m_U += oneHalf*oneThird*m_settings.dt*m_rhs; 32 33 // stage 2 34 m_kernel->compute(tmpU + oneHalf*m_settings.dt*m_rhs, m_rhs, m_settings.t + oneHalf*m_settings.dt, data); 35 m_U += oneThird*m_settings.dt*m_rhs; 36 37 // stage 3 38 m_kernel->compute(tmpU + oneHalf*m_settings.dt*m_rhs, m_rhs, m_settings.t + oneHalf*m_settings.dt, data); 39 m_U += oneThird*m_settings.dt*m_rhs; 40 41 // stage 4 42 m_kernel->compute(tmpU + m_settings.dt*m_rhs, m_rhs, m_settings.t + m_settings.dt, data); 43 m_U += oneHalf*oneThird*m_settings.dt*m_rhs; 44 45 m_settings.t += m_settings.dt; 46 ++(m_settings.step); 47 48 if (m_settings.step % m_settings.reportGranularity == 0) 49 std::printf( 50 "Time = %e;\tStep = %zu\n", m_settings.t, m_settings.step); 51 } 52 }; 53 54 #endif /* RK4_H_LROUVHDE */