ode-toolbox

RKV56.h (7870B)

---

 /* File:   RKV56.h */
 /* Date:   Thu Feb 11 11:36:11 2016 */
 /* Author: Fabian Wermelinger */
 /* Tag:    Runge-Kutta-Verner 5th-6th order variable step method */
 /* Copyright 2016 ETH Zurich. All Rights Reserved. */
 #ifndef RKV56_H_A6VN39JN
 #define RKV56_H_A6VN39JN
 
 #include <ODETB/TimeStepper/explicit/Euler.h>
 
 #include <cassert>
 #include <cmath>
 #include <cstdlib>
 #include <iostream>
 #include <limits>
 
 template <typename Tinput, typename Trhs>
 class RKV56 : 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;
     Trhs m_rhs7;
     Trhs m_rhs8;
 
     Tinput m_output;
 
     static constexpr Real b1  = 3./40., b3 = 875./2244., b4 = 23./72., b5 = 264./1955., b7 = 125./11592., b8 = 43./616.;
     static constexpr Real c21 = 1./6., ct2 = 1./6.;
     static constexpr Real c31 = 4./75., c32 = 16./75., ct3 = 4./15.;
     static constexpr Real c41 = 5./6., c42 = -8./3., c43 = 2.5, ct4 = 2./3.;
     static constexpr Real c51 = -165./64., c52 = 55./6., c53 = -425./64., c54 = 85./96., ct5 = 5./6.;
     static constexpr Real c61 = 2.4, c62 = -8., c63 = 4015./612., c64 = -11./36., c65 = 88./255., ct6 = 1.0;
     static constexpr Real c71 = -8263./15000., c72 = 124./75., c73 = -643./680., c74 = -81./250., c75 = 2484./10625., ct7 = 1./15.;
     static constexpr Real c81 = 3501./1720., c82 = -300./43., c83 = 297275./52632., c84 = -319./2322., c85 = 24068./84065., c87 = 3850./26703., ct8 = 1.0;
     static constexpr Real e1  = -1./160., e3 = -125./17952., e4 = 1./144., e5 = -12./1955., e6 = -3./44., e7 = 125./11592., e8 = 43./616.;
 
     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);
         m_kernel->compute(m_U + dt*(c71*m_rhs + c72*m_rhs2 + c73*m_rhs3 + c74*m_rhs4 + c75*m_rhs5), m_rhs7, t + ct7*dt, data);
         m_kernel->compute(m_U + dt*(c81*m_rhs + c82*m_rhs2 + c83*m_rhs3 + c84*m_rhs4 + c85*m_rhs5 + c87*m_rhs7), m_rhs8, t + ct8*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] + e7*m_rhs7[i] + e8*m_rhs8[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:
     RKV56(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_rhs7(U.size()), m_rhs8(U.size()), m_output(U.size()) {}
     virtual ~RKV56() {}
 
     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 + b7*m_rhs7 + b8*m_rhs8);
             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 /* RKV56_H_A6VN39JN */