ceph/src/boost/boost/numeric/odeint/stepper/velocity_verlet.hpp

   1 /*
   2   [auto_generated]
   3   boost/numeric/odeint/stepper/velocity_verlet.hpp
   4
   5   [begin_description]
   6   tba.
   7   [end_description]
   8
   9   Copyright 2009-2012 Karsten Ahnert
  10   Copyright 2009-2012 Mario Mulansky
  11
  12   Distributed under the Boost Software License, Version 1.0.
  13   (See accompanying file LICENSE_1_0.txt or
  14   copy at http://www.boost.org/LICENSE_1_0.txt)
  15 */
  16
  17
  18 #ifndef BOOST_NUMERIC_ODEINT_STEPPER_VELOCITY_VERLET_HPP_DEFINED
  19 #define BOOST_NUMERIC_ODEINT_STEPPER_VELOCITY_VERLET_HPP_DEFINED
  20
  21 #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
  22 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
  23
  24 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
  25 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
  26 #include <boost/numeric/odeint/util/resizer.hpp>
  27 #include <boost/numeric/odeint/util/state_wrapper.hpp>
  28 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
  29
  30 #include <boost/numeric/odeint/util/bind.hpp>
  31 #include <boost/numeric/odeint/util/copy.hpp>
  32 #include <boost/numeric/odeint/util/resizer.hpp>
  33 // #include <boost/numeric/odeint/util/is_pair.hpp>
  34 // #include <boost/array.hpp>
  35
  36
  37
  38 namespace boost {
  39 namespace numeric {
  40 namespace odeint {
  41
  42
  43
  44 template <
  45     class Coor ,
  46     class Velocity = Coor ,
  47     class Value = double ,
  48     class Acceleration = Coor ,
  49     class Time = Value ,
  50     class TimeSq = Time ,
  51     class Algebra = typename algebra_dispatcher< Coor >::algebra_type ,
  52     class Operations = typename operations_dispatcher< Coor >::operations_type ,
  53     class Resizer = initially_resizer
  54       >
  55 class velocity_verlet : public algebra_stepper_base< Algebra , Operations >
  56 {
  57 public:
  58
  59     typedef algebra_stepper_base< Algebra , Operations > algebra_stepper_base_type;
  60     typedef typename algebra_stepper_base_type::algebra_type algebra_type;
  61     typedef typename algebra_stepper_base_type::operations_type operations_type;
  62
  63     typedef Coor coor_type;
  64     typedef Velocity velocity_type;
  65     typedef Acceleration acceleration_type;
  66     typedef std::pair< coor_type , velocity_type > state_type;
  67     typedef std::pair< velocity_type , acceleration_type > deriv_type;
  68     typedef state_wrapper< acceleration_type > wrapped_acceleration_type;
  69     typedef Value value_type;
  70     typedef Time time_type;
  71     typedef TimeSq time_square_type;
  72     typedef Resizer resizer_type;
  73     typedef stepper_tag stepper_category;
  74
  75     typedef unsigned short order_type;
  76
  77     static const order_type order_value = 1;
  78
  79     /**
  80      * \return Returns the order of the stepper.
  81      */
  82     order_type order( void ) const
  83     {
  84         return order_value;
  85     }
  86
  87
  88     velocity_verlet( const algebra_type & algebra = algebra_type() )
  89         : algebra_stepper_base_type( algebra ) , m_first_call( true )
  90         , m_a1() , m_a2() , m_current_a1( true ) { }
  91
  92
  93     template< class System , class StateInOut >
  94     void do_step( System system , StateInOut & x , time_type t , time_type dt )
  95     {
  96         do_step_v1( system , x , t , dt );
  97     }
  98
  99
 100     template< class System , class StateInOut >
 101     void do_step( System system , const StateInOut & x , time_type t , time_type dt )
 102     {
 103         do_step_v1( system , x , t , dt );
 104     }
 105
 106
 107     template< class System , class CoorIn , class VelocityIn , class AccelerationIn ,
 108                              class CoorOut , class VelocityOut , class AccelerationOut >
 109     void do_step( System system , CoorIn const & qin , VelocityIn const & pin , AccelerationIn const & ain ,
 110                   CoorOut & qout , VelocityOut & pout , AccelerationOut & aout , time_type t , time_type dt )
 111     {
 112         const value_type one = static_cast< value_type >( 1.0 );
 113         const value_type one_half = static_cast< value_type >( 0.5 );
 114
 115         algebra_stepper_base_type::m_algebra.for_each4(
 116             qout , qin , pin , ain ,
 117             typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one * dt , one_half * dt * dt ) );
 118
 119         typename odeint::unwrap_reference< System >::type & sys = system;
 120
 121         sys( qout , pin , aout , t + dt );
 122
 123         algebra_stepper_base_type::m_algebra.for_each4(
 124             pout , pin , ain , aout ,
 125             typename operations_type::template scale_sum3< value_type , time_type , time_type >( one , one_half * dt , one_half * dt ) );
 126     }
 127
 128
 129     template< class StateIn >
 130     void adjust_size( const StateIn & x )
 131     {
 132         if( resize_impl( x ) )
 133             m_first_call = true;
 134     }
 135
 136     void reset( void )
 137     {
 138         m_first_call = true;
 139     }
 140
 141
 142     /**
 143      * \fn velocity_verlet::initialize( const AccelerationIn &qin )
 144      * \brief Initializes the internal state of the stepper.
 145      * \param deriv The acceleration of x. The next call of `do_step` expects that the acceleration of `x` passed to `do_step`
 146      *              has the value of `qin`.
 147      */
 148     template< class AccelerationIn >
 149     void initialize( const AccelerationIn & ain )
 150     {
 151         // alloc a
 152         m_resizer.adjust_size( ain ,
 153                                detail::bind( &velocity_verlet::template resize_impl< AccelerationIn > ,
 154                                              detail::ref( *this ) , detail::_1 ) );
 155         boost::numeric::odeint::copy( ain , get_current_acc() );
 156         m_first_call = false;
 157     }
 158
 159
 160     template< class System , class CoorIn , class VelocityIn >
 161     void initialize( System system , const CoorIn & qin , const VelocityIn & pin , time_type t )
 162     {
 163         m_resizer.adjust_size( qin ,
 164                                detail::bind( &velocity_verlet::template resize_impl< CoorIn > ,
 165                                              detail::ref( *this ) , detail::_1 ) );
 166         initialize_acc( system , qin , pin , t );
 167     }
 168
 169     bool is_initialized( void ) const
 170     {
 171         return ! m_first_call;
 172     }
 173
 174
 175 private:
 176
 177     template< class System , class CoorIn , class VelocityIn >
 178     void initialize_acc( System system , const CoorIn & qin , const VelocityIn & pin , time_type t )
 179     {
 180         typename odeint::unwrap_reference< System >::type & sys = system;
 181         sys( qin , pin , get_current_acc() , t );
 182         m_first_call = false;
 183     }
 184
 185     template< class System , class StateInOut >
 186     void do_step_v1( System system , StateInOut & x , time_type t , time_type dt )
 187     {
 188         typedef typename odeint::unwrap_reference< StateInOut >::type state_in_type;
 189         typedef typename odeint::unwrap_reference< typename state_in_type::first_type >::type coor_in_type;
 190         typedef typename odeint::unwrap_reference< typename state_in_type::second_type >::type momentum_in_type;
 191
 192         typedef typename boost::remove_reference< coor_in_type >::type xyz_type;
 193         state_in_type & statein = x;
 194         coor_in_type & qinout = statein.first;
 195         momentum_in_type & pinout = statein.second;
 196
 197         // alloc a
 198         if( m_resizer.adjust_size( qinout ,
 199                                    detail::bind( &velocity_verlet::template resize_impl< xyz_type > ,
 200                                                  detail::ref( *this ) , detail::_1 ) )
 201          || m_first_call )
 202         {
 203             initialize_acc( system , qinout , pinout , t );
 204         }
 205
 206         // check first
 207         do_step( system , qinout , pinout , get_current_acc() , qinout , pinout , get_old_acc() , t , dt );
 208         toggle_current_acc();
 209     }
 210
 211     template< class StateIn >
 212     bool resize_impl( const StateIn & x )
 213     {
 214         bool resized = false;
 215         resized |= adjust_size_by_resizeability( m_a1 , x , typename is_resizeable< acceleration_type >::type() );
 216         resized |= adjust_size_by_resizeability( m_a2 , x , typename is_resizeable< acceleration_type >::type() );
 217         return resized;
 218     }
 219
 220     acceleration_type & get_current_acc( void )
 221     {
 222         return m_current_a1 ? m_a1.m_v : m_a2.m_v ;
 223     }
 224
 225     const acceleration_type & get_current_acc( void ) const
 226     {
 227         return m_current_a1 ? m_a1.m_v : m_a2.m_v ;
 228     }
 229
 230     acceleration_type & get_old_acc( void )
 231     {
 232         return m_current_a1 ? m_a2.m_v : m_a1.m_v ;
 233     }
 234
 235     const acceleration_type & get_old_acc( void ) const
 236     {
 237         return m_current_a1 ? m_a2.m_v : m_a1.m_v ;
 238     }
 239
 240     void toggle_current_acc( void )
 241     {
 242         m_current_a1 = ! m_current_a1;
 243     }
 244
 245     resizer_type m_resizer;
 246     bool m_first_call;
 247     wrapped_acceleration_type m_a1 , m_a2;
 248     bool m_current_a1;
 249 };
 250
 251 /**
 252  * \class velocity_verlet
 253  * \brief The Velocity-Verlet algorithm.
 254  *
 255  * <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
 256  * 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.
 257  *
 258  * \tparam Coor The type representing the coordinates.
 259  * \tparam Velocity The type representing the velocities.
 260  * \tparam Value The type value type.
 261  * \tparam Acceleration The type representing the acceleration.
 262  * \tparam Time The time representing the independent variable - the time.
 263  * \tparam TimeSq The time representing the square of the time.
 264  * \tparam Algebra The algebra.
 265  * \tparam Operations The operations type.
 266  * \tparam Resizer The resizer policy type.
 267  */
 268
 269
 270     /**
 271      * \fn velocity_verlet::velocity_verlet( const algebra_type &algebra )
 272      * \brief Constructs the velocity_verlet class. This constructor can be used as a default
 273      * constructor if the algebra has a default constructor.
 274      * \param algebra A copy of algebra is made and stored.
 275      */
 276
 277
 278     /**
 279      * \fn velocity_verlet::do_step( System system , StateInOut &x , time_type t , time_type dt )
 280      * \brief This method performs one step. It transforms the result in-place.
 281      *
 282      * It can be used like
 283      * \code
 284      * pair< coordinates , velocities > state;
 285      * stepper.do_step( sys , x , t , dt );
 286      * \endcode
 287      *
 288      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
 289      *               Second Order System concept.
 290      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
 291      * \param t The value of the time, at which the step should be performed.
 292      * \param dt The step size.
 293      */
 294
 295     /**
 296      * \fn velocity_verlet::do_step( System system , const StateInOut &x , time_type t , time_type dt )
 297      * \brief This method performs one step. It transforms the result in-place.
 298      *
 299      * It can be used like
 300      * \code
 301      * pair< coordinates , velocities > state;
 302      * stepper.do_step( sys , x , t , dt );
 303      * \endcode
 304      *
 305      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
 306      *               Second Order System concept.
 307      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
 308      * \param t The value of the time, at which the step should be performed.
 309      * \param dt The step size.
 310      */
 311
 312
 313
 314     /**
 315      * \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 )
 316      * \brief This method performs one step. It transforms the result in-place. Additionally to the other methods
 317      * the coordinates, velocities and accelerations are passed directly to do_step and they are transformed out-of-place.
 318      *
 319      * It can be used like
 320      * \code
 321      * coordinates qin , qout;
 322      * velocities pin , pout;
 323      * accelerations ain, aout;
 324      * stepper.do_step( sys , qin , pin , ain , qout , pout , aout , t , dt );
 325      * \endcode
 326      *
 327      * \param system The system function to solve, hence the r.h.s. of the ordinary differential equation. It must fulfill the
 328      *               Second Order System concept.
 329      * \param x The state of the ODE which should be solved. The state is pair of Coor and Velocity.
 330      * \param t The value of the time, at which the step should be performed.
 331      * \param dt The step size.
 332      */
 333
 334
 335     /**
 336      * \fn void velocity_verlet::adjust_size( const StateIn &x )
 337      * \brief Adjust the size of all temporaries in the stepper manually.
 338      * \param x A state from which the size of the temporaries to be resized is deduced.
 339      */
 340
 341
 342     /**
 343      * \fn velocity_verlet::reset( void )
 344      * \brief Resets the internal state of this stepper. After calling this method it is safe to use all
 345      * `do_step` method without explicitly initializing the stepper.
 346      */
 347
 348
 349
 350     /**
 351      * \fn velocity_verlet::initialize( System system , const CoorIn &qin , const VelocityIn &pin , time_type t )
 352      * \brief Initializes the internal state of the stepper.
 353      *
 354      * This method is equivalent to
 355      * \code
 356      * Acceleration a;
 357      * system( qin , pin , a , t );
 358      * stepper.initialize( a );
 359      * \endcode
 360      *
 361      * \param system The system function for the next calls of `do_step`.
 362      * \param qin The current coordinates of the ODE.
 363      * \param pin The current velocities of the ODE.
 364      * \param t The current time of the ODE.
 365      */
 366
 367
 368     /**
 369      * \fn velocity_verlet::is_initialized()
 370      * \returns Returns if the stepper is initialized.
 371     */
 372
 373
 374
 375
 376 } // namespace odeint
 377 } // namespace numeric
 378 } // namespace boost
 379
 380
 381 #endif // BOOST_NUMERIC_ODEINT_STEPPER_VELOCITY_VERLET_HPP_DEFINED