1 /* Boost lorenz.cpp test file
2
3 Copyright 2012 Karsten Ahnert
4 Copyright 2012 Mario Mulansky
5
6 This file tests the odeint library with the vexcl types
7
8 Distributed under the Boost Software License, Version 1.0.
9 (See accompanying file LICENSE_1_0.txt or
10 copy at http://www.boost.org/LICENSE_1_0.txt)
11 */
12
13 #define BOOST_TEST_MODULE odeint_vexcl
14
15 #include <vector>
16 #include <iostream>
17
18 #include <vexcl/vexcl.hpp>
19
20 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
21 #include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
22 #include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
23 #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
24 #include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
25 #include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
26 #include <boost/numeric/odeint/integrate/integrate_const.hpp>
27 #include <boost/numeric/odeint/external/vexcl/vexcl.hpp>
28
29 #include <boost/test/unit_test.hpp>
30
31 using namespace boost::unit_test;
32 namespace odeint = boost::numeric::odeint;
33
34
35 typedef vex::vector< double > vector_type;
36 typedef vex::multivector< double, 3 > state_type;
37
38
39
40 const double sigma = 10.0;
41 const double b = 8.0 / 3.0;
42
43 struct sys_func
44 {
45 const vector_type &R;
46
sys_funcsys_func47 sys_func( const vector_type &_R ) : R( _R ) { }
48
operator ()sys_func49 void operator()( const state_type &x , state_type &dxdt , double t ) const
50 {
51 dxdt(0) = -sigma * ( x(0) - x(1) );
52 dxdt(1) = R * x(0) - x(1) - x(0) * x(2);
53 dxdt(2) = - b * x(2) + x(0) * x(1);
54 }
55 };
56
57 const size_t n = 1024 ;
58 const double dt = 0.01;
59 const double t_max = 1.0;
60
61
BOOST_AUTO_TEST_CASE(integrate_const_rk4)62 BOOST_AUTO_TEST_CASE( integrate_const_rk4 )
63 {
64 vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
65
66 double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
67 std::vector<double> x( n * 3 ) , r( n );
68 for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
69
70 state_type X( ctx.queue() , n );
71 X(0) = 10.0;
72 X(1) = 10.0;
73 X(2) = 10.0;
74
75 vector_type R( ctx.queue() , r );
76
77 odeint::runge_kutta4<
78 state_type , double , state_type , double ,
79 odeint::vector_space_algebra , odeint::default_operations
80 > stepper;
81
82 odeint::integrate_const( stepper , sys_func( R ) , X , 0.0 , t_max , dt );
83
84 std::vector< double > res( 3 * n );
85 vex::copy( X(0) , res );
86 std::cout << res[0] << std::endl;
87 }
88
BOOST_AUTO_TEST_CASE(integrate_const_controlled_rk54)89 BOOST_AUTO_TEST_CASE( integrate_const_controlled_rk54 )
90 {
91 vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
92
93 double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
94 std::vector<double> x( n * 3 ) , r( n );
95 for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
96
97 state_type X( ctx.queue() , n );
98 X(0) = 10.0;
99 X(1) = 10.0;
100 X(2) = 10.0;
101
102 vector_type R( ctx.queue() , r );
103
104 typedef odeint::runge_kutta_cash_karp54<
105 state_type , double , state_type , double ,
106 odeint::vector_space_algebra , odeint::default_operations
107 > stepper_type;
108 typedef odeint::controlled_runge_kutta< stepper_type > controlled_stepper_type;
109
110 odeint::integrate_const( controlled_stepper_type() , sys_func( R ) , X , 0.0 , t_max , dt );
111
112 std::vector< double > res( 3 * n );
113 vex::copy( X(0) , res );
114 std::cout << res[0] << std::endl;
115 }
116
BOOST_AUTO_TEST_CASE(integrate_const_dense_output_dopri5)117 BOOST_AUTO_TEST_CASE( integrate_const_dense_output_dopri5 )
118 {
119 vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
120
121 double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
122 std::vector<double> x( n * 3 ) , r( n );
123 for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
124
125 state_type X( ctx.queue() , n );
126 X(0) = 10.0;
127 X(1) = 10.0;
128 X(2) = 10.0;
129
130 vector_type R( ctx.queue() , r );
131
132 typedef odeint::runge_kutta_dopri5<
133 state_type , double , state_type , double ,
134 odeint::vector_space_algebra , odeint::default_operations
135 > stepper_type;
136 typedef odeint::controlled_runge_kutta< stepper_type > controlled_stepper_type;
137 typedef odeint::dense_output_runge_kutta< controlled_stepper_type > dense_output_stepper_type;
138
139 odeint::integrate_const( dense_output_stepper_type() , sys_func( R ) , X , 0.0 , t_max , dt );
140
141 std::vector< double > res( 3 * n );
142 vex::copy( X(0) , res );
143 std::cout << res[0] << std::endl;
144 }
145
146
147
148