1 /*
2 Copyright 2011-2012 Karsten Ahnert
3 Copyright 2011-2013 Mario Mulansky
4
5 Distributed under the Boost Software License, Version 1.0.
6 (See accompanying file LICENSE_1_0.txt or
7 copy at http://www.boost.org/LICENSE_1_0.txt)
8 */
9
10 #include <iostream>
11 #include <cmath>
12 #include <utility>
13
14
15 #include <thrust/device_vector.h>
16 #include <thrust/reduce.h>
17 #include <thrust/functional.h>
18
19 #include <boost/numeric/odeint.hpp>
20
21 #include <boost/numeric/odeint/external/thrust/thrust.hpp>
22
23 #include <boost/random/mersenne_twister.hpp>
24 #include <boost/random/uniform_real.hpp>
25 #include <boost/random/variate_generator.hpp>
26
27
28 using namespace std;
29 using namespace boost::numeric::odeint;
30
31 //change this to float if your device does not support double computation
32 typedef double value_type;
33
34 //change this to host_vector< ... > of you want to run on CPU
35 typedef thrust::device_vector< value_type > state_type;
36 typedef thrust::device_vector< size_t > index_vector_type;
37 // typedef thrust::host_vector< value_type > state_type;
38 // typedef thrust::host_vector< size_t > index_vector_type;
39
40
41 const value_type sigma = 10.0;
42 const value_type b = 8.0 / 3.0;
43
44
45 //[ thrust_lorenz_parameters_define_simple_system
46 struct lorenz_system
47 {
48 struct lorenz_functor
49 {
50 template< class T >
51 __host__ __device__
operator ()lorenz_system::lorenz_functor52 void operator()( T t ) const
53 {
54 // unpack the parameter we want to vary and the Lorenz variables
55 value_type R = thrust::get< 3 >( t );
56 value_type x = thrust::get< 0 >( t );
57 value_type y = thrust::get< 1 >( t );
58 value_type z = thrust::get< 2 >( t );
59 thrust::get< 4 >( t ) = sigma * ( y - x );
60 thrust::get< 5 >( t ) = R * x - y - x * z;
61 thrust::get< 6 >( t ) = -b * z + x * y ;
62
63 }
64 };
65
lorenz_systemlorenz_system66 lorenz_system( size_t N , const state_type &beta )
67 : m_N( N ) , m_beta( beta ) { }
68
69 template< class State , class Deriv >
operator ()lorenz_system70 void operator()( const State &x , Deriv &dxdt , value_type t ) const
71 {
72 thrust::for_each(
73 thrust::make_zip_iterator( thrust::make_tuple(
74 boost::begin( x ) ,
75 boost::begin( x ) + m_N ,
76 boost::begin( x ) + 2 * m_N ,
77 m_beta.begin() ,
78 boost::begin( dxdt ) ,
79 boost::begin( dxdt ) + m_N ,
80 boost::begin( dxdt ) + 2 * m_N ) ) ,
81 thrust::make_zip_iterator( thrust::make_tuple(
82 boost::begin( x ) + m_N ,
83 boost::begin( x ) + 2 * m_N ,
84 boost::begin( x ) + 3 * m_N ,
85 m_beta.begin() ,
86 boost::begin( dxdt ) + m_N ,
87 boost::begin( dxdt ) + 2 * m_N ,
88 boost::begin( dxdt ) + 3 * m_N ) ) ,
89 lorenz_functor() );
90 }
91 size_t m_N;
92 const state_type &m_beta;
93 };
94 //]
95
96 struct lorenz_perturbation_system
97 {
98 struct lorenz_perturbation_functor
99 {
100 template< class T >
101 __host__ __device__
operator ()lorenz_perturbation_system::lorenz_perturbation_functor102 void operator()( T t ) const
103 {
104 value_type R = thrust::get< 1 >( t );
105 value_type x = thrust::get< 0 >( thrust::get< 0 >( t ) );
106 value_type y = thrust::get< 1 >( thrust::get< 0 >( t ) );
107 value_type z = thrust::get< 2 >( thrust::get< 0 >( t ) );
108 value_type dx = thrust::get< 3 >( thrust::get< 0 >( t ) );
109 value_type dy = thrust::get< 4 >( thrust::get< 0 >( t ) );
110 value_type dz = thrust::get< 5 >( thrust::get< 0 >( t ) );
111 thrust::get< 0 >( thrust::get< 2 >( t ) ) = sigma * ( y - x );
112 thrust::get< 1 >( thrust::get< 2 >( t ) ) = R * x - y - x * z;
113 thrust::get< 2 >( thrust::get< 2 >( t ) ) = -b * z + x * y ;
114 thrust::get< 3 >( thrust::get< 2 >( t ) ) = sigma * ( dy - dx );
115 thrust::get< 4 >( thrust::get< 2 >( t ) ) = ( R - z ) * dx - dy - x * dz;
116 thrust::get< 5 >( thrust::get< 2 >( t ) ) = y * dx + x * dy - b * dz;
117 }
118 };
119
lorenz_perturbation_systemlorenz_perturbation_system120 lorenz_perturbation_system( size_t N , const state_type &beta )
121 : m_N( N ) , m_beta( beta ) { }
122
123 template< class State , class Deriv >
operator ()lorenz_perturbation_system124 void operator()( const State &x , Deriv &dxdt , value_type t ) const
125 {
126 thrust::for_each(
127 thrust::make_zip_iterator( thrust::make_tuple(
128 thrust::make_zip_iterator( thrust::make_tuple(
129 boost::begin( x ) ,
130 boost::begin( x ) + m_N ,
131 boost::begin( x ) + 2 * m_N ,
132 boost::begin( x ) + 3 * m_N ,
133 boost::begin( x ) + 4 * m_N ,
134 boost::begin( x ) + 5 * m_N ) ) ,
135 m_beta.begin() ,
136 thrust::make_zip_iterator( thrust::make_tuple(
137 boost::begin( dxdt ) ,
138 boost::begin( dxdt ) + m_N ,
139 boost::begin( dxdt ) + 2 * m_N ,
140 boost::begin( dxdt ) + 3 * m_N ,
141 boost::begin( dxdt ) + 4 * m_N ,
142 boost::begin( dxdt ) + 5 * m_N ) )
143 ) ) ,
144 thrust::make_zip_iterator( thrust::make_tuple(
145 thrust::make_zip_iterator( thrust::make_tuple(
146 boost::begin( x ) + m_N ,
147 boost::begin( x ) + 2 * m_N ,
148 boost::begin( x ) + 3 * m_N ,
149 boost::begin( x ) + 4 * m_N ,
150 boost::begin( x ) + 5 * m_N ,
151 boost::begin( x ) + 6 * m_N ) ) ,
152 m_beta.begin() ,
153 thrust::make_zip_iterator( thrust::make_tuple(
154 boost::begin( dxdt ) + m_N ,
155 boost::begin( dxdt ) + 2 * m_N ,
156 boost::begin( dxdt ) + 3 * m_N ,
157 boost::begin( dxdt ) + 4 * m_N ,
158 boost::begin( dxdt ) + 5 * m_N ,
159 boost::begin( dxdt ) + 6 * m_N ) )
160 ) ) ,
161 lorenz_perturbation_functor() );
162 }
163
164 size_t m_N;
165 const state_type &m_beta;
166 };
167
168 struct lyap_observer
169 {
170 //[thrust_lorenz_parameters_observer_functor
171 struct lyap_functor
172 {
173 template< class T >
174 __host__ __device__
operator ()lyap_observer::lyap_functor175 void operator()( T t ) const
176 {
177 value_type &dx = thrust::get< 0 >( t );
178 value_type &dy = thrust::get< 1 >( t );
179 value_type &dz = thrust::get< 2 >( t );
180 value_type norm = sqrt( dx * dx + dy * dy + dz * dz );
181 dx /= norm;
182 dy /= norm;
183 dz /= norm;
184 thrust::get< 3 >( t ) += log( norm );
185 }
186 };
187 //]
188
lyap_observerlyap_observer189 lyap_observer( size_t N , size_t every = 100 )
190 : m_N( N ) , m_lyap( N ) , m_every( every ) , m_count( 0 )
191 {
192 thrust::fill( m_lyap.begin() , m_lyap.end() , 0.0 );
193 }
194
195 template< class Lyap >
fill_lyaplyap_observer196 void fill_lyap( Lyap &lyap )
197 {
198 thrust::copy( m_lyap.begin() , m_lyap.end() , lyap.begin() );
199 for( size_t i=0 ; i<lyap.size() ; ++i )
200 lyap[i] /= m_t_overall;
201 }
202
203
204 template< class State >
operator ()lyap_observer205 void operator()( State &x , value_type t )
206 {
207 if( ( m_count != 0 ) && ( ( m_count % m_every ) == 0 ) )
208 {
209 thrust::for_each(
210 thrust::make_zip_iterator( thrust::make_tuple(
211 boost::begin( x ) + 3 * m_N ,
212 boost::begin( x ) + 4 * m_N ,
213 boost::begin( x ) + 5 * m_N ,
214 m_lyap.begin() ) ) ,
215 thrust::make_zip_iterator( thrust::make_tuple(
216 boost::begin( x ) + 4 * m_N ,
217 boost::begin( x ) + 5 * m_N ,
218 boost::begin( x ) + 6 * m_N ,
219 m_lyap.end() ) ) ,
220 lyap_functor() );
221 clog << t << "\n";
222 }
223 ++m_count;
224 m_t_overall = t;
225 }
226
227 size_t m_N;
228 state_type m_lyap;
229 size_t m_every;
230 size_t m_count;
231 value_type m_t_overall;
232 };
233
234 const size_t N = 1024*2;
235 const value_type dt = 0.01;
236
237
main(int arc,char * argv[])238 int main( int arc , char* argv[] )
239 {
240 int driver_version , runtime_version;
241 cudaDriverGetVersion( &driver_version );
242 cudaRuntimeGetVersion ( &runtime_version );
243 cout << driver_version << "\t" << runtime_version << endl;
244
245
246 //[ thrust_lorenz_parameters_define_beta
247 vector< value_type > beta_host( N );
248 const value_type beta_min = 0.0 , beta_max = 56.0;
249 for( size_t i=0 ; i<N ; ++i )
250 beta_host[i] = beta_min + value_type( i ) * ( beta_max - beta_min ) / value_type( N - 1 );
251
252 state_type beta = beta_host;
253 //]
254
255 //[ thrust_lorenz_parameters_integration
256 state_type x( 6 * N );
257
258 // initialize x,y,z
259 thrust::fill( x.begin() , x.begin() + 3 * N , 10.0 );
260
261 // initial dx
262 thrust::fill( x.begin() + 3 * N , x.begin() + 4 * N , 1.0 );
263
264 // initialize dy,dz
265 thrust::fill( x.begin() + 4 * N , x.end() , 0.0 );
266
267
268 // create error stepper, can be used with make_controlled or make_dense_output
269 typedef runge_kutta_dopri5< state_type , value_type , state_type , value_type > stepper_type;
270
271
272 lorenz_system lorenz( N , beta );
273 lorenz_perturbation_system lorenz_perturbation( N , beta );
274 lyap_observer obs( N , 1 );
275
276 // calculate transients
277 integrate_adaptive( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , lorenz , std::make_pair( x.begin() , x.begin() + 3 * N ) , 0.0 , 10.0 , dt );
278
279 // calculate the Lyapunov exponents -- the main loop
280 double t = 0.0;
281 while( t < 10000.0 )
282 {
283 integrate_adaptive( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , lorenz_perturbation , x , t , t + 1.0 , 0.1 );
284 t += 1.0;
285 obs( x , t );
286 }
287
288 vector< value_type > lyap( N );
289 obs.fill_lyap( lyap );
290
291 for( size_t i=0 ; i<N ; ++i )
292 cout << beta_host[i] << "\t" << lyap[i] << "\n";
293 //]
294
295 return 0;
296 }
297