1 /*
2  [auto_generated]
3  libs/numeric/odeint/test/symplectic_steppers.cpp
4 
5  [begin_description]
6  tba.
7  [end_description]
8 
9  Copyright 2012 Karsten Ahnert
10  Copyright 2013 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 #include <boost/config.hpp>
18 #ifdef BOOST_MSVC
19     #pragma warning(disable:4996)
20 #endif
21 
22 #define BOOST_TEST_MODULE odeint_symplectic_steppers
23 
24 #define BOOST_FUSION_INVOKE_MAX_ARITY 15
25 #define BOOST_RESULT_OF_NUM_ARGS 15
26 #define FUSION_MAX_VECTOR_SIZE 15
27 
28 
29 #include <boost/test/unit_test.hpp>
30 
31 #include <boost/array.hpp>
32 #include <boost/static_assert.hpp>
33 #include <boost/type_traits/is_same.hpp>
34 
35 #include <boost/mpl/vector.hpp>
36 #include <boost/mpl/zip_view.hpp>
37 #include <boost/mpl/vector_c.hpp>
38 #include <boost/mpl/insert_range.hpp>
39 #include <boost/mpl/end.hpp>
40 #include <boost/mpl/size.hpp>
41 #include <boost/mpl/copy.hpp>
42 #include <boost/mpl/placeholders.hpp>
43 #include <boost/mpl/inserter.hpp>
44 #include <boost/mpl/at.hpp>
45 
46 #include <boost/fusion/container/vector.hpp>
47 #include <boost/fusion/include/make_vector.hpp>
48 
49 #include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
50 #include <boost/numeric/odeint/algebra/fusion_algebra.hpp>
51 #include <boost/numeric/odeint/stepper/symplectic_euler.hpp>
52 #include <boost/numeric/odeint/stepper/symplectic_rkn_sb3a_mclachlan.hpp>
53 #include <boost/numeric/odeint/stepper/symplectic_rkn_sb3a_m4_mclachlan.hpp>
54 #include <boost/numeric/odeint/integrate/integrate_const.hpp>
55 
56 #include "diagnostic_state_type.hpp"
57 #include "const_range.hpp"
58 #include "dummy_odes.hpp"
59 #include "boost_units_helpers.hpp"
60 
61 
62 using namespace boost::unit_test;
63 using namespace boost::numeric::odeint;
64 namespace mpl = boost::mpl;
65 namespace fusion = boost::fusion;
66 
67 class custom_range_algebra : public range_algebra { };
68 class custom_default_operations : public default_operations { };
69 
70 
71 template< class Coor , class Mom , class Value , class CoorDeriv , class MomDeriv , class Time ,
72           class Algebra , class Operations , class Resizer >
73 class complete_steppers : public mpl::vector<
74     symplectic_euler< Coor , Mom , Value , CoorDeriv , MomDeriv , Time ,
75                       Algebra , Operations , Resizer >
76     , symplectic_rkn_sb3a_mclachlan< Coor , Mom , Value , CoorDeriv , MomDeriv , Time ,
77                                      Algebra , Operations , Resizer >
78     , symplectic_rkn_sb3a_m4_mclachlan< Coor , Mom , Value , CoorDeriv , MomDeriv , Time ,
79                                         Algebra , Operations , Resizer >
80     > {};
81 
82 template< class Resizer >
83 class vector_steppers : public complete_steppers<
84     diagnostic_state_type , diagnostic_state_type2 , double ,
85     diagnostic_deriv_type , diagnostic_deriv_type2 , double ,
86     custom_range_algebra , custom_default_operations , Resizer
87     > { };
88 
89 
90 typedef mpl::vector< initially_resizer , always_resizer , never_resizer > resizers;
91 typedef mpl::vector_c< int , 1 , 3 , 0 > resizer_multiplicities ;
92 
93 
94 typedef mpl::copy<
95     resizers ,
96     mpl::inserter<
97         mpl::vector0<> ,
98         mpl::insert_range<
99             mpl::_1 ,
100             mpl::end< mpl::_1 > ,
101             vector_steppers< mpl::_2 >
102             >
103         >
104     >::type all_stepper_methods;
105 
106 
107 typedef mpl::size< vector_steppers< initially_resizer > >::type num_steppers;
108 typedef mpl::copy<
109     resizer_multiplicities ,
110     mpl::inserter<
111         mpl::vector0<> ,
112         mpl::insert_range<
113             mpl::_1 ,
114             mpl::end< mpl::_1 > ,
115             const_range< num_steppers , mpl::_2 >
116             >
117         >
118     >::type all_multiplicities;
119 
120 
121 
122 
123 
124 
125 
126 BOOST_AUTO_TEST_SUITE( symplectic_steppers_test )
127 
128 
BOOST_AUTO_TEST_CASE_TEMPLATE(test_assoc_types,Stepper,vector_steppers<initially_resizer>)129 BOOST_AUTO_TEST_CASE_TEMPLATE( test_assoc_types , Stepper , vector_steppers< initially_resizer > )
130 {
131     BOOST_STATIC_ASSERT_MSG(
132         ( boost::is_same< typename Stepper::coor_type , diagnostic_state_type >::value ) ,
133         "Coordinate type" );
134     BOOST_STATIC_ASSERT_MSG(
135         ( boost::is_same< typename Stepper::momentum_type , diagnostic_state_type2 >::value ) ,
136         "Momentum type" );
137     BOOST_STATIC_ASSERT_MSG(
138         ( boost::is_same< typename Stepper::coor_deriv_type , diagnostic_deriv_type >::value ) ,
139         "Coordinate deriv type" );
140     BOOST_STATIC_ASSERT_MSG(
141         ( boost::is_same< typename Stepper::momentum_deriv_type , diagnostic_deriv_type2 >::value ) ,
142         "Momentum deriv type" );
143 
144     BOOST_STATIC_ASSERT_MSG(
145         ( boost::is_same< typename Stepper::state_type , std::pair< diagnostic_state_type , diagnostic_state_type2 > >::value ) ,
146         "State type" );
147     BOOST_STATIC_ASSERT_MSG(
148         ( boost::is_same< typename Stepper::deriv_type , std::pair< diagnostic_deriv_type , diagnostic_deriv_type2 > >::value ) ,
149         "Deriv type" );
150 
151     BOOST_STATIC_ASSERT_MSG( ( boost::is_same< typename Stepper::value_type , double >::value ) , "Value type" );
152     BOOST_STATIC_ASSERT_MSG( ( boost::is_same< typename Stepper::time_type , double >::value ) , "Time type" );
153     BOOST_STATIC_ASSERT_MSG( ( boost::is_same< typename Stepper::algebra_type , custom_range_algebra >::value ) , "Algebra type" );
154     BOOST_STATIC_ASSERT_MSG( ( boost::is_same< typename Stepper::operations_type , custom_default_operations >::value ) , "Operations type" );
155 
156     BOOST_STATIC_ASSERT_MSG( ( boost::is_same< typename Stepper::resizer_type , initially_resizer >::value ) , "Resizer type" );
157     BOOST_STATIC_ASSERT_MSG( ( boost::is_same< typename Stepper::stepper_category , stepper_tag >::value ) , "Stepper category" );
158 }
159 
160 
BOOST_AUTO_TEST_CASE_TEMPLATE(test_adjust_size,Stepper,vector_steppers<initially_resizer>)161 BOOST_AUTO_TEST_CASE_TEMPLATE( test_adjust_size , Stepper , vector_steppers< initially_resizer > )
162 {
163     counter_state::init_counter();
164     counter_deriv::init_counter();
165     counter_state2::init_counter();
166     counter_deriv2::init_counter();
167 
168     {
169         Stepper stepper;
170         diagnostic_state_type x;
171         stepper.adjust_size( x );
172     }
173 
174     TEST_COUNTERS( counter_state , 0 , 0 , 0 , 0 );
175     TEST_COUNTERS( counter_state2 , 0 , 0 , 0 , 0 );
176     TEST_COUNTERS( counter_deriv , 1 , 1 , 0 , 1 );
177     TEST_COUNTERS( counter_deriv2 , 1 , 1 , 0 , 1 );
178 }
179 
180 
181 typedef mpl::zip_view< mpl::vector< all_stepper_methods , all_multiplicities > >::type zipped_steppers;
BOOST_AUTO_TEST_CASE_TEMPLATE(test_resizing,Stepper,zipped_steppers)182 BOOST_AUTO_TEST_CASE_TEMPLATE( test_resizing , Stepper , zipped_steppers )
183 {
184     typedef typename mpl::at_c< Stepper , 0 >::type stepper_type;
185     const size_t multiplicity = mpl::at_c< Stepper , 1 >::type::value;
186 
187     counter_state::init_counter();
188     counter_deriv::init_counter();
189     counter_state2::init_counter();
190     counter_deriv2::init_counter();
191 
192     {
193         stepper_type stepper;
194         std::pair< diagnostic_state_type , diagnostic_state_type2 > x;
195         stepper.do_step( constant_mom_func() , x , 0.0 , 0.1 );
196         stepper.do_step( constant_mom_func() , x , 0.0 , 0.1 );
197         stepper.do_step( constant_mom_func() , x , 0.0 , 0.1 );
198     }
199 
200     TEST_COUNTERS( counter_state , 0 , 0 , 0 , 0 );
201     TEST_COUNTERS( counter_state2 , 0 , 0 , 0 , 0 );
202     // dqdt is not needed when called with mom func only, so no resizing
203     // TEST_COUNTERS( counter_deriv , multiplicity , 1 , 0 , 1 );
204     TEST_COUNTERS( counter_deriv2 , multiplicity , 1 , 0 , 1 );
205 }
206 
207 
BOOST_AUTO_TEST_CASE_TEMPLATE(test_copying1,Stepper,vector_steppers<initially_resizer>)208 BOOST_AUTO_TEST_CASE_TEMPLATE( test_copying1 , Stepper , vector_steppers< initially_resizer > )
209 {
210     counter_state::init_counter();
211     counter_deriv::init_counter();
212     counter_state2::init_counter();
213     counter_deriv2::init_counter();
214 
215     {
216         Stepper stepper;
217         Stepper stepper2( stepper );
218     }
219 
220     TEST_COUNTERS( counter_state , 0 , 0 , 0 , 0 );
221     TEST_COUNTERS( counter_state2 , 0 , 0 , 0 , 0 );
222     TEST_COUNTERS( counter_deriv , 0 , 2 , 1 , 2 );
223     TEST_COUNTERS( counter_deriv2 , 0 , 2 , 1 , 2 );
224 }
225 
226 
BOOST_AUTO_TEST_CASE_TEMPLATE(test_copying2,Stepper,vector_steppers<initially_resizer>)227 BOOST_AUTO_TEST_CASE_TEMPLATE( test_copying2 , Stepper , vector_steppers< initially_resizer > )
228 {
229     counter_state::init_counter();
230     counter_deriv::init_counter();
231     counter_state2::init_counter();
232     counter_deriv2::init_counter();
233 
234     {
235         Stepper stepper;
236         std::pair< diagnostic_state_type , diagnostic_state_type2 > x;
237         stepper.do_step( constant_mom_func() , x , 0.0 , 0.1 );
238         Stepper stepper2( stepper );
239     }
240 
241     TEST_COUNTERS( counter_state , 0 , 0 , 0 , 0 );
242     TEST_COUNTERS( counter_state2 , 0 , 0 , 0 , 0 );
243     // dqdt is not needed when called with mom func only, so no resizing
244     //TEST_COUNTERS( counter_deriv , 1 , 2 , 1 , 2 );
245     TEST_COUNTERS( counter_deriv2 , 1 , 2 , 1 , 2 );
246 }
247 
248 
BOOST_AUTO_TEST_CASE_TEMPLATE(test_do_step_v1,Stepper,vector_steppers<initially_resizer>)249 BOOST_AUTO_TEST_CASE_TEMPLATE( test_do_step_v1 , Stepper , vector_steppers< initially_resizer > )
250 {
251     Stepper s;
252     std::pair< diagnostic_state_type , diagnostic_state_type2 > x1 , x2 , x3 , x4;
253     x1.first[0] = 1.0;
254     x1.second[0] = 2.0;
255     x2 = x3 = x4 = x1;
256     diagnostic_state_type x5_coor , x5_mom;
257     x5_coor[0] = x1.first[0];
258     x5_mom[0] = x1.second[0];
259 
260     s.do_step( constant_mom_func() , x1 , 0.0 , 0.1 );
261 
262     s.do_step( std::make_pair( default_coor_func() , constant_mom_func() ) , x2 , 0.0 , 0.1 );
263 
264     default_coor_func cf;
265     constant_mom_func mf;
266     s.do_step( std::make_pair( boost::ref( cf ) , boost::ref( mf ) ) , x3 , 0.0 , 0.1 );
267 
268     std::pair< default_coor_func , constant_mom_func > pf;
269     s.do_step( boost::ref( pf ) , x4 , 0.0 , 0.1 );
270 
271     s.do_step( constant_mom_func() , std::make_pair( boost::ref( x5_coor ) , boost::ref( x5_mom ) ) , 0.0 , 0.1 );
272 
273     // checking for absolute values is not possible here, since the steppers are to different
274     BOOST_CHECK_CLOSE( x1.first[0] , x2.first[0] , 1.0e-14 );
275     BOOST_CHECK_CLOSE( x2.first[0] , x3.first[0] , 1.0e-14 );
276     BOOST_CHECK_CLOSE( x3.first[0] , x4.first[0] , 1.0e-14 );
277     BOOST_CHECK_CLOSE( x4.first[0] , x5_coor[0] , 1.0e-14 );
278 
279     BOOST_CHECK_CLOSE( x1.second[0] , x2.second[0] , 1.0e-14 );
280     BOOST_CHECK_CLOSE( x2.second[0] , x3.second[0] , 1.0e-14 );
281     BOOST_CHECK_CLOSE( x3.second[0] , x4.second[0] , 1.0e-14 );
282     BOOST_CHECK_CLOSE( x4.second[0] , x5_mom[0] , 1.0e-14 );
283 }
284 
BOOST_AUTO_TEST_CASE_TEMPLATE(test_do_step_range,Stepper,vector_steppers<initially_resizer>)285 BOOST_AUTO_TEST_CASE_TEMPLATE( test_do_step_range , Stepper , vector_steppers< initially_resizer > )
286 {
287     Stepper s;
288     diagnostic_state_type q , p ;
289     q[0] = 1.0;
290     p[0] = 2.0;
291 
292     std::vector< double > x;
293     x.push_back( 1.0 );
294     x.push_back( 2.0 );
295     s.do_step( constant_mom_func() ,
296                std::make_pair( x.begin() , x.begin() + 1 ) ,
297                std::make_pair( x.begin() + 1 , x.begin() + 2 ) ,
298                0.0 , 0.1 );
299 
300     s.do_step( constant_mom_func() , q , p , 0.0 , 0.1 );
301 
302     BOOST_CHECK_CLOSE( q[0] , x[0] , 1.0e-14 );
303     BOOST_CHECK_CLOSE( p[0] , x[1] , 1.0e-14 );
304 }
305 
BOOST_AUTO_TEST_CASE_TEMPLATE(test_do_step_v2,Stepper,vector_steppers<initially_resizer>)306 BOOST_AUTO_TEST_CASE_TEMPLATE( test_do_step_v2 , Stepper , vector_steppers< initially_resizer > )
307 {
308     Stepper s;
309     diagnostic_state_type q , p ;
310     q[0] = 1.0;
311     p[0] = 2.0;
312     diagnostic_state_type q2 = q , p2 = p;
313 
314 
315     s.do_step( constant_mom_func() , q , p , 0.0 , 0.1 );
316     s.do_step( constant_mom_func() , std::make_pair( boost::ref( q2 ) , boost::ref( p2 ) ) , 0.0 , 0.1 );
317 
318     BOOST_CHECK_CLOSE( q[0] , q2[0] , 1.0e-14 );
319     BOOST_CHECK_CLOSE( p[0] , p2[0] , 1.0e-14 );
320 }
321 
BOOST_AUTO_TEST_CASE_TEMPLATE(test_do_step_v3,Stepper,vector_steppers<initially_resizer>)322 BOOST_AUTO_TEST_CASE_TEMPLATE( test_do_step_v3 , Stepper , vector_steppers< initially_resizer > )
323 {
324     Stepper s;
325     std::pair< diagnostic_state_type , diagnostic_state_type2 > x_in , x_out;
326     x_in.first[0] = 1.0;
327     x_in.second[0] = 2.0;
328     diagnostic_state_type q2 , p2;
329     q2[0] = x_in.first[0];
330     p2[0] = x_in.second[0];
331 
332 
333     s.do_step( constant_mom_func() , x_in , 0.0 , x_out , 0.1 );
334     s.do_step( constant_mom_func() , std::make_pair( boost::ref( q2 ) , boost::ref( p2 ) ) , 0.0 , 0.1 );
335 
336     BOOST_CHECK_CLOSE( x_in.first[0] , 1.0 , 1.0e-14 );
337     BOOST_CHECK_CLOSE( x_in.second[0] , 2.0 , 1.0e-14 );
338     BOOST_CHECK_CLOSE( x_out.first[0] , q2[0] , 1.0e-14 );
339     BOOST_CHECK_CLOSE( x_out.second[0] , p2[0] , 1.0e-14 );
340 }
341 
342 
343 typedef double vector_space;
344 typedef complete_steppers< vector_space , vector_space , double ,
345                            vector_space , vector_space , double ,
346                            vector_space_algebra , default_operations , initially_resizer > vector_space_steppers;
BOOST_AUTO_TEST_CASE_TEMPLATE(test_with_vector_space_algebra,Stepper,vector_space_steppers)347 BOOST_AUTO_TEST_CASE_TEMPLATE( test_with_vector_space_algebra , Stepper , vector_space_steppers )
348 {
349     Stepper s;
350     std::pair< vector_space , vector_space > x;
351     s.do_step( constant_mom_func_vector_space_1d() , x , 0.0 , 0.1 );
352 
353     s.do_step( std::make_pair( default_coor_func_vector_space_1d() , constant_mom_func_vector_space_1d() ) , x , 0.0 , 0.1 );
354 }
355 
356 
357 typedef boost::fusion::vector< length_type > coor_type;
358 typedef boost::fusion::vector< velocity_type > mom_type;
359 typedef boost::fusion::vector< acceleration_type > acc_type;
360 typedef complete_steppers< coor_type , mom_type , double ,
361                            mom_type , acc_type , time_type,
362                            fusion_algebra , default_operations , initially_resizer > boost_unit_steppers;
BOOST_AUTO_TEST_CASE_TEMPLATE(test_with_boost_units,Stepper,boost_unit_steppers)363 BOOST_AUTO_TEST_CASE_TEMPLATE( test_with_boost_units , Stepper , boost_unit_steppers )
364 {
365     namespace fusion = boost::fusion;
366     namespace si = boost::units::si;
367     Stepper s;
368 
369     coor_type q = fusion::make_vector( 1.0 * si::meter );
370     mom_type p = fusion::make_vector( 2.0 * si::meter_per_second );
371     time_type t = 0.0 * si::second;
372     time_type dt = 0.1 * si::second;
373 
374     coor_type q1 = q , q2 = q;
375     mom_type p1 = p , p2 = p;
376 
377     s.do_step( oscillator_mom_func_units() , std::make_pair( boost::ref( q ) , boost::ref( p ) ) , t , dt );
378 
379     s.do_step( std::make_pair( oscillator_coor_func_units() , oscillator_mom_func_units() ) ,
380                std::make_pair( boost::ref( q1 ) , boost::ref( p1 ) ) , t , dt );
381 
382     s.do_step( oscillator_mom_func_units() , q2 , p2 , t , dt );
383 
384     BOOST_CHECK_CLOSE( ( fusion::at_c< 0 >( q  ).value() ) , ( fusion::at_c< 0 >( q1 ).value() ) , 1.0e-14 );
385     BOOST_CHECK_CLOSE( ( fusion::at_c< 0 >( q1 ).value() ) , ( fusion::at_c< 0 >( q2 ).value() ) , 1.0e-14 );
386     BOOST_CHECK_CLOSE( ( fusion::at_c< 0 >( p  ).value() ) , ( fusion::at_c< 0 >( p1 ).value() ) , 1.0e-14 );
387     BOOST_CHECK_CLOSE( ( fusion::at_c< 0 >( p1 ).value() ) , ( fusion::at_c< 0 >( p2 ).value() ) , 1.0e-14 );
388 }
389 
390 
391 BOOST_AUTO_TEST_SUITE_END()
392