Chris@102: /*
Chris@102:   [auto_generated]
Chris@102:   boost/numeric/odeint/stepper/velocity_verlet.hpp
Chris@102: 
Chris@102:   [begin_description]
Chris@102:   tba.
Chris@102:   [end_description]
Chris@102: 
Chris@102:   Copyright 2009-2012 Karsten Ahnert
Chris@102:   Copyright 2009-2012 Mario Mulansky
Chris@102: 
Chris@102:   Distributed under the Boost Software License, Version 1.0.
Chris@102:   (See accompanying file LICENSE_1_0.txt or
Chris@102:   copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@102: */
Chris@102: 
Chris@102: 
Chris@102: #ifndef BOOST_NUMERIC_ODEINT_STEPPER_VELOCITY_VERLET_HPP_DEFINED
Chris@102: #define BOOST_NUMERIC_ODEINT_STEPPER_VELOCITY_VERLET_HPP_DEFINED
Chris@102: 
Chris@102: #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
Chris@102: #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
Chris@102: 
Chris@102: #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
Chris@102: #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
Chris@102: #include <boost/numeric/odeint/util/resizer.hpp>
Chris@102: #include <boost/numeric/odeint/util/state_wrapper.hpp>
Chris@102: #include <boost/numeric/odeint/util/unwrap_reference.hpp>
Chris@102: 
Chris@102: #include <boost/numeric/odeint/util/bind.hpp>
Chris@102: #include <boost/numeric/odeint/util/copy.hpp>
Chris@102: #include <boost/numeric/odeint/util/resizer.hpp>
Chris@102: // #include <boost/numeric/odeint/util/is_pair.hpp>
Chris@102: // #include <boost/array.hpp>
Chris@102: 
Chris@102: 
Chris@102: 
Chris@102: namespace boost {
Chris@102: namespace numeric {
Chris@102: namespace odeint {
Chris@102: 
Chris@102: 
Chris@102: 
Chris@102: template <
Chris@102:     class Coor ,
Chris@102:     class Velocity = Coor ,
Chris@102:     class Value = double ,
Chris@102:     class Acceleration = Coor ,
Chris@102:     class Time = Value ,
Chris@102:     class TimeSq = Time ,
Chris@102:     class Algebra = typename algebra_dispatcher< Coor >::algebra_type ,
Chris@102:     class Operations = typename operations_dispatcher< Coor >::operations_type ,
Chris@102:     class Resizer = initially_resizer
Chris@102:       >
Chris@102: class velocity_verlet : public algebra_stepper_base< Algebra , Operations >
Chris@102: {
Chris@102: public:
Chris@102: 
Chris@102:     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
Chris@102:     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
Chris@102:     typedef typename algebra_stepper_base_type::operations_type operations_type;
Chris@102: 
Chris@102:     typedef Coor coor_type;
Chris@102:     typedef Velocity velocity_type;
Chris@102:     typedef Acceleration acceleration_type;
Chris@102:     typedef std::pair< coor_type , velocity_type > state_type;
Chris@102:     typedef std::pair< velocity_type , acceleration_type > deriv_type;
Chris@102:     typedef state_wrapper< acceleration_type > wrapped_acceleration_type;
Chris@102:     typedef Value value_type;
Chris@102:     typedef Time time_type;
Chris@102:     typedef TimeSq time_square_type;
Chris@102:     typedef Resizer resizer_type;
Chris@102:     typedef stepper_tag stepper_category;
Chris@102: 
Chris@102:     typedef unsigned short order_type;
Chris@102: 
Chris@102:     static const order_type order_value = 1;
Chris@102: 
Chris@102:     /**
Chris@102:      * \return Returns the order of the stepper.
Chris@102:      */
Chris@102:     order_type order( void ) const
Chris@102:     {
Chris@102:         return order_value;
Chris@102:     }
Chris@102: 
Chris@102: 
Chris@102:     velocity_verlet( const algebra_type & algebra = algebra_type() )
Chris@102:         : algebra_stepper_base_type( algebra ) , m_first_call( true )
Chris@102:         , m_a1() , m_a2() , m_current_a1( true ) { }
Chris@102: 
Chris@102: 
Chris@102:     template< class System , class StateInOut >
Chris@102:     void do_step( System system , StateInOut & x , time_type t , time_type dt )
Chris@102:     {
Chris@102:         do_step_v1( system , x , t , dt );
Chris@102:     }
Chris@102: 
Chris@102:     
Chris@102:     template< class System , class StateInOut >
Chris@102:     void do_step( System system , const StateInOut & x , time_type t , time_type dt )
Chris@102:     {
Chris@102:         do_step_v1( system , x , t , dt );
Chris@102:     }
Chris@102: 
Chris@102:     
Chris@102:     template< class System , class CoorIn , class VelocityIn , class AccelerationIn ,
Chris@102:                              class CoorOut , class VelocityOut , class AccelerationOut >
Chris@102:     void do_step( System system , CoorIn const & qin , VelocityIn const & pin , AccelerationIn const & ain ,
Chris@102:                   CoorOut & qout , VelocityOut & pout , AccelerationOut & aout , time_type t , time_type dt )
Chris@102:     {
Chris@102:         const value_type one = static_cast< value_type >( 1.0 );
Chris@102:         const value_type one_half = static_cast< value_type >( 0.5 );
Chris@102: 
Chris@102:         algebra_stepper_base_type::m_algebra.for_each4(
Chris@102:             qout , qin , pin , ain ,
Chris@102:             typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one * dt , one_half * dt * dt ) );
Chris@102: 
Chris@102:         typename odeint::unwrap_reference< System >::type & sys = system;
Chris@102: 
Chris@102:         sys( qout , pin , aout , t + dt );
Chris@102: 
Chris@102:         algebra_stepper_base_type::m_algebra.for_each4(
Chris@102:             pout , pin , ain , aout ,
Chris@102:             typename operations_type::template scale_sum3< value_type , time_type , time_type >( one , one_half * dt , one_half * dt ) );
Chris@102:     }
Chris@102: 
Chris@102: 
Chris@102:     template< class StateIn >
Chris@102:     void adjust_size( const StateIn & x )
Chris@102:     {
Chris@102:         if( resize_impl( x ) )
Chris@102:             m_first_call = true;
Chris@102:     }
Chris@102: 
Chris@102:     void reset( void )
Chris@102:     {
Chris@102:         m_first_call = true;
Chris@102:     }
Chris@102: 
Chris@102:     
Chris@102:     /**
Chris@102:      * \fn velocity_verlet::initialize( const AccelerationIn &qin )
Chris@102:      * \brief Initializes the internal state of the stepper.
Chris@102:      * \param deriv The acceleration of x. The next call of `do_step` expects that the acceleration of `x` passed to `do_step`
Chris@102:      *              has the value of `qin`.
Chris@102:      */
Chris@102:     template< class AccelerationIn >
Chris@102:     void initialize( const AccelerationIn & ain )
Chris@102:     {
Chris@102:         // alloc a
Chris@102:         m_resizer.adjust_size( ain ,
Chris@102:                                detail::bind( &velocity_verlet::template resize_impl< AccelerationIn > ,
Chris@102:                                              detail::ref( *this ) , detail::_1 ) );
Chris@102:         boost::numeric::odeint::copy( ain , get_current_acc() );
Chris@102:         m_first_call = false;
Chris@102:     }
Chris@102: 
Chris@102: 
Chris@102:     template< class System , class CoorIn , class VelocityIn >
Chris@102:     void initialize( System system , const CoorIn & qin , const VelocityIn & pin , time_type t )
Chris@102:     {
Chris@102:         m_resizer.adjust_size( qin ,
Chris@102:                                detail::bind( &velocity_verlet::template resize_impl< CoorIn > ,
Chris@102:                                              detail::ref( *this ) , detail::_1 ) );
Chris@102:         initialize_acc( system , qin , pin , t );
Chris@102:     }
Chris@102: 
Chris@102:     bool is_initialized( void ) const
Chris@102:     {
Chris@102:         return ! m_first_call;
Chris@102:     }
Chris@102: 
Chris@102: 
Chris@102: private:
Chris@102:     
Chris@102:     template< class System , class CoorIn , class VelocityIn >
Chris@102:     void initialize_acc( System system , const CoorIn & qin , const VelocityIn & pin , time_type t )
Chris@102:     {
Chris@102:         typename odeint::unwrap_reference< System >::type & sys = system;
Chris@102:         sys( qin , pin , get_current_acc() , t );
Chris@102:         m_first_call = false;
Chris@102:     }
Chris@102:     
Chris@102:     template< class System , class StateInOut >
Chris@102:     void do_step_v1( System system , StateInOut & x , time_type t , time_type dt )
Chris@102:     {
Chris@102:         typedef typename odeint::unwrap_reference< StateInOut >::type state_in_type;
Chris@102:         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
Chris@102:         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
Chris@102:         
Chris@102:         typedef typename boost::remove_reference< coor_in_type >::type xyz_type;
Chris@102:         state_in_type & statein = x;
Chris@102:         coor_in_type & qinout = statein.first;
Chris@102:         momentum_in_type & pinout = statein.second;
Chris@102: 
Chris@102:         // alloc a
Chris@102:         if( m_resizer.adjust_size( qinout ,
Chris@102:                                    detail::bind( &velocity_verlet::template resize_impl< xyz_type > ,
Chris@102:                                                  detail::ref( *this ) , detail::_1 ) )
Chris@102:          || m_first_call )
Chris@102:         {
Chris@102:             initialize_acc( system , qinout , pinout , t );
Chris@102:         }
Chris@102: 
Chris@102:         // check first
Chris@102:         do_step( system , qinout , pinout , get_current_acc() , qinout , pinout , get_old_acc() , t , dt );
Chris@102:         toggle_current_acc();
Chris@102:     }
Chris@102: 
Chris@102:     template< class StateIn >
Chris@102:     bool resize_impl( const StateIn & x )
Chris@102:     {
Chris@102:         bool resized = false;
Chris@102:         resized |= adjust_size_by_resizeability( m_a1 , x , typename is_resizeable< acceleration_type >::type() );
Chris@102:         resized |= adjust_size_by_resizeability( m_a2 , x , typename is_resizeable< acceleration_type >::type() );
Chris@102:         return resized;
Chris@102:     }
Chris@102: 
Chris@102:     acceleration_type & get_current_acc( void )
Chris@102:     {
Chris@102:         return m_current_a1 ? m_a1.m_v : m_a2.m_v ;
Chris@102:     }
Chris@102: 
Chris@102:     const acceleration_type & get_current_acc( void ) const
Chris@102:     {
Chris@102:         return m_current_a1 ? m_a1.m_v : m_a2.m_v ;
Chris@102:     }
Chris@102: 
Chris@102:     acceleration_type & get_old_acc( void )
Chris@102:     {
Chris@102:         return m_current_a1 ? m_a2.m_v : m_a1.m_v ;
Chris@102:     }
Chris@102: 
Chris@102:     const acceleration_type & get_old_acc( void ) const
Chris@102:     {
Chris@102:         return m_current_a1 ? m_a2.m_v : m_a1.m_v ;
Chris@102:     }
Chris@102: 
Chris@102:     void toggle_current_acc( void )
Chris@102:     {
Chris@102:         m_current_a1 = ! m_current_a1;
Chris@102:     }
Chris@102: 
Chris@102:     resizer_type m_resizer;
Chris@102:     bool m_first_call;
Chris@102:     wrapped_acceleration_type m_a1 , m_a2;
Chris@102:     bool m_current_a1;
Chris@102: };
Chris@102: 
Chris@102: /**
Chris@102:  * \class velocity_verlet
Chris@102:  * \brief The Velocity-Verlet algorithm.
Chris@102:  *
Chris@102:  * <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
Chris@102:  * 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.
Chris@102:  * 
Chris@102:  * \tparam Coor The type representing the coordinates.
Chris@102:  * \tparam Velocity The type representing the velocities.
Chris@102:  * \tparam Value The type value type.
Chris@102:  * \tparam Acceleration The type representing the acceleration.
Chris@102:  * \tparam Time The time representing the independent variable - the time.
Chris@102:  * \tparam TimeSq The time representing the square of the time.
Chris@102:  * \tparam Algebra The algebra.
Chris@102:  * \tparam Operations The operations type.
Chris@102:  * \tparam Resizer The resizer policy type.
Chris@102:  */
Chris@102: 
Chris@102: 
Chris@102:     /**
Chris@102:      * \fn velocity_verlet::velocity_verlet( const algebra_type &algebra )
Chris@102:      * \brief Constructs the velocity_verlet class. This constructor can be used as a default
Chris@102:      * constructor if the algebra has a default constructor. 
Chris@102:      * \param algebra A copy of algebra is made and stored.
Chris@102:      */
Chris@102: 
Chris@102:     
Chris@102:     /**
Chris@102:      * \fn velocity_verlet::do_step( System system , StateInOut &x , time_type t , time_type dt )
Chris@102:      * \brief This method performs one step. It transforms the result in-place.
Chris@102:      * 
Chris@102:      * It can be used like
Chris@102:      * \code
Chris@102:      * pair< coordinates , velocities > state;
Chris@102:      * stepper.do_step( sys , x , t , dt );
Chris@102:      * \endcode
Chris@102:      *
Chris@102:      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
Chris@102:      *               Second Order System concept.
Chris@102:      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
Chris@102:      * \param t The value of the time, at which the step should be performed.
Chris@102:      * \param dt The step size.
Chris@102:      */
Chris@102: 
Chris@102:     /**
Chris@102:      * \fn velocity_verlet::do_step( System system , const StateInOut &x , time_type t , time_type dt )
Chris@102:      * \brief This method performs one step. It transforms the result in-place.
Chris@102:      * 
Chris@102:      * It can be used like
Chris@102:      * \code
Chris@102:      * pair< coordinates , velocities > state;
Chris@102:      * stepper.do_step( sys , x , t , dt );
Chris@102:      * \endcode
Chris@102:      *
Chris@102:      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
Chris@102:      *               Second Order System concept.
Chris@102:      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
Chris@102:      * \param t The value of the time, at which the step should be performed.
Chris@102:      * \param dt The step size.
Chris@102:      */    
Chris@102: 
Chris@102:     
Chris@102: 
Chris@102:     /**
Chris@102:      * \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 )
Chris@102:      * \brief This method performs one step. It transforms the result in-place. Additionally to the other methods
Chris@102:      * the coordinates, velocities and accelerations are passed directly to do_step and they are transformed out-of-place.
Chris@102:      * 
Chris@102:      * It can be used like
Chris@102:      * \code
Chris@102:      * coordinates qin , qout;
Chris@102:      * velocities pin , pout;
Chris@102:      * accelerations ain, aout;
Chris@102:      * stepper.do_step( sys , qin , pin , ain , qout , pout , aout , t , dt );
Chris@102:      * \endcode
Chris@102:      *
Chris@102:      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
Chris@102:      *               Second Order System concept.
Chris@102:      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
Chris@102:      * \param t The value of the time, at which the step should be performed.
Chris@102:      * \param dt The step size.
Chris@102:      */
Chris@102: 
Chris@102:     
Chris@102:     /**
Chris@102:      * \fn void velocity_verlet::adjust_size( const StateIn &x )
Chris@102:      * \brief Adjust the size of all temporaries in the stepper manually.
Chris@102:      * \param x A state from which the size of the temporaries to be resized is deduced.
Chris@102:      */
Chris@102: 
Chris@102: 
Chris@102:     /**
Chris@102:      * \fn velocity_verlet::reset( void )
Chris@102:      * \brief Resets the internal state of this stepper. After calling this method it is safe to use all
Chris@102:      * `do_step` method without explicitly initializing the stepper.
Chris@102:      */
Chris@102: 
Chris@102:     
Chris@102: 
Chris@102:     /**
Chris@102:      * \fn velocity_verlet::initialize( System system , const CoorIn &qin , const VelocityIn &pin , time_type t )
Chris@102:      * \brief Initializes the internal state of the stepper.
Chris@102:      *
Chris@102:      * This method is equivalent to 
Chris@102:      * \code
Chris@102:      * Acceleration a;
Chris@102:      * system( qin , pin , a , t );
Chris@102:      * stepper.initialize( a );
Chris@102:      * \endcode
Chris@102:      *
Chris@102:      * \param system The system function for the next calls of `do_step`.
Chris@102:      * \param qin The current coordinates of the ODE.
Chris@102:      * \param pin The current velocities of the ODE.
Chris@102:      * \param t The current time of the ODE.
Chris@102:      */
Chris@102:     
Chris@102:     
Chris@102:     /**
Chris@102:      * \fn velocity_verlet::is_initialized()
Chris@102:      * \returns Returns if the stepper is initialized.
Chris@102:     */
Chris@102:     
Chris@102:     
Chris@102:     
Chris@102:     
Chris@102: } // namespace odeint
Chris@102: } // namespace numeric
Chris@102: } // namespace boost
Chris@102: 
Chris@102: 
Chris@102: #endif // BOOST_NUMERIC_ODEINT_STEPPER_VELOCITY_VERLET_HPP_DEFINED