Problème programme Fortran90

Version imprimable

Bonjour à tous,
J'ai pour mes études un programme en fortran90 à faire,
Je dois simuler des particules avec une trajectoire rectiligne aléatoire et calculer le temps de la 1ère collision entre 2 particules ainsi que les nouvelles trajectoires de ces deux particules.
Le programme doit être en 2 dimensions, cependant j'ai un problème lors du calcul du temps de collision

Voici mon 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
program tp6
  implicit none
  integer n, t, tmax, i, j, icol, jcol, k, p
  parameter(p=2,n=p**3,tmax=1000)
  real(8) tfin,L,r(1:2,0:3,1:n),rij(1:2,1:n,1:n),h2,rij2(1:n,1:n),&
       x(1:2,1:n),y(1:2,1:n),pi,h1,v(1:2,1:n),v0,tcol,delta, tij
 
  parameter(tfin=1.0,h2=(tfin/tmax)**2,h1=sqrt(h2),&
       L=2**(1.d0/6)*p,pi=2*asin(1.d0),v0=0.5,delta=3.0)
 
  forall(k=1:2,i=1:n) r(k,0,i)=modulo((i-1)/p**(k-1),p)
  r=r*L/p; r=r+0.5
  call random_number(x); call random_number(y)
  v=v0*cos(2*pi*x)*sqrt(-2*log(1-y))
  forall(k=1:2) v(k,:)=v(k,:)-sum(v(k,:))/n
  r(:,1,:)=r(:,0,:)+h1*v(:,:)
  ! Nouvelles trajectoires : (v1' = v2 , v2' = v1 )
 
 
 
 
  ! Temps de collision :
  do t=1, tmax
     do i = 1,n;
        do j=i+1,n
           tij = min((delta - r(:,i) + r(:,j) ) / (v(1,i)-v(1,j)) , (-delta -r(:,i)+r(:,j))/(v(1,i) - v(1,j)))  ! On prend la valeur de t la plus petite. Solution d'équation du second degré.   !! Il y a des erreurs de code sur celle ligne" Rank mismatch in array reference at (1) (2/3)" 
 
           if (( tij>0) .and. (tij<tcol)) then  ! Temps de collision des deux premières particules 
              tcol = tij
              icol=i; jcol=j
           endif
        enddo
     enddo
  enddo
end program tp6

J'ai le message d'erreur suivant :
tij = min((delta - r(:,i) + r(:,j) ) / (v(1,i)-v(1,j)) , (-delta -r(:,i)+r(:,j))/(v(1,i) - v(1,j)))
1
Error: Rank mismatch in array reference at (1) (2/3)

Je n'arrive pas à résoudre le problème, j'ai essayé avec plusieurs choses dans les parenthèses de v( ) et r( ) mais toujours un problème..
Étant étudiant en physique et n'ayant pas eu de cours de fortran90, serait il possible que quelqu'un m'explique le problème pour que je puisse le corriger ?
J'ai réussi à résoudre ce problème en 1 dimension mais j'ai un problème lorsqu'il faut passer à la dimension supérieure
Code en 1 D :

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
program projet2
integer n, t, tmax, i, j, icol, jcol
parameter (tmax=10, n=6)
real r(1:n), rij(1:n,1:n), delta, v0, vi, tij, v(1:n), x(1:n), y(1:n), tcol, pi ! il faut que le delta soit > 0 pour que la collision ai lieu a t > 0
! x et y sont les nombres aléatoires qu'on va générer
 
parameter (delta=3.0, v0=0.5, pi=2*asin(1.0)) ! delta > 0 
r=2*delta*(/ (i, i=1,n) /)
r=r-sum(r)/n
call random_number(x) ; call random_number (y)
v = v0*cos(2*pi*x)*sqrt(-2*log(1-y))
v=v-sum(v)/n 
do t=1, tmax
   tcol = r(n)/v0
   do i=1,n;
      do j=i+1,n
      tij = min((delta-r(i) + r(j) )/(v(i)-v(j)) , (-delta-r(i)+r(j))/(v(i)-v(j))) ! On prend la valeur de t la plus petite. Solution d'équation du second degré. 
      if (( tij>0) .and. (tij<tcol)) then   ! On définit le temps de colision des deux premières particules. 
         tcol = tij
         icol=i; jcol=j
      endif
   enddo
