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

   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__
  68         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
  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
  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
 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 }