Discussion :

[Nearos]

Bonjour,

Je suis nouveau au sein du forum, donc excusez-moi d'avance si le problème a pu être posé auparavant. Cela fait pas mal de temps que je suis coincé dans le plot d'un programme. Il s'agit d'un système couplé où le but est d'afficher les valeurs d'une concentration et d'une densité en utilisant des différences finies.
Le schéma est représenté en bas ( en screen ).

Voici mon code pour le cas " theta = 1 ":

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
program theta0
  implicit none
 
  real*8, dimension(:,:), allocatable :: R,T1,delta,phi,C,A,B,delta1
  real*8 :: c0,c1,c2,alpha,beta,r0,dt,dx,Tmax,l,g
  integer :: nb_iter,m,i,j
 
 
 
  print*, ' m = '
  read*, m
!100
  print*, ' nb_iter = '
  read*, nb_iter
!qu'importe
 
 
  allocate(R(0:nb_iter,0:m))
  allocate(delta(0:nb_iter-1,1:m-1))
  allocate(T1(1:nb_iter,1:m-1)) 
  allocate(phi(0:nb_iter,0:m-1))
  allocate(C(0:nb_iter,0:m-1))
  allocate(B(0:nb_iter,0:m-1))
  allocate(delta1(0:nb_iter-1,1:m-1))
 
 
 
! Constantes
 
  c0 = 10.
 
  alpha = 0.01
 
  beta = 0.1
 
  r0 = 1.
 
  Tmax = 0.5
 
 
 
 
 
! Conditions initiales
 
  C(0,1:m-1) = c0
  phi(0,1:m-1) = alpha*c0 + beta
  R(0,1:m-1) = 0.
  B(0,1:m-1) = 0.
 
! Conditions aux limites en 0
 
  R(0:nb_iter,0) = r0
  do i = 0,nb_iter
    C(i,0) = c0*exp(-r0*i/10)
    phi(i,0) = alpha*C(i,0) + beta
    B(i,0) = r0/phi(i,0)
  end do
 
 
! Conditions aux limites en m
 
  R(0:nb_iter,m) = R(0:nb_iter,m-1)
 
 
 
dt = Tmax/nb_iter
 
dx = 1./m
 
!B(0:nb_iter,0:m) = R(0:nb_iter,0:m)/phi(0:nb_iter,0:m)
 
 
 
do i = 1,nb_iter
  do j = 2,m
 
  C(i,j-1) = C(i-1,j-1)*exp(-dt*R(i-1,j-1))
 
  phi(i,j-1) = alpha*C(i,j-1) + beta
 
  delta(i-1,j-1) = (phi(i-1,j-1)+phi(i-1,j))*(B(i-1,j)-B(i-1,j-1))-(phi(i-1,j-2)+phi(i-1,j-1))*(B(i-1,j-1)-B(i-1,j-2))
 
  delta1(i-1,j-1) = (dt/dx**2)*delta(i-1,j-1)
 
  T1(i,j-1) = 1./(1-dt*C(i,j-1))
 
  R(i,j-1) = T1(i,j-1)*(R(i-1,j-1)+delta1(i-1,j-1))
 
  B(i,j-1) = R(i,j-1)/phi(i,j-1)
 
  end do
end do
 
 
 
open(unit=10,file='theta1.dat')
 
do i = 0,m-1
  write(10,*) i*dx, R(0,i), R(1,i), R(5,i), R(10,i), R(20,i), R(50,i), R(300,i), C(0,i), C(1,i), C(3,i), C(5,i), C(10,i), C(20,i)
 
end do
 
 
 
end program theta0
 
 
 
end program theta0

La densité étant la représentative de "C", et la concentration celle de "Rou", " phi " vaut " alpha*C + beta ", donnée non explicitée dans les schémas, la notation C_n_m signifie qu'on se trouve à la position m au moment n.

En les affichant en fixant le temps tout en faisant varier l'espace, je devrais obtenir les deux courbes affichées en bas de page ( celle de gauche représente la concentration "R", et celle de droite la densité "C" ), ( veuillez ne pas tenir compte des courbes plottées en utilisant les élèments finies ).

J'obtiens pour ma part des courbes totalement hors contextes, malgré que mon code semble cohérent ( des discontinuités incompréhensibles par moment ..).

C'est dans ce contexte que je suis venu chercher de l'aide, ou potentiellement, choper des tuiles pour pouvoir avancer.

Merci pour le temps que vous aurez consacré.

[Ehouarn]

Bonjour,

Si tu compiles ton code avec des options de débogage, alors tu trouveras vite le problème (et c'est vraiment LA bonne habitude à prendre):
Par exemple:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
gfortran -ffpe-trap=invalid,zero,overflow -fbounds-check -g3 -O0 -fstack-protector-all -fbacktrace theta0.f90 -o theta0

Et à l'exécution (en prenant au pif m=50 et nb_iter=100), ça plante avec:
At line 82 of file theta0.f90
Fortran runtime error: Array reference out of bounds for array 'phi', upper bound of dimension 2 exceeded (50 > 49)