1 /*
2  Copyright 2011-2013 Mario Mulansky
3  Copyright 2011 Karsten Ahnert
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 /*
11  * This example shows how to use odeint on CUDA devices with thrust.
12  * Note that we require at least Version 3.2 of the nVidia CUDA SDK
13  * and the thrust library should be installed in the CUDA include
14  * folder.
15  *
16  * As example we use a chain of phase oscillators with nearest neighbour
17  * coupling, as described in:
18  *
19  * Avis H. Cohen, Philip J. Holmes and Richard H. Rand:
20  * JOURNAL OF MATHEMATICAL BIOLOGY Volume 13, Number 3, 345-369,
21  *
22  */
23 
24 #include <iostream>
25 #include <cmath>
26 
27 #include <thrust/device_vector.h>
28 #include <thrust/iterator/permutation_iterator.h>
29 #include <thrust/iterator/counting_iterator.h>
30 
31 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
32 #include <boost/numeric/odeint/integrate/integrate_const.hpp>
33 #include <boost/numeric/odeint/external/thrust/thrust.hpp>
34 
35 using namespace std;
36 
37 using namespace boost::numeric::odeint;
38 
39 
40 //change this to float if your device does not support double computation
41 typedef double value_type;
42 
43 
44 //[ thrust_phase_chain_system
45 //change this to host_vector< ... > if you want to run on CPU
46 typedef thrust::device_vector< value_type > state_type;
47 typedef thrust::device_vector< size_t > index_vector_type;
48 //typedef thrust::host_vector< value_type > state_type;
49 //typedef thrust::host_vector< size_t > index_vector_type;
50 
51 //<-
52 /*
53  * This implements the rhs of the dynamical equation:
54  * \phi'_0 = \omega_0 + sin( \phi_1 - \phi_0 )
55  * \phi'_i  = \omega_i + sin( \phi_i+1 - \phi_i ) + sin( \phi_i - \phi_i-1 )
56  * \phi'_N-1 = \omega_N-1 + sin( \phi_N-1 - \phi_N-2 )
57  */
58 //->
59 class phase_oscillators
60 {
61 
62 public:
63 
64     struct sys_functor
65     {
66         template< class Tuple >
67         __host__ __device__
operator ()phase_oscillators::sys_functor68         void operator()( Tuple t )  // this functor works on tuples of values
69         {
70             // first, unpack the tuple into value, neighbors and omega
71             const value_type phi = thrust::get<0>(t);
72             const value_type phi_left = thrust::get<1>(t);  // left neighbor
73             const value_type phi_right = thrust::get<2>(t); // right neighbor
74             const value_type omega = thrust::get<3>(t);
75             // the dynamical equation
76             thrust::get<4>(t) = omega + sin( phi_right - phi ) + sin( phi - phi_left );
77         }
78     };
79 
phase_oscillators(const state_type & omega)80     phase_oscillators( const state_type &omega )
81         : m_omega( omega ) , m_N( omega.size() ) , m_prev( omega.size() ) , m_next( omega.size() )
82     {
83         // build indices pointing to left and right neighbours
84         thrust::counting_iterator<size_t> c( 0 );
85         thrust::copy( c , c+m_N-1 , m_prev.begin()+1 );
86         m_prev[0] = 0; // m_prev = { 0 , 0 , 1 , 2 , 3 , ... , N-1 }
87 
88         thrust::copy( c+1 , c+m_N , m_next.begin() );
89         m_next[m_N-1] = m_N-1; // m_next = { 1 , 2 , 3 , ... , N-1 , N-1 }
90     }
91 
operator ()(const state_type & x,state_type & dxdt,const value_type dt)92     void operator() ( const state_type &x , state_type &dxdt , const value_type dt )
93     {
94         thrust::for_each(
95                 thrust::make_zip_iterator(
96                         thrust::make_tuple(
97                                 x.begin() ,
98                                 thrust::make_permutation_iterator( x.begin() , m_prev.begin() ) ,
99                                 thrust::make_permutation_iterator( x.begin() , m_next.begin() ) ,
100                                 m_omega.begin() ,
101                                 dxdt.begin()
102                                 ) ),
103                 thrust::make_zip_iterator(
104                         thrust::make_tuple(
105                                 x.end() ,
106                                 thrust::make_permutation_iterator( x.begin() , m_prev.end() ) ,
107                                 thrust::make_permutation_iterator( x.begin() , m_next.end() ) ,
108                                 m_omega.end() ,
109                                 dxdt.end()) ) ,
110                 sys_functor()
111                 );
112     }
113 
114 private:
115 
116     const state_type &m_omega;
117     const size_t m_N;
118     index_vector_type m_prev;
119     index_vector_type m_next;
120 };
121 //]
122 
123 const size_t N = 32768;
124 const value_type pi = 3.1415926535897932384626433832795029;
125 const value_type epsilon = 6.0 / ( N * N ); // should be < 8/N^2 to see phase locking
126 const value_type dt = 0.1;
127 
main(int arc,char * argv[])128 int main( int arc , char* argv[] )
129 {
130     //[ thrust_phase_chain_integration
131     // create initial conditions and omegas on host:
132     vector< value_type > x_host( N );
133     vector< value_type > omega_host( N );
134     for( size_t i=0 ; i<N ; ++i )
135     {
136         x_host[i] = 2.0 * pi * drand48();
137         omega_host[i] = ( N - i ) * epsilon; // decreasing frequencies
138     }
139 
140     // copy to device
141     state_type x = x_host;
142     state_type omega = omega_host;
143 
144     // create stepper
145     runge_kutta4< state_type , value_type , state_type , value_type > stepper;
146 
147     // create phase oscillator system function
148     phase_oscillators sys( omega );
149 
150     // integrate
151     integrate_const( stepper , sys , x , 0.0 , 10.0 , dt );
152 
153     thrust::copy( x.begin() , x.end() , std::ostream_iterator< value_type >( std::cout , "\n" ) );
154     std::cout << std::endl;
155     //]
156 }
157