enddo
enddo
end program projet2

Merci d'avance pour votre réponse

10/04/2017, 07h49
Ehouarn

Bonjour,

Citation:

Envoyé par Kerbor

J'ai le message d'erreur suivant :
tij = min((delta - r(:,i) + r(:,j) ) / (v(1,i)-v(1,j)) , (-delta -r(:,i)+r(:,j))/(v(1,i) - v(1,j)))
1
Error: Rank mismatch in array reference at (1) (2/3)

Ben oui, tu déclare un tableau r(1:2,0:3,1:n) à trois dimensions, alors que dans la ligne de code ci-dessus tu fais comme si il n'en avait que deux.

Merci pour la réponse !
Le programme compile ! :D

J'ai ceci au final

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
program projet
  implicit none
  integer n, t, tmax, i, j, icol, jcol, k, p    ! n = nombre de particules 
  parameter(p=2,n=p**3,tmax=1000)
  real(8) tfin,L,r(1:2,1:n),h2, x(1:2,1:n),y(1:2,1:n),pi,h1,v(1:2,1:n),&  
       v0,tcol,delta, tij
 
  parameter(tfin=1.0,h2=(tfin/tmax)**2,h1=sqrt(h2),&
       L=2**(1.d0/6)*p,pi=2*asin(1.d0),v0=0.5,delta=3.0)
 
  forall(k=1:2,i=1:n) r(k,i)=modulo((i-1)/p**(k-1),p)
  r=r*L/p; r=r+0.5
  call random_number(x); call random_number(y)
  v=v0*cos(2*pi*x)*sqrt(-2*log(1-y))
  forall(k=1:2) v(k,:)=v(k,:)-sum(v(k,:))/n
  r(:,:)=r(:,:)+h1*v(:,:)
  ! Temps de collision :
  do t= 1,tmax 
     do i = 1,n;
        do j=i+1,n
           tij = min((delta - r(2,i) + r(2,j) ) / (v(2,i)-v(2,j)) , (-delta -r(2,i)+r(2,j))/(v(2,i) - v(2,j)))
           ! On prend la valeur de t la plus petite. Solution d'équation du second degré.          
           if (( tij>0) .and. (tij<tcol)) then  ! Temps de collision des deux premières particules 
              tcol = tij
              icol=i; jcol=j
           endif 
        enddo
     enddo 
  enddo
  ! Nouvelles trajectoires : (v1' = v2 , v2' = v1 )
  if (tcol == tij ) then 
     v(1:2,i) = v(1:2,j)
     v(1:2,j) = v(1:2,i)
  endif
 
end program projet

15/04/2017, 10h03
Ehouarn
Bonjour,

Très bien si ça compile, mais est-ce que ça donne des bons résultats?
Sans avoir décortiqué ton code, j'ai comme un doute quand je vois ça:
Citation:
Envoyé par Kerbor

Code:

1 2 3 4 5 6 ! Nouvelles trajectoires : (v1' = v2 , v2' = v1 ) if (tcol == tij ) then v(1:2,i) = v(1:2,j) v(1:2,j) = v(1:2,i) endif
On dirait que tu tentes une permutation entre v(1:2,i) et v(1:2,j) alors qu'en fait tu ne fait que mettre v(1:2,i) à v(1:2,j), donc de fait égaliser les deux.
15/04/2017, 12h29
Kerbor
Bonjour !

Oui j'y ai pensé aussi.. Je ne sais pas comment faire pour écrire que la vitesse de la particule 1 après collision est égale à la vitesse de la particule 2 avant collision et inversement
Code:

