odeint/stepper/velocity_verlet.hpp

0001 /*
0002   [auto_generated]
0003   boost/numeric/odeint/stepper/velocity_verlet.hpp
0004
0005   [begin_description]
0006   tba.
0007   [end_description]
0008
0009   Copyright 2009-2012 Karsten Ahnert
0010   Copyright 2009-2012 Mario Mulansky
0011
0012   Distributed under the Boost Software License, Version 1.0.
0013   (See accompanying file LICENSE_1_0.txt or
0014   copy at http://www.boost.org/LICENSE_1_0.txt)
0015 */
0016
0017
0018 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_VELOCITY_VERLET_HPP_DEFINED
0019 #define BOOST_NUMERIC_ODEINT_STEPPER_VELOCITY_VERLET_HPP_DEFINED
0020
0021 #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
0022 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
0023
0024 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
0025 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
0026 #include <boost/numeric/odeint/util/resizer.hpp>
0027 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0028 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
0029
0030 #include <boost/numeric/odeint/util/bind.hpp>
0031 #include <boost/numeric/odeint/util/copy.hpp>
0032 #include <boost/numeric/odeint/util/resizer.hpp>
0033 // #include <boost/numeric/odeint/util/is_pair.hpp>
0034 // #include <boost/array.hpp>
0035
0036
0037
0038 namespace boost {
0039 namespace numeric {
0040 namespace odeint {
0041
0042
0043
0044 template <
0045     class Coor ,
0046     class Velocity = Coor ,
0047     class Value = double ,
0048     class Acceleration = Coor ,
0049     class Time = Value ,
0050     class TimeSq = Time ,
0051     class Algebra = typename algebra_dispatcher< Coor >::algebra_type ,
0052     class Operations = typename operations_dispatcher< Coor >::operations_type ,
0053     class Resizer = initially_resizer
0054       >
0055 class velocity_verlet : public algebra_stepper_base< Algebra , Operations >
0056 {
0057 public:
0058
0059     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
0060     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
0061     typedef typename algebra_stepper_base_type::operations_type operations_type;
0062
0063     typedef Coor coor_type;
0064     typedef Velocity velocity_type;
0065     typedef Acceleration acceleration_type;
0066     typedef std::pair< coor_type , velocity_type > state_type;
0067     typedef std::pair< velocity_type , acceleration_type > deriv_type;
0068     typedef state_wrapper< acceleration_type > wrapped_acceleration_type;
0069     typedef Value value_type;
0070     typedef Time time_type;
0071     typedef TimeSq time_square_type;
0072     typedef Resizer resizer_type;
0073     typedef stepper_tag stepper_category;
0074
0075     typedef unsigned short order_type;
0076
0077     static const order_type order_value = 1;
0078
0079     /**
0080      * \return Returns the order of the stepper.
0081      */
0082     order_type order( void ) const
0083     {
0084         return order_value;
0085     }
0086
0087
0088     velocity_verlet( const algebra_type & algebra = algebra_type() )
0089         : algebra_stepper_base_type( algebra ) , m_first_call( true )
0090         , m_a1() , m_a2() , m_current_a1( true ) { }
0091
0092
0093     template< class System , class StateInOut >
0094     void do_step( System system , StateInOut & x , time_type t , time_type dt )
0095     {
0096         do_step_v1( system , x , t , dt );
0097     }
0098
0099
0100     template< class System , class StateInOut >
0101     void do_step( System system , const StateInOut & x , time_type t , time_type dt )
0102     {
0103         do_step_v1( system , x , t , dt );
0104     }
0105
0106
0107     template< class System , class CoorIn , class VelocityIn , class AccelerationIn ,
0108                              class CoorOut , class VelocityOut , class AccelerationOut >
0109     void do_step( System system , CoorIn const & qin , VelocityIn const & pin , AccelerationIn const & ain ,
0110                   CoorOut & qout , VelocityOut & pout , AccelerationOut & aout , time_type t , time_type dt )
0111     {
0112         const value_type one = static_cast< value_type >( 1.0 );
0113         const value_type one_half = static_cast< value_type >( 0.5 );
0114
0115         algebra_stepper_base_type::m_algebra.for_each4(
0116             qout , qin , pin , ain ,
0117             typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one * dt , one_half * dt * dt ) );
0118
0119         typename odeint::unwrap_reference< System >::type & sys = system;
0120
0121         sys( qout , pin , aout , t + dt );
0122
0123         algebra_stepper_base_type::m_algebra.for_each4(
0124             pout , pin , ain , aout ,
0125             typename operations_type::template scale_sum3< value_type , time_type , time_type >( one , one_half * dt , one_half * dt ) );
0126     }
0127
0128
0129     template< class StateIn >
0130     void adjust_size( const StateIn & x )
0131     {
0132         if( resize_impl( x ) )
0133             m_first_call = true;
0134     }
0135
0136     void reset( void )
0137     {
0138         m_first_call = true;
0139     }
0140
0141
0142     /**
0143      * \fn velocity_verlet::initialize( const AccelerationIn &qin )
0144      * \brief Initializes the internal state of the stepper.
0145      * \param deriv The acceleration of x. The next call of `do_step` expects that the acceleration of `x` passed to `do_step`
0146      *              has the value of `qin`.
0147      */
0148     template< class AccelerationIn >
0149     void initialize( const AccelerationIn & ain )
0150     {
0151         // alloc a
0152         m_resizer.adjust_size( ain ,
0153                                detail::bind( &velocity_verlet::template resize_impl< AccelerationIn > ,
0154                                              detail::ref( *this ) , detail::_1 ) );
0155         boost::numeric::odeint::copy( ain , get_current_acc() );
0156         m_first_call = false;
0157     }
0158
0159
0160     template< class System , class CoorIn , class VelocityIn >
0161     void initialize( System system , const CoorIn & qin , const VelocityIn & pin , time_type t )
0162     {
0163         m_resizer.adjust_size( qin ,
0164                                detail::bind( &velocity_verlet::template resize_impl< CoorIn > ,
0165                                              detail::ref( *this ) , detail::_1 ) );
0166         initialize_acc( system , qin , pin , t );
0167     }
0168
0169     bool is_initialized( void ) const
0170     {
0171         return ! m_first_call;
0172     }
0173
0174
0175 private:
0176
0177     template< class System , class CoorIn , class VelocityIn >
0178     void initialize_acc( System system , const CoorIn & qin , const VelocityIn & pin , time_type t )
0179     {
0180         typename odeint::unwrap_reference< System >::type & sys = system;
0181         sys( qin , pin , get_current_acc() , t );
0182         m_first_call = false;
0183     }
0184
0185     template< class System , class StateInOut >
0186     void do_step_v1( System system , StateInOut & x , time_type t , time_type dt )
0187     {
0188         typedef typename odeint::unwrap_reference< StateInOut >::type state_in_type;
0189         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
0190         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
0191
0192         typedef typename boost::remove_reference< coor_in_type >::type xyz_type;
0193         state_in_type & statein = x;
0194         coor_in_type & qinout = statein.first;
0195         momentum_in_type & pinout = statein.second;
0196
0197         // alloc a
0198         if( m_resizer.adjust_size( qinout ,
0199                                    detail::bind( &velocity_verlet::template resize_impl< xyz_type > ,
0200                                                  detail::ref( *this ) , detail::_1 ) )
0201          || m_first_call )
0202         {
0203             initialize_acc( system , qinout , pinout , t );
0204         }
0205
0206         // check first
0207         do_step( system , qinout , pinout , get_current_acc() , qinout , pinout , get_old_acc() , t , dt );
0208         toggle_current_acc();
0209     }
0210
0211     template< class StateIn >
0212     bool resize_impl( const StateIn & x )
0213     {
0214         bool resized = false;
0215         resized |= adjust_size_by_resizeability( m_a1 , x , typename is_resizeable< acceleration_type >::type() );
0216         resized |= adjust_size_by_resizeability( m_a2 , x , typename is_resizeable< acceleration_type >::type() );
0217         return resized;
0218     }
0219
0220     acceleration_type & get_current_acc( void )
0221     {
0222         return m_current_a1 ? m_a1.m_v : m_a2.m_v ;
0223     }
0224
0225     const acceleration_type & get_current_acc( void ) const
0226     {
0227         return m_current_a1 ? m_a1.m_v : m_a2.m_v ;
0228     }
0229
0230     acceleration_type & get_old_acc( void )
0231     {
0232         return m_current_a1 ? m_a2.m_v : m_a1.m_v ;
0233     }
0234
0235     const acceleration_type & get_old_acc( void ) const
0236     {
0237         return m_current_a1 ? m_a2.m_v : m_a1.m_v ;
0238     }
0239
0240     void toggle_current_acc( void )
0241     {
0242         m_current_a1 = ! m_current_a1;
0243     }
0244
0245     resizer_type m_resizer;
0246     bool m_first_call;
0247     wrapped_acceleration_type m_a1 , m_a2;
0248     bool m_current_a1;
0249 };
0250
0251 /**
0252  * \class velocity_verlet
0253  * \brief The Velocity-Verlet algorithm.
0254  *
0255  * <a href="http://en.wikipedia.org/wiki/Verlet_integration" >The Velocity-Verlet algorithm</a> is a method for simulation of molecular dynamics systems. It solves the ODE
0256  * a=f(r,v',t)  where r are the coordinates, v are the velocities and a are the accelerations, hence v = dr/dt, a=dv/dt.
0257  *
0258  * \tparam Coor The type representing the coordinates.
0259  * \tparam Velocity The type representing the velocities.
0260  * \tparam Value The type value type.
0261  * \tparam Acceleration The type representing the acceleration.
0262  * \tparam Time The time representing the independent variable - the time.
0263  * \tparam TimeSq The time representing the square of the time.
0264  * \tparam Algebra The algebra.
0265  * \tparam Operations The operations type.
0266  * \tparam Resizer The resizer policy type.
0267  */
0268
0269
0270     /**
0271      * \fn velocity_verlet::velocity_verlet( const algebra_type &algebra )
0272      * \brief Constructs the velocity_verlet class. This constructor can be used as a default
0273      * constructor if the algebra has a default constructor.
0274      * \param algebra A copy of algebra is made and stored.
0275      */
0276
0277
0278     /**
0279      * \fn velocity_verlet::do_step( System system , StateInOut &x , time_type t , time_type dt )
0280      * \brief This method performs one step. It transforms the result in-place.
0281      *
0282      * It can be used like
0283      * \code
0284      * pair< coordinates , velocities > state;
0285      * stepper.do_step( sys , x , t , dt );
0286      * \endcode
0287      *
0288      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
0289      *               Second Order System concept.
0290      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
0291      * \param t The value of the time, at which the step should be performed.
0292      * \param dt The step size.
0293      */
0294
0295     /**
0296      * \fn velocity_verlet::do_step( System system , const StateInOut &x , time_type t , time_type dt )
0297      * \brief This method performs one step. It transforms the result in-place.
0298      *
0299      * It can be used like
0300      * \code
0301      * pair< coordinates , velocities > state;
0302      * stepper.do_step( sys , x , t , dt );
0303      * \endcode
0304      *
0305      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
0306      *               Second Order System concept.
0307      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
0308      * \param t The value of the time, at which the step should be performed.
0309      * \param dt The step size.
0310      */
0311
0312
0313
0314     /**
0315      * \fn velocity_verlet::do_step( System system , CoorIn const & qin , VelocityIn const & pin , AccelerationIn const & ain , CoorOut & qout , VelocityOut & pout , AccelerationOut & aout , time_type t , time_type dt )
0316      * \brief This method performs one step. It transforms the result in-place. Additionally to the other methods
0317      * the coordinates, velocities and accelerations are passed directly to do_step and they are transformed out-of-place.
0318      *
0319      * It can be used like
0320      * \code
0321      * coordinates qin , qout;
0322      * velocities pin , pout;
0323      * accelerations ain, aout;
0324      * stepper.do_step( sys , qin , pin , ain , qout , pout , aout , t , dt );
0325      * \endcode
0326      *
0327      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
0328      *               Second Order System concept.
0329      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
0330      * \param t The value of the time, at which the step should be performed.
0331      * \param dt The step size.
0332      */
0333
0334
0335     /**
0336      * \fn void velocity_verlet::adjust_size( const StateIn &x )
0337      * \brief Adjust the size of all temporaries in the stepper manually.
0338      * \param x A state from which the size of the temporaries to be resized is deduced.
0339      */
0340
0341
0342     /**
0343      * \fn velocity_verlet::reset( void )
0344      * \brief Resets the internal state of this stepper. After calling this method it is safe to use all
0345      * `do_step` method without explicitly initializing the stepper.
0346      */
0347
0348
0349
0350     /**
0351      * \fn velocity_verlet::initialize( System system , const CoorIn &qin , const VelocityIn &pin , time_type t )
0352      * \brief Initializes the internal state of the stepper.
0353      *
0354      * This method is equivalent to
0355      * \code
0356      * Acceleration a;
0357      * system( qin , pin , a , t );
0358      * stepper.initialize( a );
0359      * \endcode
0360      *
0361      * \param system The system function for the next calls of `do_step`.
0362      * \param qin The current coordinates of the ODE.
0363      * \param pin The current velocities of the ODE.
0364      * \param t The current time of the ODE.
0365      */
0366
0367
0368     /**
0369      * \fn velocity_verlet::is_initialized()
0370      * \returns Returns if the stepper is initialized.
0371     */
0372
0373
0374
0375
0376 } // namespace odeint
0377 } // namespace numeric
0378 } // namespace boost
0379
0380
0381 #endif // BOOST_NUMERIC_ODEINT_STEPPER_VELOCITY_VERLET_HPP_DEFINED