1 /*
2 * phase_chain.cpp
3 *
4 * Example of MPI 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_mpi_header
19 #include <boost/numeric/odeint.hpp>
20 #include <boost/numeric/odeint/external/mpi/mpi.hpp>
21 //]
22
23 using namespace std;
24 using namespace boost::numeric::odeint;
25 using boost::timer::cpu_timer;
26 using boost::math::double_constants::pi;
27
28 //[phase_chain_state
29 typedef mpi_state< vector<double> > state_type;
30 //]
31
32 //[phase_chain_mpi_rhs
33 struct phase_chain_mpi_state
34 {
phase_chain_mpi_statephase_chain_mpi_state35 phase_chain_mpi_state( double gamma = 0.5 )
36 : m_gamma( gamma ) { }
37
operator ()phase_chain_mpi_state38 void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
39 {
40 const size_t M = x().size();
41 const bool have_left = x.world.rank() > 0,
42 have_right = x.world.rank() < x.world.size()-1;
43 double x_left, x_right;
44 boost::mpi::request r_left, r_right;
45 if( have_left )
46 {
47 x.world.isend( x.world.rank()-1 , 0 , x().front() ); // send to x_right
48 r_left = x.world.irecv( x.world.rank()-1 , 0 , x_left ); // receive from x().back()
49 }
50 if( have_right )
51 {
52 x.world.isend( x.world.rank()+1 , 0 , x().back() ); // send to x_left
53 r_right = x.world.irecv( x.world.rank()+1 , 0 , x_right ); // receive from x().front()
54 }
55 for(size_t m = 1 ; m < M-1 ; ++m)
56 {
57 dxdt()[m] = coupling_func( x()[m+1] - x()[m] ) +
58 coupling_func( x()[m-1] - x()[m] );
59 }
60 dxdt()[0] = coupling_func( x()[1] - x()[0] );
61 if( have_left )
62 {
63 r_left.wait();
64 dxdt()[0] += coupling_func( x_left - x().front() );
65 }
66 dxdt()[M-1] = coupling_func( x()[M-2] - x()[M-1] );
67 if( have_right )
68 {
69 r_right.wait();
70 dxdt()[M-1] += coupling_func( x_right - x().back() );
71 }
72 }
73
coupling_funcphase_chain_mpi_state74 double coupling_func( double x ) const
75 {
76 return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
77 }
78
79 double m_gamma;
80 };
81 //]
82
83
main(int argc,char ** argv)84 int main( int argc , char **argv )
85 {
86 //[phase_chain_mpi_init
87 boost::mpi::environment env( argc , argv );
88 boost::mpi::communicator world;
89
90 const size_t N = 131101;
91 vector<double> x;
92 if( world.rank() == 0 )
93 {
94 x.resize( N );
95 boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
96 boost::random::mt19937 engine( 0 );
97 generate( x.begin() , x.end() , boost::bind( distribution , engine ) );
98 }
99
100 state_type x_split( world );
101 split( x , x_split );
102 //]
103
104
105 cpu_timer timer;
106 //[phase_chain_mpi_integrate
107 integrate_n_steps( runge_kutta4<state_type>() , phase_chain_mpi_state( 1.2 ) ,
108 x_split , 0.0 , 0.01 , 100 );
109 unsplit( x_split , x );
110 //]
111
112 if( world.rank() == 0 )
113 {
114 double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
115 std::cerr << run_time << "s" << std::endl;
116 // copy(x.begin(), x.end(), ostream_iterator<double>(cout, "\n"));
117 }
118
119 return 0;
120 }
121