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 <array>
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, [this](auto&& arg) { return this->resize_impl<AccelerationIn>(std::forward<decltype(arg)>(arg)); });
0153         boost::numeric::odeint::copy( ain , get_current_acc() );
0154         m_first_call = false;
0155     }
0156
0157
0158     template< class System , class CoorIn , class VelocityIn >
0159     void initialize( System system , const CoorIn & qin , const VelocityIn & pin , time_type t )
0160     {
0161         m_resizer.adjust_size(qin, [this](auto&& arg) { return this->resize_impl<CoorIn>(std::forward<decltype(arg)>(arg)); });
0162         initialize_acc( system , qin , pin , t );
0163     }
0164
0165     bool is_initialized( void ) const
0166     {
0167         return ! m_first_call;
0168     }
0169
0170
0171 private:
0172
0173     template< class System , class CoorIn , class VelocityIn >
0174     void initialize_acc( System system , const CoorIn & qin , const VelocityIn & pin , time_type t )
0175     {
0176         typename odeint::unwrap_reference< System >::type & sys = system;
0177         sys( qin , pin , get_current_acc() , t );
0178         m_first_call = false;
0179     }
0180
0181     template< class System , class StateInOut >
0182     void do_step_v1( System system , StateInOut & x , time_type t , time_type dt )
0183     {
0184         typedef typename odeint::unwrap_reference< StateInOut >::type state_in_type;
0185         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
0186         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
0187
0188         typedef typename boost::remove_reference< coor_in_type >::type xyz_type;
0189         state_in_type & statein = x;
0190         coor_in_type & qinout = statein.first;
0191         momentum_in_type & pinout = statein.second;
0192
0193         // alloc a
0194         if( m_resizer.adjust_size(qinout, [this](auto&& arg) { return this->resize_impl<xyz_type>(std::forward<decltype(arg)>(arg)); })
0195            || m_first_call )
0196         {
0197             initialize_acc( system , qinout , pinout , t );
0198         }
0199
0200         // check first
0201         do_step( system , qinout , pinout , get_current_acc() , qinout , pinout , get_old_acc() , t , dt );
0202         toggle_current_acc();
0203     }
0204
0205     template< class StateIn >
0206     bool resize_impl( const StateIn & x )
0207     {
0208         bool resized = false;
0209         resized |= adjust_size_by_resizeability( m_a1 , x , typename is_resizeable< acceleration_type >::type() );
0210         resized |= adjust_size_by_resizeability( m_a2 , x , typename is_resizeable< acceleration_type >::type() );
0211         return resized;
0212     }
0213
0214     acceleration_type & get_current_acc( void )
0215     {
0216         return m_current_a1 ? m_a1.m_v : m_a2.m_v ;
0217     }
0218
0219     const acceleration_type & get_current_acc( void ) const
0220     {
0221         return m_current_a1 ? m_a1.m_v : m_a2.m_v ;
0222     }
0223
0224     acceleration_type & get_old_acc( void )
0225     {
0226         return m_current_a1 ? m_a2.m_v : m_a1.m_v ;
0227     }
0228
0229     const acceleration_type & get_old_acc( void ) const
0230     {
0231         return m_current_a1 ? m_a2.m_v : m_a1.m_v ;
0232     }
0233
0234     void toggle_current_acc( void )
0235     {
0236         m_current_a1 = ! m_current_a1;
0237     }
0238
0239     resizer_type m_resizer;
0240     bool m_first_call;
0241     wrapped_acceleration_type m_a1 , m_a2;
0242     bool m_current_a1;
0243 };
0244
0245 /**
0246  * \class velocity_verlet
0247  * \brief The Velocity-Verlet algorithm.
0248  *
0249  * <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
0250  * 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.
0251  *
0252  * \tparam Coor The type representing the coordinates.
0253  * \tparam Velocity The type representing the velocities.
0254  * \tparam Value The type value type.
0255  * \tparam Acceleration The type representing the acceleration.
0256  * \tparam Time The time representing the independent variable - the time.
0257  * \tparam TimeSq The time representing the square of the time.
0258  * \tparam Algebra The algebra.
0259  * \tparam Operations The operations type.
0260  * \tparam Resizer The resizer policy type.
0261  */
0262
0263
0264     /**
0265      * \fn velocity_verlet::velocity_verlet( const algebra_type &algebra )
0266      * \brief Constructs the velocity_verlet class. This constructor can be used as a default
0267      * constructor if the algebra has a default constructor.
0268      * \param algebra A copy of algebra is made and stored.
0269      */
0270
0271
0272     /**
0273      * \fn velocity_verlet::do_step( System system , StateInOut &x , time_type t , time_type dt )
0274      * \brief This method performs one step. It transforms the result in-place.
0275      *
0276      * It can be used like
0277      * \code
0278      * pair< coordinates , velocities > state;
0279      * stepper.do_step( sys , x , t , dt );
0280      * \endcode
0281      *
0282      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
0283      *               Second Order System concept.
0284      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
0285      * \param t The value of the time, at which the step should be performed.
0286      * \param dt The step size.
0287      */
0288
0289     /**
0290      * \fn velocity_verlet::do_step( System system , const StateInOut &x , time_type t , time_type dt )
0291      * \brief This method performs one step. It transforms the result in-place.
0292      *
0293      * It can be used like
0294      * \code
0295      * pair< coordinates , velocities > state;
0296      * stepper.do_step( sys , x , t , dt );
0297      * \endcode
0298      *
0299      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
0300      *               Second Order System concept.
0301      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
0302      * \param t The value of the time, at which the step should be performed.
0303      * \param dt The step size.
0304      */
0305
0306
0307
0308     /**
0309      * \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 )
0310      * \brief This method performs one step. It transforms the result in-place. Additionally to the other methods
0311      * the coordinates, velocities and accelerations are passed directly to do_step and they are transformed out-of-place.
0312      *
0313      * It can be used like
0314      * \code
0315      * coordinates qin , qout;
0316      * velocities pin , pout;
0317      * accelerations ain, aout;
0318      * stepper.do_step( sys , qin , pin , ain , qout , pout , aout , t , dt );
0319      * \endcode
0320      *
0321      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
0322      *               Second Order System concept.
0323      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
0324      * \param t The value of the time, at which the step should be performed.
0325      * \param dt The step size.
0326      */
0327
0328
0329     /**
0330      * \fn void velocity_verlet::adjust_size( const StateIn &x )
0331      * \brief Adjust the size of all temporaries in the stepper manually.
0332      * \param x A state from which the size of the temporaries to be resized is deduced.
0333      */
0334
0335
0336     /**
0337      * \fn velocity_verlet::reset( void )
0338      * \brief Resets the internal state of this stepper. After calling this method it is safe to use all
0339      * `do_step` method without explicitly initializing the stepper.
0340      */
0341
0342
0343
0344     /**
0345      * \fn velocity_verlet::initialize( System system , const CoorIn &qin , const VelocityIn &pin , time_type t )
0346      * \brief Initializes the internal state of the stepper.
0347      *
0348      * This method is equivalent to
0349      * \code
0350      * Acceleration a;
0351      * system( qin , pin , a , t );
0352      * stepper.initialize( a );
0353      * \endcode
0354      *
0355      * \param system The system function for the next calls of `do_step`.
0356      * \param qin The current coordinates of the ODE.
0357      * \param pin The current velocities of the ODE.
0358      * \param t The current time of the ODE.
0359      */
0360
0361
0362     /**
0363      * \fn velocity_verlet::is_initialized()
0364      * \returns Returns if the stepper is initialized.
0365     */
0366
0367
0368
0369
0370 } // namespace odeint
0371 } // namespace numeric
0372 } // namespace boost
0373
0374
0375 #endif // BOOST_NUMERIC_ODEINT_STEPPER_VELOCITY_VERLET_HPP_DEFINED