1 /* 2 [auto_generated] 3 boost/numeric/odeint/stepper/detail/generic_rk_algorithm.hpp 4 5 [begin_description] 6 Implementation of the generic Runge-Kutta method. 7 [end_description] 8 9 Copyright 2011-2013 Mario Mulansky 10 Copyright 2011-2012 Karsten Ahnert 11 Copyright 2012 Christoph Koke 12 13 Distributed under the Boost Software License, Version 1.0. 14 (See accompanying file LICENSE_1_0.txt or 15 copy at http://www.boost.org/LICENSE_1_0.txt) 16 */ 17 18 19 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_GENERIC_RK_ALGORITHM_HPP_INCLUDED 20 #define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_GENERIC_RK_ALGORITHM_HPP_INCLUDED 21 22 #include <boost/static_assert.hpp> 23 24 #include <boost/mpl/vector.hpp> 25 #include <boost/mpl/push_back.hpp> 26 #include <boost/mpl/for_each.hpp> 27 #include <boost/mpl/range_c.hpp> 28 #include <boost/mpl/copy.hpp> 29 #include <boost/mpl/size_t.hpp> 30 31 #include <boost/fusion/algorithm.hpp> 32 #include <boost/fusion/iterator.hpp> 33 #include <boost/fusion/mpl.hpp> 34 #include <boost/fusion/sequence.hpp> 35 36 #include <boost/array.hpp> 37 38 #include <boost/numeric/odeint/algebra/range_algebra.hpp> 39 #include <boost/numeric/odeint/algebra/default_operations.hpp> 40 #include <boost/numeric/odeint/stepper/detail/generic_rk_call_algebra.hpp> 41 #include <boost/numeric/odeint/stepper/detail/generic_rk_operations.hpp> 42 #include <boost/numeric/odeint/util/bind.hpp> 43 44 namespace boost { 45 namespace numeric { 46 namespace odeint { 47 namespace detail { 48 49 template< class T , class Constant > 50 struct array_wrapper 51 { 52 typedef const typename boost::array< T , Constant::value > type; 53 }; 54 55 template< class T , size_t i > 56 struct stage 57 { 58 T c; 59 boost::array< T , i > a; 60 }; 61 62 63 template< class T , class Constant > 64 struct stage_wrapper 65 { 66 typedef stage< T , Constant::value > type; 67 }; 68 69 70 template< 71 size_t StageCount, 72 class Value , 73 class Algebra , 74 class Operations 75 > 76 class generic_rk_algorithm { 77 78 public: 79 typedef mpl::range_c< size_t , 1 , StageCount > stage_indices; 80 81 typedef typename boost::fusion::result_of::as_vector 82 < 83 typename boost::mpl::copy 84 < 85 stage_indices , 86 boost::mpl::inserter 87 < 88 boost::mpl::vector0< > , 89 boost::mpl::push_back< boost::mpl::_1 , array_wrapper< Value , boost::mpl::_2 > > 90 > 91 >::type 92 >::type coef_a_type; 93 94 typedef boost::array< Value , StageCount > coef_b_type; 95 typedef boost::array< Value , StageCount > coef_c_type; 96 97 typedef typename boost::fusion::result_of::as_vector 98 < 99 typename boost::mpl::push_back 100 < 101 typename boost::mpl::copy 102 < 103 stage_indices, 104 boost::mpl::inserter 105 < 106 boost::mpl::vector0<> , 107 boost::mpl::push_back< boost::mpl::_1 , stage_wrapper< Value , boost::mpl::_2 > > 108 > 109 >::type , 110 stage< Value , StageCount > 111 >::type 112 >::type stage_vector_base; 113 114 115 struct stage_vector : public stage_vector_base 116 { 117 struct do_insertion 118 { 119 stage_vector_base &m_base; 120 const coef_a_type &m_a; 121 const coef_c_type &m_c; 122 do_insertionboost::numeric::odeint::detail::generic_rk_algorithm::stage_vector::do_insertion123 do_insertion( stage_vector_base &base , const coef_a_type &a , const coef_c_type &c ) 124 : m_base( base ) , m_a( a ) , m_c( c ) { } 125 126 template< class Index > operator ()boost::numeric::odeint::detail::generic_rk_algorithm::stage_vector::do_insertion127 void operator()( Index ) const 128 { 129 //boost::fusion::at< Index >( m_base ) = stage< double , Index::value+1 , intermediate_stage >( m_c[ Index::value ] , boost::fusion::at< Index >( m_a ) ); 130 boost::fusion::at< Index >( m_base ).c = m_c[ Index::value ]; 131 boost::fusion::at< Index >( m_base ).a = boost::fusion::at< Index >( m_a ); 132 } 133 }; 134 135 struct print_butcher 136 { 137 const stage_vector_base &m_base; 138 std::ostream &m_os; 139 print_butcherboost::numeric::odeint::detail::generic_rk_algorithm::stage_vector::print_butcher140 print_butcher( const stage_vector_base &base , std::ostream &os ) 141 : m_base( base ) , m_os( os ) 142 { } 143 144 template<class Index> operator ()boost::numeric::odeint::detail::generic_rk_algorithm::stage_vector::print_butcher145 void operator()(Index) const { 146 m_os << boost::fusion::at<Index>(m_base).c << " | "; 147 for( size_t i=0 ; i<Index::value ; ++i ) 148 m_os << boost::fusion::at<Index>(m_base).a[i] << " "; 149 m_os << std::endl; 150 } 151 }; 152 153 stage_vectorboost::numeric::odeint::detail::generic_rk_algorithm::stage_vector154 stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c ) 155 { 156 typedef boost::mpl::range_c< size_t , 0 , StageCount-1 > indices; 157 boost::mpl::for_each< indices >( do_insertion( *this , a , c ) ); 158 boost::fusion::at_c< StageCount - 1 >( *this ).c = c[ StageCount - 1 ]; 159 boost::fusion::at_c< StageCount - 1 >( *this ).a = b; 160 } 161 printboost::numeric::odeint::detail::generic_rk_algorithm::stage_vector162 void print( std::ostream &os ) const 163 { 164 typedef boost::mpl::range_c< size_t , 0 , StageCount > indices; 165 boost::mpl::for_each< indices >( print_butcher( *this , os ) ); 166 } 167 }; 168 169 170 171 template< class System , class StateIn , class StateTemp , class DerivIn , class Deriv , 172 class StateOut , class Time > 173 struct calculate_stage 174 { 175 Algebra &algebra; 176 System &system; 177 const StateIn &x; 178 StateTemp &x_tmp; 179 StateOut &x_out; 180 const DerivIn &dxdt; 181 Deriv *F; 182 Time t; 183 Time dt; 184 calculate_stageboost::numeric::odeint::detail::generic_rk_algorithm::calculate_stage185 calculate_stage( Algebra &_algebra , System &_system , const StateIn &_x , const DerivIn &_dxdt , StateOut &_out , 186 StateTemp &_x_tmp , Deriv *_F , Time _t , Time _dt ) 187 : algebra( _algebra ) , system( _system ) , x( _x ) , x_tmp( _x_tmp ) , x_out( _out) , dxdt( _dxdt ) , F( _F ) , t( _t ) , dt( _dt ) 188 {} 189 190 191 template< typename T , size_t stage_number > operator ()boost::numeric::odeint::detail::generic_rk_algorithm::calculate_stage192 void inline operator()( stage< T , stage_number > const &stage ) const 193 //typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , intermediate_stage >::type const &stage ) const 194 { 195 if( stage_number > 1 ) 196 { 197 #ifdef BOOST_MSVC 198 #pragma warning( disable : 4307 34 ) 199 #endif 200 system( x_tmp , F[stage_number-2].m_v , t + stage.c * dt ); 201 #ifdef BOOST_MSVC 202 #pragma warning( default : 4307 34 ) 203 #endif 204 } 205 //std::cout << stage_number-2 << ", t': " << t + stage.c * dt << std::endl; 206 207 if( stage_number < StageCount ) 208 detail::template generic_rk_call_algebra< stage_number , Algebra >()( algebra , x_tmp , x , dxdt , F , 209 detail::generic_rk_scale_sum< stage_number , Operations , Value , Time >( stage.a , dt) ); 210 // algebra_type::template for_eachn<stage_number>( x_tmp , x , dxdt , F , 211 // typename operations_type::template scale_sumn< stage_number , time_type >( stage.a , dt ) ); 212 else 213 detail::template generic_rk_call_algebra< stage_number , Algebra >()( algebra , x_out , x , dxdt , F , 214 detail::generic_rk_scale_sum< stage_number , Operations , Value , Time >( stage.a , dt ) ); 215 // algebra_type::template for_eachn<stage_number>( x_out , x , dxdt , F , 216 // typename operations_type::template scale_sumn< stage_number , time_type >( stage.a , dt ) ); 217 } 218 219 }; 220 generic_rk_algorithm(const coef_a_type & a,const coef_b_type & b,const coef_c_type & c)221 generic_rk_algorithm( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c ) 222 : m_stages( a , b , c ) 223 { } 224 225 template< class System , class StateIn , class DerivIn , class Time , class StateOut , class StateTemp , class Deriv > do_step(Algebra & algebra,System system,const StateIn & in,const DerivIn & dxdt,Time t,StateOut & out,Time dt,StateTemp & x_tmp,Deriv F[StageCount-1]) const226 void inline do_step( Algebra &algebra , System system , const StateIn &in , const DerivIn &dxdt , 227 Time t , StateOut &out , Time dt , 228 StateTemp &x_tmp , Deriv F[StageCount-1] ) const 229 { 230 typedef typename odeint::unwrap_reference< System >::type unwrapped_system_type; 231 unwrapped_system_type &sys = system; 232 boost::fusion::for_each( m_stages , calculate_stage< 233 unwrapped_system_type , StateIn , StateTemp , DerivIn , Deriv , StateOut , Time > 234 ( algebra , sys , in , dxdt , out , x_tmp , F , t , dt ) ); 235 } 236 237 private: 238 stage_vector m_stages; 239 }; 240 241 242 } 243 } 244 } 245 } 246 247 #endif // BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_GENERIC_RK_ALGORITHM_HPP_INCLUDED 248