[Odeint] Utiliser odeint avec openmp et un integrateur adaptatif
(je m'excuse par avance, ecrivant sur un clavier anglais je n'ai pas acces aux accents)
Bonjour,
j'essaie d'utiliser odeint avec openmp (pour paralleliser le calcul). Tout marche tres bien lorsque j'utilise des pas constants, mais j'obtiens lorsque je veux utiliser des pas adaptatifs des erreurs que je ne parviens pas a comprendre.
Voila les parties essentielles du code :
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
|
using namespace std;
using namespace boost::numeric::odeint;
const int jmax = 10;
typedef double value_type;
typedef boost::array<value_type ,2*(jmax+1) > state_type;
//The step function
void rhs( const state_type A , state_type &dAdt , const value_type t )
{
value_type RHStemp0;
value_type RHStemp1;
//We will write the RHS of the equations a big sum
#pragma omp parallel for schedule(runtime)
for(int j = 0; j < jmax+1 ; j++ ) //Real part
{
RHStemp0 = value_type(0.0);
RHStemp1 = value_type(0.0);
for (int k = 0; k< jmax+1 ;k++)
{
for (int l = max(0,j+k-jmax); l < 1 + min(jmax,j+k);l++)
{
RHStemp0 = RHStemp0 + S[j*SIZE_S*SIZE_S + k*SIZE_S + l]*(-A[k+jmax+1]*A[l]*A[j+k-l] + A[k]*A[l+jmax+1]*A[j+k-l]
+ A[k]*A[l]*A[j+k-l+jmax+1] +A[k+jmax+1]*A[l+jmax+1]*A[j+k-l+jmax+1]);
RHStemp1 = RHStemp1 + S[j*SIZE_S*SIZE_S + k*SIZE_S + l]*(A[k]*A[l]*A[j+k-l] - A[k]*A[l+jmax+1]*A[j+k-l+jmax+1]
+ A[k+jmax+1]*A[l]*A[j+k-l+jmax+1] +A[k+jmax+1]*A[l+jmax+1]*A[j+k-l]);
}
}
dAdt[j] = (-1/(value_type((2*(2*j+3)))))*RHStemp0;
dAdt[j+jmax+1] = (1/(value_type((2*(2*j+3)))))*RHStemp1;
}
}
int main()
{
const state_type initial = loadInitialData(); //Initial condition
omp_set_num_threads(jmax+1);
int chunk_size = jmax/omp_get_max_threads();
omp_set_schedule( omp_sched_dynamic, chunk_size );
//I define my controlled error steppers
typedef runge_kutta_fehlberg78< state_type , value_type ,
state_type , value_type,openmp_range_algebra> error_stepper_type;
typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
controlled_stepper_type controlled_stepper;
int steps = integrate_adaptive( controlled_stepper ,rhs ,
initial, TINITIAL , TFINAL,INITIAL_STEP , push_back_state_and_time( A_vec , times ) );
} |
Je ne montre pas ici la definition de toutes les variables, mais je doute qu'elles soient la source du probleme puisque tout marche tres bien si je supprime l'option openmp_range_algebra de la definition de error_stepper_type. Le programme tourne aussi si j'utilise openmp_range_algebra avec un stepper de taille constante, comme le scheme Runge Kutta d'ordre 4.
Cependant avec mon code, j'obtiens l'erreur suivante :
Code:
1
2
|
invalid conversion from 'boost::range_iterator<const boost::array<double, 22ull>, void>::type {aka const double*}' to 'boost::range_iterator<boost::array<double, 22ull>, void>::type {aka double*}' [-fpermissive]| |
Il me semble aue le compilateur essaie de modifier une quantite qui serait constante. Cette erreur apparait dans le fichier openmp_range_algebra.hpp, dans le code suivant :
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
|
template< class S >
static typename norm_result_type< S >::type norm_inf( const S &s )
{
using std::max;
using std::abs;
typedef typename norm_result_type< S >::type result_type;
result_type init = static_cast< result_type >( 0 );
const size_t len = boost::size(s);
typename boost::range_iterator<S>::type beg = boost::begin(s);
#pragma omp parallel for reduction(max: init) schedule(dynamic)
for( size_t i = 0 ; i < len ; ++i )
init = max( init , abs( beg[i] ) );
return init;
} |
J'espere avoir ete clair, j'aimerais reussir a utiliser un stepper adaptatif avec la parallelisation du code.
Merci beaucoup pour votre attention et votre aide.
