1 /*
2 * Copyright 2012 Karsten Ahnert
3 * Copyright 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 <vector>
12
13 #include <vexcl/vexcl.hpp>
14
15 #include <boost/numeric/odeint.hpp>
16 //[ vexcl_includes
17 #include <boost/numeric/odeint/external/vexcl/vexcl.hpp>
18 //]
19
20 namespace odeint = boost::numeric::odeint;
21
22 //[ vexcl_state_types
23 typedef vex::vector< double > vector_type;
24 typedef vex::multivector< double, 3 > state_type;
25 //]
26
27
28 //[ vexcl_system
29 const double sigma = 10.0;
30 const double b = 8.0 / 3.0;
31
32 struct sys_func
33 {
34 const vector_type &R;
35
sys_funcsys_func36 sys_func( const vector_type &_R ) : R( _R ) { }
37
operator ()sys_func38 void operator()( const state_type &x , state_type &dxdt , double t ) const
39 {
40 dxdt(0) = -sigma * ( x(0) - x(1) );
41 dxdt(1) = R * x(0) - x(1) - x(0) * x(2);
42 dxdt(2) = - b * x(2) + x(0) * x(1);
43 }
44 };
45 //]
46
47
main(int argc,char ** argv)48 int main( int argc , char **argv )
49 {
50 using namespace std;
51 using namespace odeint;
52
53 //[ vexcl_main
54 // setup the opencl context
55 vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
56 std::cout << ctx << std::endl;
57
58 // set up number of system, time step and integration time
59 const size_t n = 1024 * 1024;
60 const double dt = 0.01;
61 const double t_max = 1000.0;
62
63 // initialize R
64 double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
65 std::vector<double> x( n * 3 ) , r( n );
66 for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
67 vector_type R( ctx.queue() , r );
68
69 // initialize the state of the lorenz ensemble
70 state_type X(ctx.queue(), n);
71 X(0) = 10.0;
72 X(1) = 10.0;
73 X(2) = 10.0;
74
75 // create a stepper
76 runge_kutta4< state_type > stepper;
77
78 // solve the system
79 integrate_const( stepper , sys_func( R ) , X , 0.0 , t_max , dt );
80 //]
81
82 std::vector< double > res( 3 * n );
83 vex::copy( X(0) , res );
84 for( size_t i=0 ; i<n ; ++i )
85 cout << r[i] << "\t" << res[i] << "\t" << "\n";
86 }
87