ceph/src/boost/libs/numeric/odeint/include/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp

   1 /*
   2  [auto_generated]
   3  boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp
   4
   5  [begin_description]
   6  Base class for symplectic Runge-Kutta-Nystrom steppers.
   7  [end_description]
   8
   9  Copyright 2011-2013 Karsten Ahnert
  10  Copyright 2011-2013 Mario Mulansky
  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_BASE_SYMPLECTIC_RKN_STEPPER_BASE_HPP_INCLUDED
  20 #define BOOST_NUMERIC_ODEINT_STEPPER_BASE_SYMPLECTIC_RKN_STEPPER_BASE_HPP_INCLUDED
  21
  22 #include <boost/array.hpp>
  23
  24 #include <boost/numeric/odeint/util/bind.hpp>
  25 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
  26
  27 #include <boost/numeric/odeint/util/copy.hpp>
  28 #include <boost/numeric/odeint/util/is_pair.hpp>
  29
  30 #include <boost/numeric/odeint/util/state_wrapper.hpp>
  31 #include <boost/numeric/odeint/util/resizer.hpp>
  32
  33 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
  34
  35 #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
  36
  37
  38
  39
  40 namespace boost {
  41 namespace numeric {
  42 namespace odeint {
  43
  44
  45 template<
  46 size_t NumOfStages ,
  47 unsigned short Order ,
  48 class Coor ,
  49 class Momentum ,
  50 class Value ,
  51 class CoorDeriv ,
  52 class MomentumDeriv ,
  53 class Time ,
  54 class Algebra ,
  55 class Operations ,
  56 class Resizer
  57 >
  58 class symplectic_nystroem_stepper_base : public algebra_stepper_base< Algebra , Operations >
  59 {
  60
  61 public:
  62
  63     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
  64     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
  65     typedef typename algebra_stepper_base_type::operations_type operations_type;
  66
  67     const static size_t num_of_stages = NumOfStages;
  68     typedef Coor coor_type;
  69     typedef Momentum momentum_type;
  70     typedef std::pair< coor_type , momentum_type > state_type;
  71     typedef CoorDeriv coor_deriv_type;
  72     typedef state_wrapper< coor_deriv_type> wrapped_coor_deriv_type;
  73     typedef MomentumDeriv momentum_deriv_type;
  74     typedef state_wrapper< momentum_deriv_type > wrapped_momentum_deriv_type;
  75     typedef std::pair< coor_deriv_type , momentum_deriv_type > deriv_type;
  76     typedef Value value_type;
  77     typedef Time time_type;
  78     typedef Resizer resizer_type;
  79     typedef stepper_tag stepper_category;
  80
  81     #ifndef DOXYGEN_SKIP
  82     typedef symplectic_nystroem_stepper_base< NumOfStages , Order , Coor , Momentum , Value ,
  83             CoorDeriv , MomentumDeriv , Time , Algebra , Operations , Resizer > internal_stepper_base_type;
  84     #endif
  85     typedef unsigned short order_type;
  86
  87     static const order_type order_value = Order;
  88
  89     typedef boost::array< value_type , num_of_stages > coef_type;
  90
  91     symplectic_nystroem_stepper_base( const coef_type &coef_a , const coef_type &coef_b , const algebra_type &algebra = algebra_type() )
  92         : algebra_stepper_base_type( algebra ) , m_coef_a( coef_a ) , m_coef_b( coef_b ) ,
  93           m_dqdt_resizer() , m_dpdt_resizer() , m_dqdt() , m_dpdt()
  94     { }
  95
  96
  97     order_type order( void ) const
  98     {
  99         return order_value;
 100     }
 101
 102     /*
 103      * Version 1 : do_step( system , x , t , dt )
 104      *
 105      * This version does not solve the forwarding problem, boost.range can not be used.
 106      */
 107     template< class System , class StateInOut >
 108     void do_step( System system , const StateInOut &state , time_type t , time_type dt )
 109     {
 110         typedef typename odeint::unwrap_reference< System >::type system_type;
 111         do_step_impl( system , state , t , state , dt , typename is_pair< system_type >::type() );
 112     }
 113
 114     /**
 115      * \brief Same function as above. It differs only in a different const specifier in order
 116      * to solve the forwarding problem, can be used with Boost.Range.
 117      */
 118     template< class System , class StateInOut >
 119     void do_step( System system , StateInOut &state , time_type t , time_type dt )
 120     {
 121         typedef typename odeint::unwrap_reference< System >::type system_type;
 122         do_step_impl( system , state , t , state , dt , typename is_pair< system_type >::type() );
 123     }
 124
 125
 126
 127
 128     /*
 129      * Version 2 : do_step( system , q , p , t , dt );
 130      *
 131      * For Convenience
 132      *
 133      * The two overloads are needed in order to solve the forwarding problem.
 134      */
 135     template< class System , class CoorInOut , class MomentumInOut >
 136     void do_step( System system , CoorInOut &q , MomentumInOut &p , time_type t , time_type dt )
 137     {
 138         do_step( system , std::make_pair( detail::ref( q ) , detail::ref( p ) ) , t , dt );
 139     }
 140
 141     /**
 142      * \brief Same function as do_step( system , q , p , t , dt ). It differs only in a different const specifier in order
 143      * to solve the forwarding problem, can be called with Boost.Range.
 144      */
 145     template< class System , class CoorInOut , class MomentumInOut >
 146     void do_step( System system , const CoorInOut &q , const MomentumInOut &p , time_type t , time_type dt )
 147     {
 148         do_step( system , std::make_pair( detail::ref( q ) , detail::ref( p ) ) , t , dt );
 149     }
 150
 151
 152
 153
 154
 155     /*
 156      * Version 3 : do_step( system , in , t , out , dt )
 157      *
 158      * The forwarding problem is not solved in this version
 159      */
 160     template< class System , class StateIn , class StateOut >
 161     void do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
 162     {
 163         typedef typename odeint::unwrap_reference< System >::type system_type;
 164         do_step_impl( system , in , t , out , dt , typename is_pair< system_type >::type() );
 165     }
 166
 167
 168     template< class StateType >
 169     void adjust_size( const StateType &x )
 170     {
 171         resize_dqdt( x );
 172         resize_dpdt( x );
 173     }
 174
 175     /** \brief Returns the coefficients a. */
 176     const coef_type& coef_a( void ) const { return m_coef_a; }
 177
 178     /** \brief Returns the coefficients b. */
 179     const coef_type& coef_b( void ) const { return m_coef_b; }
 180
 181 private:
 182
 183     // stepper for systems with function for dq/dt = f(p) and dp/dt = -f(q)
 184     template< class System , class StateIn , class StateOut >
 185     void do_step_impl( System system , const StateIn &in , time_type /* t */ , StateOut &out , time_type dt , boost::mpl::true_ )
 186     {
 187         typedef typename odeint::unwrap_reference< System >::type system_type;
 188         typedef typename odeint::unwrap_reference< typename system_type::first_type >::type coor_deriv_func_type;
 189         typedef typename odeint::unwrap_reference< typename system_type::second_type >::type momentum_deriv_func_type;
 190         system_type &sys = system;
 191         coor_deriv_func_type &coor_func = sys.first;
 192         momentum_deriv_func_type &momentum_func = sys.second;
 193
 194         typedef typename odeint::unwrap_reference< StateIn >::type state_in_type;
 195         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
 196         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
 197         const state_in_type &state_in = in;
 198         const coor_in_type &coor_in = state_in.first;
 199         const momentum_in_type &momentum_in = state_in.second;
 200
 201         typedef typename odeint::unwrap_reference< StateOut >::type state_out_type;
 202         typedef typename odeint::unwrap_reference< typename state_out_type::first_type >::type coor_out_type;
 203         typedef typename odeint::unwrap_reference< typename state_out_type::second_type >::type momentum_out_type;
 204         state_out_type &state_out = out;
 205         coor_out_type &coor_out = state_out.first;
 206         momentum_out_type &momentum_out = state_out.second;
 207
 208         m_dqdt_resizer.adjust_size( coor_in , detail::bind( &internal_stepper_base_type::template resize_dqdt< coor_in_type > , detail::ref( *this ) , detail::_1 ) );
 209         m_dpdt_resizer.adjust_size( momentum_in , detail::bind( &internal_stepper_base_type::template resize_dpdt< momentum_in_type > , detail::ref( *this ) , detail::_1 ) );
 210
 211         // ToDo: check sizes?
 212
 213         for( size_t l=0 ; l<num_of_stages ; ++l )
 214         {
 215             if( l == 0 )
 216             {
 217                 coor_func( momentum_in , m_dqdt.m_v );
 218                 this->m_algebra.for_each3( coor_out , coor_in , m_dqdt.m_v ,
 219                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
 220                 momentum_func( coor_out , m_dpdt.m_v );
 221                 this->m_algebra.for_each3( momentum_out , momentum_in , m_dpdt.m_v ,
 222                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
 223             }
 224             else
 225             {
 226                 coor_func( momentum_out , m_dqdt.m_v );
 227                 this->m_algebra.for_each3( coor_out , coor_out , m_dqdt.m_v ,
 228                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
 229                 momentum_func( coor_out , m_dpdt.m_v );
 230                 this->m_algebra.for_each3( momentum_out , momentum_out , m_dpdt.m_v ,
 231                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
 232             }
 233         }
 234     }
 235
 236
 237     // stepper for systems with only function dp /dt = -f(q), dq/dt = p, time not required but still expected for compatibility reasons
 238     template< class System , class StateIn , class StateOut >
 239     void do_step_impl( System system , const StateIn &in , time_type  /* t */ , StateOut &out , time_type dt , boost::mpl::false_ )
 240     {
 241         typedef typename odeint::unwrap_reference< System >::type momentum_deriv_func_type;
 242         momentum_deriv_func_type &momentum_func = system;
 243
 244         typedef typename odeint::unwrap_reference< StateIn >::type state_in_type;
 245         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
 246         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
 247         const state_in_type &state_in = in;
 248         const coor_in_type &coor_in = state_in.first;
 249         const momentum_in_type &momentum_in = state_in.second;
 250
 251         typedef typename odeint::unwrap_reference< StateOut >::type state_out_type;
 252         typedef typename odeint::unwrap_reference< typename state_out_type::first_type >::type coor_out_type;
 253         typedef typename odeint::unwrap_reference< typename state_out_type::second_type >::type momentum_out_type;
 254         state_out_type &state_out = out;
 255         coor_out_type &coor_out = state_out.first;
 256         momentum_out_type &momentum_out = state_out.second;
 257
 258
 259         // m_dqdt not required when called with momentum_func only - don't resize
 260         // m_dqdt_resizer.adjust_size( coor_in , detail::bind( &internal_stepper_base_type::template resize_dqdt< coor_in_type > , detail::ref( *this ) , detail::_1 ) );
 261         m_dpdt_resizer.adjust_size( momentum_in , detail::bind( &internal_stepper_base_type::template resize_dpdt< momentum_in_type > , detail::ref( *this ) , detail::_1 ) );
 262
 263
 264         // ToDo: check sizes?
 265
 266         // step 0
 267         this->m_algebra.for_each3( coor_out  , coor_in , momentum_in ,
 268                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[0] * dt ) );
 269         momentum_func( coor_out , m_dpdt.m_v );
 270         this->m_algebra.for_each3( momentum_out , momentum_in , m_dpdt.m_v ,
 271                                            typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[0] * dt ) );
 272
 273         for( size_t l=1 ; l<num_of_stages ; ++l )
 274         {
 275             this->m_algebra.for_each3( coor_out , coor_out , momentum_out ,
 276                         typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_a[l] * dt ) );
 277             momentum_func( coor_out , m_dpdt.m_v );
 278             this->m_algebra.for_each3( momentum_out , momentum_out , m_dpdt.m_v ,
 279                                        typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , m_coef_b[l] * dt ) );
 280         }
 281     }
 282
 283     template< class StateIn >
 284     bool resize_dqdt( const StateIn &x )
 285     {
 286         return adjust_size_by_resizeability( m_dqdt , x , typename is_resizeable<coor_deriv_type>::type() );
 287     }
 288
 289     template< class StateIn >
 290     bool resize_dpdt( const StateIn &x )
 291     {
 292         return adjust_size_by_resizeability( m_dpdt , x , typename is_resizeable<momentum_deriv_type>::type() );
 293     }
 294
 295
 296     const coef_type m_coef_a;
 297     const coef_type m_coef_b;
 298
 299     resizer_type m_dqdt_resizer;
 300     resizer_type m_dpdt_resizer;
 301     wrapped_coor_deriv_type m_dqdt;
 302     wrapped_momentum_deriv_type m_dpdt;
 303
 304 };
 305
 306 /********* DOXYGEN *********/
 307
 308 /**
 309  * \class symplectic_nystroem_stepper_base
 310  * \brief Base class for all symplectic steppers of Nystroem type.
 311  *
 312  * This class is the base class for the symplectic Runge-Kutta-Nystroem steppers. Symplectic steppers are usually
 313  * used to solve Hamiltonian systems and they conserve the phase space volume, see
 314  * <a href="http://en.wikipedia.org/wiki/Symplectic_integrator">en.wikipedia.org/wiki/Symplectic_integrator</a>.
 315  * Furthermore, the energy is conserved
 316  * in average. In detail this class of steppers can be used to solve separable Hamiltonian systems which can be written
 317  * in the form H(q,p) = H1(p) + H2(q). q is usually called the coordinate, while p is the momentum. The equations of motion
 318  * are dq/dt = dH1/dp, dp/dt = -dH2/dq.
 319  *
 320  * ToDo : add formula for solver and explanation of the coefficients
 321  *
 322  * symplectic_nystroem_stepper_base uses odeints algebra and operation system. Step size and error estimation are not
 323  * provided for this class of solvers. It derives from algebra_stepper_base. Several `do_step` variants are provided:
 324  *
 325  * - `do_step( sys , x , t , dt )` - The classical `do_step` method. The sys can be either a pair of function objects
 326  *    for the coordinate or the momentum part or one function object for the momentum part. `x` is a pair of coordinate
 327  *    and momentum. The state is updated in-place.
 328  * - `do_step( sys , q , p , t , dt )` - This method is similar to the method above with the difference that the coordinate
 329  *    and the momentum are passed explicitly and not packed into a pair.
 330  * - `do_step( sys , x_in , t , x_out , dt )` - This method transforms the state out-of-place. `x_in` and `x_out` are here pairs
 331  *    of coordinate and momentum.
 332  *
 333  * \tparam NumOfStages Number of stages.
 334  * \tparam Order The order of the stepper.
 335  * \tparam Coor The type representing the coordinates q.
 336  * \tparam Momentum The type representing the coordinates p.
 337  * \tparam Value The basic value type. Should be something like float, double or a high-precision type.
 338  * \tparam CoorDeriv The type representing the time derivative of the coordinate dq/dt.
 339  * \tparam MomemtnumDeriv The type representing the time derivative of the momentum dp/dt.
 340  * \tparam Time The type representing the time t.
 341  * \tparam Algebra The algebra.
 342  * \tparam Operations The operations.
 343  * \tparam Resizer The resizer policy.
 344  */
 345
 346     /**
 347      * \fn symplectic_nystroem_stepper_base::symplectic_nystroem_stepper_base( const coef_type &coef_a , const coef_type &coef_b , const algebra_type &algebra )
 348      * \brief Constructs a symplectic_nystroem_stepper_base class. The parameters of the specific Nystroem method and the
 349      * algebra have to be passed.
 350      * \param coef_a The coefficients a.
 351      * \param coef_b The coefficients b.
 352      * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
 353      */
 354
 355     /**
 356      * \fn symplectic_nystroem_stepper_base::order( void ) const
 357      * \return Returns the order of the stepper.
 358      */
 359
 360     /**
 361      * \fn symplectic_nystroem_stepper_base::do_step( System system , const StateInOut &state , time_type t , time_type dt )
 362      * \brief This method performs one step. The system can be either a pair of two function object
 363      * describing the momentum part and the coordinate part or one function object describing only
 364      * the momentum part. In this case the coordinate is assumed to be trivial dq/dt = p. The state
 365      * is updated in-place.
 366      *
 367      * \note boost::ref or std::ref can be used for the system as well as for the state. So, it is correct
 368      * to write `stepper.do_step( make_pair( std::ref( fq ) , std::ref( fp ) ) , make_pair( std::ref( q ) , std::ref( p ) ) , t , dt )`.
 369      *
 370      * \note This method solves the forwarding problem.
 371      *
 372      * \param system The system, can be represented as a pair of two function object or one function object. See above.
 373      * \param state The state of the ODE. It is a pair of Coor and Momentum. The state is updated in-place, therefore, the
 374      * new value of the state will be written into this variable.
 375      * \param t The time of the ODE. It is not advanced by this method.
 376      * \param dt The time step.
 377      */
 378
 379     /**
 380      * \fn symplectic_nystroem_stepper_base::do_step( System system , CoorInOut &q , MomentumInOut &p , time_type t , time_type dt )
 381      * \brief This method performs one step. The system can be either a pair of two function object
 382      * describing the momentum part and the coordinate part or one function object describing only
 383      * the momentum part. In this case the coordinate is assumed to be trivial dq/dt = p. The state
 384      * is updated in-place.
 385      *
 386      * \note boost::ref or std::ref can be used for the system. So, it is correct
 387      * to write `stepper.do_step( make_pair( std::ref( fq ) , std::ref( fp ) ) , q , p , t , dt )`.
 388      *
 389      * \note This method solves the forwarding problem.
 390      *
 391      * \param system The system, can be represented as a pair of two function object or one function object. See above.
 392      * \param q The coordinate of the ODE. It is updated in-place. Therefore, the new value of the coordinate will be written
 393      * into this variable.
 394      * \param p The momentum of the ODE. It is updated in-place. Therefore, the new value of the momentum will be written info
 395      * this variable.
 396      * \param t The time of the ODE. It is not advanced by this method.
 397      * \param dt The time step.
 398      */
 399
 400     /**
 401      * \fn symplectic_nystroem_stepper_base::do_step( System system , const StateIn &in , time_type t , StateOut &out , time_type dt )
 402      * \brief This method performs one step. The system can be either a pair of two function object
 403      * describing the momentum part and the coordinate part or one function object describing only
 404      * the momentum part. In this case the coordinate is assumed to be trivial dq/dt = p. The state
 405      * is updated out-of-place.
 406      *
 407      * \note boost::ref or std::ref can be used for the system. So, it is correct
 408      * to write `stepper.do_step( make_pair( std::ref( fq ) , std::ref( fp ) ) , x_in , t , x_out , dt )`.
 409      *
 410      * \note This method NOT solve the forwarding problem.
 411      *
 412      * \param system The system, can be represented as a pair of two function object or one function object. See above.
 413      * \param in The state of the ODE, which is a pair of coordinate and momentum. The state is updated out-of-place, therefore the
 414      * new value is written into out
 415      * \param t The time of the ODE. It is not advanced by this method.
 416      * \param out The new state of the ODE.
 417      * \param dt The time step.
 418      */
 419
 420     /**
 421      * \fn symplectic_nystroem_stepper_base::adjust_size( const StateType &x )
 422      * \brief Adjust the size of all temporaries in the stepper manually.
 423      * \param x A state from which the size of the temporaries to be resized is deduced.
 424      */
 425
 426 } // namespace odeint
 427 } // namespace numeric
 428 } // namespace boost
 429
 430
 431 #endif // BOOST_NUMERIC_ODEINT_STEPPER_BASE_SYMPLECTIC_RKN_STEPPER_BASE_HPP_INCLUDED