]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 |