ode-toolbox

RKF45.h (7164B)

---

 /* File:   RKF45.h */
 /* Date:   Wed Feb 10 17:33:42 2016 */
 /* Author: Fabian Wermelinger */
 /* Tag:    Runge-Kutta-Fehlberg 4th-5th order variable step method */
 /* Copyright 2016 ETH Zurich. All Rights Reserved. */
 #ifndef RKF45_H_UWN6WM1V
 #define RKF45_H_UWN6WM1V
 
 #include <ODETB/TimeStepper/explicit/Euler.h>
 
 #include <cassert>
 #include <cmath>
 #include <cstdlib>
 #include <limits>
 
 template <typename Tinput, typename Trhs>
 class RKF45 : public Euler<Tinput, Trhs>
 {
     using Euler<Tinput, Trhs>::m_settings;
     using Euler<Tinput, Trhs>::m_U;
     using Euler<Tinput, Trhs>::m_rhs;
     using Euler<Tinput, Trhs>::m_kernel;
 
     Real m_errold;
     const Real m_alpha;
     const Real m_beta;
     bool m_reject;
     bool m_firstPass;
 
     Trhs m_rhs2;
     Trhs m_rhs3;
     Trhs m_rhs4;
     Trhs m_rhs5;
     Trhs m_rhs6;
 
     Tinput m_output;
 
     static constexpr Real b1  = 16./135., b3 = 6656./12825., b4 = 28561./56430., b5 = -9./50., b6 = 2./55.;
     static constexpr Real c21 = 0.25, ct2 = 0.25;
     static constexpr Real c31 = 3./32., c32 = 9./32., ct3 = 3./8.;
     static constexpr Real c41 = 1932./2197., c42 = -7200./2197., c43 = 7296./2197., ct4 = 12./13.;
     static constexpr Real c51 = 439./216., c52 = -8., c53 = 3680./513., c54 = -845./4104., ct5 = 1.0;
     static constexpr Real c61 = -8./27., c62 = 2., c63 = -3544./2565., c64 = 1859./4104., c65 = -11./40., ct6 = 0.5;
     static constexpr Real e1 = 1./360., e3 = -128./4275., e4 = -2197./75240., e5 = 1./50., e6 = 2./55.;
 
     void _computeRHS(const Real t, const Real dt, void const* const data)
     {
         if (!m_reject) m_kernel->compute(m_U, m_rhs, t, data);
         m_kernel->compute(m_U + dt*(c21*m_rhs), m_rhs2, t + ct2*dt, data);
         m_kernel->compute(m_U + dt*(c31*m_rhs + c32*m_rhs2), m_rhs3, t + ct3*dt, data);
         m_kernel->compute(m_U + dt*(c41*m_rhs + c42*m_rhs2 + c43*m_rhs3), m_rhs4, t + ct4*dt, data);
         m_kernel->compute(m_U + dt*(c51*m_rhs + c52*m_rhs2 + c53*m_rhs3 + c54*m_rhs4), m_rhs5, t + ct5*dt, data);
         m_kernel->compute(m_U + dt*(c61*m_rhs + c62*m_rhs2 + c63*m_rhs3 + c64*m_rhs4 + c65*m_rhs5), m_rhs6, t + ct6*dt, data);
     }
 
     inline Real _error() const
     {
         assert(m_U.size() == m_rhs.size());
         Real maxerr = 0.0;
         // using inf-norm (could also do L2-norm)
         for (size_t i = 0; i < m_U.size(); ++i)
         {
             const typename Trhs::DataType E = e1*m_rhs[i] + e3*m_rhs3[i] + e4*m_rhs4[i] + e5*m_rhs5[i] + e6*m_rhs6[i];
             for (int j = 0; j < Trhs::DataType::SIZE; ++j)
             {
                 const Real IepsI = std::abs(E[j]);
                 const Real scale = m_settings.aTol + std::abs(m_U[i][j])*m_settings.rTol;
                 maxerr = std::max(maxerr, IepsI/scale);
             }
         }
         assert(!isnan(maxerr));
         return maxerr;
     }
 
     bool _inBound(const Real t, Real& dt, void const* const data, Real& dt_next)
     {
         const Real err = _error();
 
         const Real dtMin = m_settings.dtMin;
         if (dt <= dtMin) {
             dt_next = 1.1 * dtMin;
             m_errold = std::max(err, static_cast<Real>(1.0e-4));
             m_reject = false;
             m_firstPass = true;
             return true;
         }
 
         if (dt < std::numeric_limits<Real>::epsilon())
         {
             // if you request the impossible
             dt_next = 1.5*std::numeric_limits<Real>::epsilon();
             m_errold = std::max(err, static_cast<Real>(1.0e-4));
             m_reject = false;
             m_firstPass = true;
             return true;
         }
 
         const Real safety = m_settings.safety;
         const Real minScale = m_settings.minScale;
         const Real maxScale = m_settings.maxScale;
         Real scale;
         if (err <= 1.0)
         {
             if (err == 0.0)
                 scale = maxScale;
             else
             {
                 scale = safety*std::pow(err,-m_alpha)*std::pow(m_errold,m_beta);
                 scale = (scale < minScale) ? minScale : ((scale > maxScale) ? maxScale : scale);
             }
             if (m_reject)
                 dt_next = dt*std::min(scale, static_cast<Real>(1.0));
             else
                 dt_next = dt*scale;
             m_errold = std::max(err, static_cast<Real>(1.0e-4));
             m_reject = false;
             m_firstPass = true;
             return true;
         }
         else
         {
             scale = safety*std::pow(err,-m_alpha);
             scale = std::max(scale, minScale);
             dt *= scale;
             m_reject = true;
             m_firstPass = false;
             _computeRHS(t, dt, data);
             return false;
         }
     }
 
     void _dumpOutput(const void * const data)
     {
         const Real t  = m_settings.t;
         const Real dt = m_settings.dt;
         Real& tDump   = m_settings.tDump;
 
         const Tinput Uold = m_output;
         const Trhs rhsOld = m_rhs;
 
         m_kernel->compute(m_U, m_rhs, t+dt, data);
         m_firstPass = false;
 
         while(tDump <= (t+dt))
         {
             const Real phi = (tDump - t)/dt;
             // Hermite 3rd order
             m_output = (1.0-phi)*Uold + phi*m_U + phi*(phi-1.0)*((1.0-2.0*phi)*(m_U - Uold) + (phi-1.0)*dt*rhsOld + phi*dt*m_rhs);
             m_kernel->write(m_settings.step, t+phi*dt, phi*dt, m_output, data);
             tDump += m_settings.dtDump;
         }
     }
 
 public:
     RKF45(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) :
         Euler<Tinput,Trhs>(U,S,kern), m_errold(1.0e-4), m_alpha(0.25*(S.alpha-0.75*S.beta)), m_beta(0.25*S.beta), m_reject(false), m_firstPass(true),
         m_rhs2(U.size()), m_rhs3(U.size()), m_rhs4(U.size()), m_rhs5(U.size()), m_rhs6(U.size()), m_output(U.size()) {}
     virtual ~RKF45() {}
 
     virtual void step(void const* const data=nullptr)
     {
         size_t& step       = m_settings.step;
         Real& t            = m_settings.t;
         Real& dt           = m_settings.dt;
         const Real& tFinal = m_settings.tFinal;
         Real& tDump        = m_settings.tDump;
         if (!m_settings.bRestart && m_settings.dtDump > 0.0) tDump = t + m_settings.dtDump;
 
         while(t < tFinal)
         {
             Real dt_next;
             dt = (tFinal-t) < dt ? (tFinal-t) : dt;
             _computeRHS(t, dt, data);
             while (!_inBound(t, dt, data, dt_next)) {}
             m_output = m_U;
             m_U += dt*(b1*m_rhs + b3*m_rhs3 + b4*m_rhs4 + b5*m_rhs5 + b6*m_rhs6);
             if (m_settings.dtDump > 0.0)
             {
                 if (tDump <= (t+dt)) _dumpOutput(data);
             }
             else
                 if ((step+1) % m_settings.writeGranularity == 0)
                     m_kernel->write((step+1), t+dt, dt, m_U, data);
             t += dt;
             dt = dt_next;
             ++step;
 
             if (m_settings.step % m_settings.reportGranularity == 0)
                 std::printf(
                     "Time = %e;\tStep = %zu\n", m_settings.t, m_settings.step);
 
             if (m_settings.restartstep > 0 && step % m_settings.restartstep == 0)
                 serialize<Tinput>(m_U, m_settings);
         }
     }
 };
 
 #endif /* RKF45_H_UWN6WM1V */