ceph/src/boost/libs/numeric/odeint/examples/thrust/lorenz_parameters.cu

   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__
  52         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
  66     lorenz_system( size_t N , const state_type &beta )
  67     : m_N( N ) , m_beta( beta ) { }
  68
  69     template< class State , class Deriv >
  70     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__
 102         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
 120     lorenz_perturbation_system( size_t N , const state_type &beta )
 121     : m_N( N ) , m_beta( beta ) { }
 122
 123     template< class State , class Deriv >
 124     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__
 175         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
 189     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 >
 196     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 >
 205     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
 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 }