stepper/detail/generic_rk_algorithm.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/stepper/detail/generic_rk_algorithm.hpp
0004
0005  [begin_description]
0006  Implementation of the generic Runge-Kutta method.
0007  [end_description]
0008
0009  Copyright 2011-2013 Mario Mulansky
0010  Copyright 2011-2012 Karsten Ahnert
0011  Copyright 2012 Christoph Koke
0012
0013  Distributed under the Boost Software License, Version 1.0.
0014  (See accompanying file LICENSE_1_0.txt or
0015  copy at http://www.boost.org/LICENSE_1_0.txt)
0016  */
0017
0018
0019 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_GENERIC_RK_ALGORITHM_HPP_INCLUDED
0020 #define BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_GENERIC_RK_ALGORITHM_HPP_INCLUDED
0021
0022 #include <boost/static_assert.hpp>
0023
0024 #include <boost/mpl/vector.hpp>
0025 #include <boost/mpl/push_back.hpp>
0026 #include <boost/mpl/for_each.hpp>
0027 #include <boost/mpl/range_c.hpp>
0028 #include <boost/mpl/copy.hpp>
0029 #include <boost/mpl/size_t.hpp>
0030
0031 #include <boost/fusion/algorithm.hpp>
0032 #include <boost/fusion/iterator.hpp>
0033 #include <boost/fusion/mpl.hpp>
0034 #include <boost/fusion/sequence.hpp>
0035
0036 #include <boost/array.hpp>
0037
0038 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
0039 #include <boost/numeric/odeint/algebra/default_operations.hpp>
0040 #include <boost/numeric/odeint/stepper/detail/generic_rk_call_algebra.hpp>
0041 #include <boost/numeric/odeint/stepper/detail/generic_rk_operations.hpp>
0042 #include <boost/numeric/odeint/util/bind.hpp>
0043
0044 namespace boost {
0045 namespace numeric {
0046 namespace odeint {
0047 namespace detail {
0048
0049 template< class T , class Constant >
0050 struct array_wrapper
0051 {
0052     typedef const typename boost::array< T , Constant::value > type;
0053 };
0054
0055 template< class T , size_t i >
0056 struct stage
0057 {
0058     T c;
0059     boost::array< T , i > a;
0060 };
0061
0062
0063 template< class T , class Constant >
0064 struct stage_wrapper
0065 {
0066     typedef stage< T , Constant::value > type;
0067 };
0068
0069
0070 template<
0071 size_t StageCount,
0072 class Value ,
0073 class Algebra ,
0074 class Operations
0075 >
0076 class generic_rk_algorithm {
0077
0078 public:
0079     typedef mpl::range_c< size_t , 1 , StageCount > stage_indices;
0080
0081     typedef typename boost::fusion::result_of::as_vector
0082             <
0083             typename boost::mpl::copy
0084             <
0085             stage_indices ,
0086             boost::mpl::inserter
0087             <
0088             boost::mpl::vector0< > ,
0089             boost::mpl::push_back< boost::mpl::_1 , array_wrapper< Value , boost::mpl::_2 > >
0090     >
0091     >::type
0092     >::type coef_a_type;
0093
0094     typedef boost::array< Value , StageCount > coef_b_type;
0095     typedef boost::array< Value , StageCount > coef_c_type;
0096
0097     typedef typename boost::fusion::result_of::as_vector
0098             <
0099             typename boost::mpl::push_back
0100             <
0101             typename boost::mpl::copy
0102             <
0103             stage_indices,
0104             boost::mpl::inserter
0105             <
0106             boost::mpl::vector0<> ,
0107             boost::mpl::push_back< boost::mpl::_1 , stage_wrapper< Value , boost::mpl::_2 > >
0108     >
0109     >::type ,
0110     stage< Value , StageCount >
0111     >::type
0112     >::type stage_vector_base;
0113
0114
0115     struct stage_vector : public stage_vector_base
0116     {
0117         struct do_insertion
0118         {
0119             stage_vector_base &m_base;
0120             const coef_a_type &m_a;
0121             const coef_c_type &m_c;
0122
0123             do_insertion( stage_vector_base &base , const coef_a_type &a , const coef_c_type &c )
0124             : m_base( base ) , m_a( a ) , m_c( c ) { }
0125
0126             template< class Index >
0127             void operator()( Index ) const
0128             {
0129                 //boost::fusion::at< Index >( m_base ) = stage< double , Index::value+1 , intermediate_stage >( m_c[ Index::value ] , boost::fusion::at< Index >( m_a ) );
0130                 boost::fusion::at< Index >( m_base ).c  = m_c[ Index::value ];
0131                 boost::fusion::at< Index >( m_base ).a = boost::fusion::at< Index >( m_a );
0132             }
0133         };
0134
0135         struct print_butcher
0136         {
0137             const stage_vector_base &m_base;
0138             std::ostream &m_os;
0139
0140             print_butcher( const stage_vector_base &base , std::ostream &os )
0141             : m_base( base ) , m_os( os )
0142             { }
0143
0144             template<class Index>
0145             void operator()(Index) const {
0146                 m_os << boost::fusion::at<Index>(m_base).c << " | ";
0147                 for( size_t i=0 ; i<Index::value ; ++i )
0148                     m_os << boost::fusion::at<Index>(m_base).a[i] << " ";
0149                 m_os << std::endl;
0150             }
0151         };
0152
0153
0154         stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
0155         {
0156             typedef boost::mpl::range_c< size_t , 0 , StageCount-1 > indices;
0157             boost::mpl::for_each< indices >( do_insertion( *this , a , c ) );
0158             boost::fusion::at_c< StageCount - 1 >( *this ).c = c[ StageCount - 1 ];
0159             boost::fusion::at_c< StageCount - 1 >( *this ).a = b;
0160         }
0161
0162         void print( std::ostream &os ) const
0163         {
0164             typedef boost::mpl::range_c< size_t , 0 , StageCount > indices;
0165             boost::mpl::for_each< indices >( print_butcher( *this , os ) );
0166         }
0167     };
0168
0169
0170
0171     template< class System , class StateIn , class StateTemp , class DerivIn , class Deriv ,
0172     class StateOut , class Time >
0173     struct calculate_stage
0174     {
0175         Algebra &algebra;
0176         System &system;
0177         const StateIn &x;
0178         StateTemp &x_tmp;
0179         StateOut &x_out;
0180         const DerivIn &dxdt;
0181         Deriv *F;
0182         Time t;
0183         Time dt;
0184
0185         calculate_stage( Algebra &_algebra , System &_system , const StateIn &_x , const DerivIn &_dxdt , StateOut &_out ,
0186                 StateTemp &_x_tmp , Deriv *_F , Time _t , Time _dt )
0187         : algebra( _algebra ) , system( _system ) , x( _x ) , x_tmp( _x_tmp ) , x_out( _out) , dxdt( _dxdt ) , F( _F ) , t( _t ) , dt( _dt )
0188         {}
0189
0190
0191         template< typename T , size_t stage_number >
0192         void inline operator()( stage< T , stage_number > const &stage ) const
0193         //typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , intermediate_stage >::type const &stage ) const
0194         {
0195             if( stage_number > 1 )
0196             {
0197 #ifdef BOOST_MSVC
0198 #pragma warning( disable : 4307 34 )
0199 #endif
0200                 system( x_tmp , F[stage_number-2].m_v , t + stage.c * dt );
0201 #ifdef BOOST_MSVC
0202 #pragma warning( default : 4307 34 )
0203 #endif
0204             }
0205             //std::cout << stage_number-2 << ", t': " << t + stage.c * dt << std::endl;
0206
0207             if( stage_number < StageCount )
0208                 detail::template generic_rk_call_algebra< stage_number , Algebra >()( algebra , x_tmp , x , dxdt , F ,
0209                         detail::generic_rk_scale_sum< stage_number , Operations , Value , Time >( stage.a , dt) );
0210             //                  algebra_type::template for_eachn<stage_number>( x_tmp , x , dxdt , F ,
0211             //                          typename operations_type::template scale_sumn< stage_number , time_type >( stage.a , dt ) );
0212             else
0213                 detail::template generic_rk_call_algebra< stage_number , Algebra >()( algebra , x_out , x , dxdt , F ,
0214                         detail::generic_rk_scale_sum< stage_number , Operations , Value , Time >( stage.a , dt ) );
0215             //                algebra_type::template for_eachn<stage_number>( x_out , x , dxdt , F ,
0216             //                            typename operations_type::template scale_sumn< stage_number , time_type >( stage.a , dt ) );
0217         }
0218
0219     };
0220
0221     generic_rk_algorithm( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
0222     : m_stages( a , b , c )
0223     { }
0224
0225     template< class System , class StateIn , class DerivIn , class Time , class StateOut , class StateTemp , class Deriv >
0226     void inline do_step( Algebra &algebra , System system , const StateIn &in , const DerivIn &dxdt ,
0227             Time t , StateOut &out , Time dt ,
0228             StateTemp &x_tmp , Deriv F[StageCount-1] ) const
0229     {
0230         typedef typename odeint::unwrap_reference< System >::type unwrapped_system_type;
0231         unwrapped_system_type &sys = system;
0232         boost::fusion::for_each( m_stages , calculate_stage<
0233                 unwrapped_system_type , StateIn , StateTemp , DerivIn , Deriv , StateOut , Time >
0234         ( algebra , sys , in , dxdt , out , x_tmp , F , t , dt ) );
0235     }
0236
0237 private:
0238     stage_vector m_stages;
0239 };
0240
0241
0242 }
0243 }
0244 }
0245 }
0246
0247 #endif // BOOST_NUMERIC_ODEINT_STEPPER_DETAIL_GENERIC_RK_ALGORITHM_HPP_INCLUDED