1 2 --> v1' = v2 v2' = v1
Il faudrait que je définisse une nouvelle variable du coup ?

Citation:

real vaprès(1:2,1:n)

vaprès(1:2,i) = v (1:2,j)
vaprès(1:2,j) = v (1:2,i)

Avec ceci je pense que le problème serait résolu ?

Oui il faut passer par une variable intermédiaire, du style:

Code:

1
2
3
4
5
6
7
8
 
real v_tmp(1:2)
 
...
! permutation des contenus de v(1:2,i) et v(1:2,j)
v_tmp(1:2) = v(1:2,i)
v(1:2,i) = v(1:2,j)
v(1:2,j) = v_tmp(1:2)

15/04/2017, 17h27
Kerbor
Oh d'accord mais du coup le programme ne fait pas
v(1:2,j) = v_tmp(1:2) = v(1:2,i) ? Ce qui reviendrait au même problème :/

Si ça garde le v(1:2,i) initial dans le v_col(1:2) ?
Code:

1 2 3 4 5 6 ! Nouvelles trajectoires : (v1' = v2 , v2' = v1 ) if (tcol == tij ) then v_col(1:2)=v(1:2,i) v(1:2,i) = v(1:2,j) v(1:2,j) = v_col(1:2) endif
16/04/2017, 09h20
Ehouarn

Je ne comprends pas ton interrogation.
Ton code fait la même permutation circulaire que ce que j'ai proposé.

Je me suis embrouillé pardon :)

Je pense avoir fini avec ce programme ! Merci beaucoup pour ton aide !

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
program projet
  implicit none
  integer n, t, tmax, i, j, icol, jcol, k, p    ! n = nombre de particules 
  parameter(p=2,n=p**3,tmax=1000)
  real(8) tfin,L,r(1:2,1:n),h2, x(1:2,1:n),y(1:2,1:n),pi,h1,v(1:2,1:n),&
       v0,tcol,delta, tij, v_col(1:2)
 
  parameter(tfin=1.0,h2=(tfin/tmax)**2,h1=sqrt(h2),&
       L=2**(1.d0/6)*p,pi=2*asin(1.d0),v0=0.5,delta=3.0)
 
  forall(k=1:2,i=1:n) r(k,i)=modulo((i-1)/p**(k-1),p)
  r=r*L/p; r=r+0.5
  call random_number(x); call random_number(y)
  v=v0*cos(2*pi*x)*sqrt(-2*log(1-y))
  forall(k=1:2) v(k,:)=v(k,:)-sum(v(k,:))/n
  r(:,:)=r(:,:)+h1*v(:,:)
 
  ! Temps de collision :
  do t= 1,tmax 
     do i = 1,n-1
        do j=i+1,n
           tij = min((delta - r(2,i) + r(2,j) ) / (v(2,i)-v(2,j)) , (-delta -r(2,i)+r(2,j))/(v(2,i) - v(2,j))) ! On prend la valeur de t la plus petite. Solution d'équation du second degré.          
           if (( tij>0) .and. (tij<tcol)) then  ! Temps de collision des deux premières particules 
              tcol = tij
              icol=i; jcol=j  ! Coordonnées de la collision
           endif          
       enddo
     enddo
  enddo
 
  ! Nouvelles trajectoires : (v1' = v2 , v2' = v1 )
  if (tcol == tij ) then
     v_col(2)=v(2,i) 
     v(2,i) = v(2,j)
     v(2,j) = v_col(2)
  endif
 
  do t= 1,tmax 
     do i = 1,n-1
        do j=i+1,n
           write(6+i,*) v(2,i)
        enddo
     enddo
  enddo
 
end program projet
 
! plot 'fort.7' w l,'fort.8' w l,'fort.9' w l,'fort.10' w l,'fort.11' w l,'fort.12' w l,'fort.13' w l,'fort.14' w l, 'fort.15' w l,'fort.16' w l,'fort.17' w l,'fort.18' w l,'fort.19' w l,'fort.20' w l,'fort.21' w l