1 /* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble_simple.cpp
2
3 Copyright 2013 Karsten Ahnert
4 Copyright 2013 Mario Mulansky
5 Copyright 2013 Pascal Germroth
6
7 Parallelized Lorenz ensembles
8
9 Distributed under the Boost Software License, Version 1.0.
10 (See accompanying file LICENSE_1_0.txt or
11 copy at http://www.boost.org/LICENSE_1_0.txt)
12 */
13
14 #include <omp.h>
15 #include <vector>
16 #include <iostream>
17 #include <boost/numeric/odeint.hpp>
18 #include <boost/numeric/odeint/external/openmp/openmp.hpp>
19 #include "point_type.hpp"
20
21 using namespace std;
22
23 typedef vector<double> vector_type;
24 typedef point<double, 3> point_type;
25 typedef vector<point_type> state_type;
26
27 const double sigma = 10.0;
28 const double b = 8.0 / 3.0;
29
30 struct sys_func {
31 const vector_type &R;
sys_funcsys_func32 sys_func( const vector_type &_R ) : R( _R ) { }
33
operator ()sys_func34 void operator()( const state_type &x , state_type &dxdt , double t ) const {
35 const size_t n = x.size();
36 # pragma omp parallel for
37 for(size_t i = 0 ; i < n ; i++) {
38 const point_type &xi = x[i];
39 point_type &dxdti = dxdt[i];
40 dxdti[0] = -sigma * (xi[0] - xi[1]);
41 dxdti[1] = R[i] * xi[0] - xi[1] - xi[0] * xi[2];
42 dxdti[2] = -b * xi[2] + xi[0] * xi[1];
43 }
44 }
45 };
46
47
main()48 int main() {
49 using namespace boost::numeric::odeint;
50
51 const size_t n = 1024;
52 vector_type R(n);
53 const double Rmin = 0.1, Rmax = 50.0;
54 # pragma omp parallel for
55 for(size_t i = 0 ; i < n ; i++)
56 R[i] = Rmin + (Rmax - Rmin) / (n - 1) * i;
57
58 state_type X(n, point_type(10, 10, 10));
59
60 const double t_max = 10.0, dt = 0.01;
61
62 // Simple stepper with constant step size
63 // typedef runge_kutta4<state_type, double, state_type, double,
64 // openmp_range_algebra> stepper;
65
66 // integrate_const(stepper(), sys_func(R), X, 0.0, t_max, dt);
67
68 // Controlled stepper with variable step size
69 typedef runge_kutta_fehlberg78<state_type, double, state_type, double,
70 openmp_range_algebra> error_stepper_type;
71 typedef controlled_runge_kutta<error_stepper_type> controlled_stepper_type;
72 controlled_stepper_type controlled_stepper;
73
74 integrate_adaptive(controlled_stepper, sys_func(R), X, 0.0, t_max, dt);
75
76 copy( X.begin(), X.end(), ostream_iterator<point_type>(cout, "\n") );
77
78 return 0;
79 }
80