libs/numeric/odeint/examples/openmp/lorenz_ensemble_nested.cpp

   1 /* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble_nested.cpp
   2
   3  Copyright 2013 Karsten Ahnert
   4  Copyright 2013 Pascal Germroth
   5  Copyright 2013 Mario Mulansky
   6
   7  Parallelized Lorenz ensembles using nested omp algebra
   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 <iterator>
  18 #include <boost/numeric/odeint.hpp>
  19 #include <boost/numeric/odeint/external/openmp/openmp.hpp>
  20 #include <boost/lexical_cast.hpp>
  21 #include "point_type.hpp"
  22
  23 using namespace std;
  24 using namespace boost::numeric::odeint;
  25
  26 typedef point<double, 3> point_type;
  27 typedef vector< point_type > state_type;
  28
  29 const double sigma = 10.0;
  30 const double b = 8.0 / 3.0;
  31
  32
  33 struct sys_func {
  34     const vector<double> &R;
  35     sys_func( vector<double> &R ) : R(R) {}
  36
  37     void operator()( const state_type &x , state_type &dxdt , double t ) const {
  38 #       pragma omp parallel for
  39         for(size_t i = 0 ; i < x.size() ; i++) {
  40             dxdt[i][0] = -sigma * (x[i][0] - x[i][1]);
  41             dxdt[i][1] = R[i] * x[i][0] - x[i][1] - x[i][0] * x[i][2];
  42             dxdt[i][2] = -b * x[i][2] + x[i][0] * x[i][1];
  43         }
  44     }
  45 };
  46
  47
  48 int main(int argc, char **argv) {
  49     size_t n = 1024;
  50     if(argc > 1) n = boost::lexical_cast<size_t>(argv[1]);
  51
  52     vector<double> 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 state( n , point_type(10, 10, 10) );
  59
  60     typedef runge_kutta4< state_type, double , state_type , double ,
  61                           openmp_nested_algebra<vector_space_algebra> > stepper;
  62
  63     const double t_max = 10.0, dt = 0.01;
  64
  65     integrate_const(
  66         stepper(),
  67         sys_func(R),
  68         state,
  69         0.0, t_max, dt
  70     );
  71
  72     std::copy( state.begin(), state.end(), ostream_iterator<point_type>(cout, "\n") );
  73
  74     return 0;
  75 }