1 /*
2 * phase_chain.cpp
3 *
4 * Example of OMP parallelization with odeint
5 *
6 * Copyright 2013 Karsten Ahnert
7 * Copyright 2013 Mario Mulansky
8 * Copyright 2013 Pascal Germroth
9 * Distributed under the Boost Software License, Version 1.0. (See
10 * accompanying file LICENSE_1_0.txt or copy at
11 * http://www.boost.org/LICENSE_1_0.txt)
12 */
13
14 #include <iostream>
15 #include <vector>
16 #include <boost/random.hpp>
17 #include <boost/timer/timer.hpp>
18 //[phase_chain_openmp_header
19 #include <omp.h>
20 #include <boost/numeric/odeint.hpp>
21 #include <boost/numeric/odeint/external/openmp/openmp.hpp>
22 //]
23
24 using namespace std;
25 using namespace boost::numeric::odeint;
26 using boost::timer::cpu_timer;
27 using boost::math::double_constants::pi;
28
29 //[phase_chain_vector_state
30 typedef std::vector< double > state_type;
31 //]
32
33 //[phase_chain_rhs
34 struct phase_chain
35 {
phase_chainphase_chain36 phase_chain( double gamma = 0.5 )
37 : m_gamma( gamma ) { }
38
operator ()phase_chain39 void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
40 {
41 const size_t N = x.size();
42 #pragma omp parallel for schedule(runtime)
43 for(size_t i = 1 ; i < N - 1 ; ++i)
44 {
45 dxdt[i] = coupling_func( x[i+1] - x[i] ) +
46 coupling_func( x[i-1] - x[i] );
47 }
48 dxdt[0 ] = coupling_func( x[1 ] - x[0 ] );
49 dxdt[N-1] = coupling_func( x[N-2] - x[N-1] );
50 }
51
coupling_funcphase_chain52 double coupling_func( double x ) const
53 {
54 return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
55 }
56
57 double m_gamma;
58 };
59 //]
60
61
main(int argc,char ** argv)62 int main( int argc , char **argv )
63 {
64 //[phase_chain_init
65 size_t N = 131101;
66 state_type x( N );
67 boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
68 boost::random::mt19937 engine( 0 );
69 generate( x.begin() , x.end() , boost::bind( distribution , engine ) );
70 //]
71
72 //[phase_chain_stepper
73 typedef runge_kutta4<
74 state_type , double ,
75 state_type , double ,
76 openmp_range_algebra
77 > stepper_type;
78 //]
79
80 //[phase_chain_scheduling
81 int chunk_size = N/omp_get_max_threads();
82 omp_set_schedule( omp_sched_static , chunk_size );
83 //]
84
85 cpu_timer timer;
86 //[phase_chain_integrate
87 integrate_n_steps( stepper_type() , phase_chain( 1.2 ) ,
88 x , 0.0 , 0.01 , 100 );
89 //]
90 double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
91 std::cerr << run_time << "s" << std::endl;
92 // copy(x.begin(), x.end(), ostream_iterator<double>(cout, "\n"));
93
94 return 0;
95 }
96