ceph/src/boost/boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp

   1 /*
   2  boost/numeric/odeint/stepper/detail/pid_step_adjuster.hpp
   3
   4  [begin_description]
   5  Implementation of the stepsize controller for the controlled adams bashforth moulton stepper.
   6  [end_description]
   7
   8  Copyright 2017 Valentin Noah Hartmann
   9
  10  Distributed under the Boost Software License, Version 1.0.
  11  (See accompanying file LICENSE_1_0.txt or
  12  copy at http://www.boost.org/LICENSE_1_0.txt)
  13  */
  14
  15 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_HPP_INCLUDED
  16 #define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_PID_STEP_ADJUSTER_HPP_INCLUDED
  17
  18 #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
  19 #include <boost/numeric/odeint/stepper/detail/pid_step_adjuster_coefficients.hpp>
  20
  21 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
  22 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
  23
  24 #include <math.h>
  25
  26 namespace boost {
  27 namespace numeric {
  28 namespace odeint {
  29 namespace detail {
  30
  31 template<
  32 class Value = double,
  33 class Time = double
  34 >
  35 struct pid_op
  36 {
  37 public:
  38     typedef Value value_type;
  39     typedef Time time_type;
  40
  41     const double beta1;
  42     const double beta2;
  43     const double beta3;
  44     const double alpha1;
  45     const double alpha2;
  46
  47     const time_type dt1;
  48     const time_type dt2;
  49     const time_type dt3;
  50
  51     const size_t m_steps;
  52
  53     pid_op(const size_t steps, const double _dt1, const double _dt2, const double _dt3,
  54         const double b1 = 1, const double b2 = 0, const double b3 = 0, const double a1 = 0, const double a2 = 0)
  55     :beta1(b1), beta2(b2), beta3(b3), alpha1(a1), alpha2(a2),
  56     dt1(_dt1), dt2(_dt2), dt3(_dt3),
  57     m_steps(steps)
  58     {};
  59
  60     template<class T1, class T2, class T3, class T4>
  61     void operator()(T1 &t1, const T2 &t2, const T3 &t3, const T4 &t4)
  62     {
  63         using std::abs;
  64
  65         t1 = adapted_pow(abs(t2), -beta1/(m_steps + 1)) *
  66             adapted_pow(abs(t3), -beta2/(m_steps + 1)) *
  67             adapted_pow(abs(t4), -beta3/(m_steps + 1)) *
  68             adapted_pow(abs(dt1/dt2), -alpha1/(m_steps + 1))*
  69             adapted_pow(abs(dt2/dt3), -alpha2/(m_steps + 1));
  70
  71         t1 = 1/t1;
  72     };
  73
  74     template<class T1, class T2>
  75     void operator()(T1 &t1, const T2 &t2)
  76     {
  77         using std::abs;
  78
  79         t1 = adapted_pow(abs(t2), -beta1/(m_steps + 1));
  80
  81         t1 = 1/t1;
  82     };
  83
  84 private:
  85     template<class T>
  86     inline value_type adapted_pow(T base, double exp)
  87     {
  88         if(exp == 0)
  89         {
  90             return 1;
  91         }
  92         else if (exp > 0)
  93         {
  94             return pow(base, exp);
  95         }
  96         else
  97         {
  98             return 1/pow(base, -exp);
  99         }
 100     };
 101 };
 102
 103 template<
 104 class State,
 105 class Value = double,
 106 class Deriv = State,
 107 class Time = double,
 108 class Algebra = typename algebra_dispatcher< State >::algebra_type,
 109 class Operations = typename operations_dispatcher< Deriv >::operations_type,
 110 size_t Type = BASIC
 111 >
 112 struct pid_step_adjuster
 113 {
 114 public:
 115     static double threshold() { return 0.9; };
 116
 117     typedef State state_type;
 118     typedef Value value_type;
 119     typedef Deriv deriv_type;
 120     typedef Time time_type;
 121
 122     typedef Algebra algebra_type;
 123     typedef Operations operations_type;
 124
 125     typedef rotating_buffer<state_type, 3> error_storage_type;
 126     typedef rotating_buffer<time_type, 3> time_storage_type;
 127     typedef pid_step_adjuster_coefficients<Type> coeff_type;
 128
 129     pid_step_adjuster(double abs_tol = 1e-6, double rel_tol = 1e-6, time_type dtmax = 1.0)
 130     :m_dtmax(dtmax), m_error_storage(), m_dt_storage(), m_init(0),
 131     m_abs_tol(abs_tol), m_rel_tol(rel_tol)
 132     {};
 133
 134     time_type adjust_stepsize(const size_t steps, time_type dt, state_type &err, const state_type &x, const deriv_type &dxdt)
 135     {
 136         using std::abs;
 137         m_algebra.for_each3( err , x , dxdt ,
 138                 typename operations_type::template rel_error< value_type >( m_abs_tol , m_rel_tol , 1.0 , 1.0 * abs(get_unit_value( dt )) ) );
 139
 140         m_error_storage[0] = err;
 141         m_dt_storage[0] = dt;
 142
 143         if(m_init >= 2)
 144         {
 145             m_algebra.for_each4(err, m_error_storage[0], m_error_storage[1], m_error_storage[2],
 146                 pid_op<>(steps, m_dt_storage[0], m_dt_storage[1], m_dt_storage[2],
 147                     m_coeff[0], m_coeff[1], m_coeff[2], m_coeff[3], m_coeff[4]));
 148         }
 149         else
 150         {
 151             m_algebra.for_each2(err, m_error_storage[0],
 152                 pid_op<>(steps, m_dt_storage[0], m_dt_storage[1], m_dt_storage[2], 0.7));
 153         }
 154
 155         value_type ratio = 1 / m_algebra.norm_inf(err);
 156
 157         value_type kappa = 1.0;
 158         ratio = 1.0 + kappa*atan((ratio - 1) / kappa);
 159
 160         if(ratio*dt >= m_dtmax)
 161         {
 162             ratio = m_dtmax / dt;
 163         }
 164
 165         if(ratio >= threshold())
 166         {
 167             m_error_storage.rotate();
 168             m_dt_storage.rotate();
 169
 170             ++m_init;
 171         }
 172         else
 173         {
 174             m_init = 0;
 175         }
 176
 177         return dt * static_cast<time_type>(ratio);
 178     };
 179
 180 private:
 181     algebra_type m_algebra;
 182
 183     time_type m_dtmax;
 184     error_storage_type m_error_storage;
 185     time_storage_type m_dt_storage;
 186
 187     size_t m_init;
 188     double m_abs_tol;
 189     double m_rel_tol;
 190
 191     coeff_type m_coeff;
 192 };
 193
 194 } // detail
 195 } // odeint
 196 } // numeric
 197 } // boost
 198
 199 #endif