Interpolation linéaire et grille

**abs11** · 10/01/2012, 01h10

Bonsoir
je veux interpolé les valeurs d'un champs, donné sur une grille régulière, dans une autre grille raffinée.
Pour ceci, j'ai essayé d'employer la subroutine de l'interpolation linéaire donnée dans ce site, mais je n'arrive pas à trouver des résultats.
Voici le l'exemple que j'ai utilisé

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
Program main
	implicit real*8(A-H,O-Z)
	Dimension xp(50),yp(50)
	dimension x1(100),y1(100)
	dimension u(100,100)
	dimension a(100),b(100)
	integer,dimension(2,2):: POSXY
	real(kind=8),dimension(size(xp)):: TEMPMAILLX
	real(kind=8),dimension(size(yp)):: TEMPMAILLY
	real(kind=8)::INTERPOL
	real(kind=8) :: HI,HJ,INTERPOL1,INTERPOL2
 
	write(6,*) 'read the value of n,m,nn,mm'
	read(5,*) n,m,nn,mm
 
	open(2,file='uregul.dat')
	open(3,file='INTERPOL.dat')		
	u(:,:)=0.0
	call grid(n,m,nn,mm,xp,yp,x1,y1)

Le champs donné

Code :

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

1
2
3
4
5
6
7
 
        do j=1,m
	do i=1,n
	u(i,j)=0.1*(1.-((2.*yp(j)/20.)-1.)**2.)
	write(2,*)xp(i),yp(j),u(i,j)
	end do
	end do

l'interpolation

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
 
do j=1,mm
	do i=1,nn
 
	TEMPMAILLX=xp
	!POSXY(1,1)=minloc(abs(x1(i)-TEMPMAILLX),dim=1)
	!if(minval(abs(x1(i)-TEMPMAILLX))/=0.0) TEMPMAILLX(POSXY(1,1))=
  !   &      TEMPMAILLX(POSXY(1,1))-1.
 
		POSXY(1,1)=minloc(abs(x1(i)-xp),dim=1)
	if(minval(abs(x1(i)-xp))/=0.0) xp(POSXY(1,1))=
     &      xp(POSXY(1,1))-1.
 
	POSXY(2,1)=minloc(abs(x1(i)-TEMPMAILLX),dim=1)
 
	TEMPMAILLY=yp
	POSXY(1,2)=minloc(abs(y1(j)-TEMPMAILLY),dim=1)
	if(minval(abs(y1(j)-TEMPMAILLY))/=0.0) TEMPMAILLY(POSXY(1,2))=
     &	  TEMPMAILLY(POSXY(1,2))+0.5
	POSXY(2,2)=minloc(abs(y1(j)-TEMPMAILLY),dim=1)
 
	if(minval(abs(x1(i)-TEMPMAILLX))/=0.0) HI=-(xp(POSXY(1,1))-
     &      x1(i))/(xp(POSXY(2,1))-xp(POSXY(1,1)))
	if(minval(abs(y1(j)-TEMPMAILLY))/=0.0) HJ=-(yp(POSXY(1,2))-
     &      y1(j))/(yp(POSXY(2,2))-yp(POSXY(1,2)))
 
 
	INTERPOL1=(1-HI)*u(POSXY(1,1),POSXY(1,2))+HI*u(POSXY(2,1)
     &      ,POSXY(1,2))
 
	INTERPOL2=(1-HI)*u(POSXY(1,1),POSXY(2,2))+HI*u(POSXY(2,1)
     &      ,POSXY(2,2))
 
	INTERPOL=(1-HJ)*INTERPOL1+HJ*INTERPOL2
 
 
      !do i=1,nn
	!do j=1,mm
	!write(3,*)x1(i),y1(j),INTERPOL
	end do
	end do
 
stop
end

la grille régulière et grille raffinée

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
 
 
 
 
	subroutine grid(n,m,nn,mm,xp,yp,x1,y1)
	implicit real*8(A-H,O-Z)
	Dimension xp(50), yp(50)
	dimension x1(100),y1(100)
 
*****generation des x ******
	do 21 i=1,n
	alpha= float(i-1)
	xp(i)=alpha
 21	 continue
	write(6,102)
102	format(15x,'		Horizontal coordinate		')
	write(6,100) (xp(j),j=1,n)
100	format(1x,6f12.6)
 
****	génération des cordonées y*****
 
	do 31 j=1,m
 
	alpha= float(j-1)
	yp(j)=alpha
31	continue
	write(6,103)
103	format(15x,'		longitudinal coordinate		')
	write(6,100)(yp(j),j=1,m)
 
 
***** maillage raffiné( insertion de nouveaux noeuds)******
 
	tdx=(xp(n)-xp(1))/float(nn-1)
	tdy=(yp(m)-yp(1))/float(mm-1)
	do 41 i=1,nn
41	x1(i)=float(i-1)*tdx
	do 51 j=1,mm
51	y1(j)=float(j-1)*tdy
 
 
	return 
	end

Merci de votre aide je suis nouveau sur fortran et je n'arrive pas à résoudre ce problème
Merci encore

**Marlan** · 11/01/2012, 09h06

Salut,

Pourquoi ne gardes tu pas la routine d'interpolation comme une routine?

Un exemple d'utilisation :

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
 
PROGRAM Raff
IMPLICIT NONE
REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: TAB, TAB_RAFF !Un tableau et !un tableau raffiné correspondant aux valeurs de champs 2D sur le maillage défini !ci-dessous
REAL(KIND=8), ALLOCATABLE, DIMENSION(:) :: MX, MY, MX_RAFF, MY_RAFF !Points de grille et points de grille raffinés
INTEGER :: NI=10, NJ=11, NI_RAFF=20, NJ_RAFF=22 !Nombre de points de !grille associés au maillage et à la version raffinée
INTEGER :: I,J,K !Les très classiques indices de boucles
 
INTERFACE
subroutine INTERPOL_BARI(X,CHAMP,MAILLX,MAILLY,INTERPOL)
!==========INTERPOL_BARI===============
!SUBROUTINE permettant l'interpolation lineaire en un point donne
!afin de calculer la valeur d'un vecteur en ce meme point.
!Entree : X, vecteur position, MAILL : maillage, CHAMP : les valeur du champ
!dont on souhaite connaitre la valeur en X.
!Sortie : INTERPOL: la valeur interpolee du champ a la position X
!Marlan, 2011
implicit none
	real(kind=8),dimension(:),intent(in) :: X, MAILLX, MAILLY
	real(kind=8),dimension(:,:),intent(in) :: CHAMP
	real(kind=8),intent(out)::INTERPOL
END SUBROUTINE INTERPOL_BARI
END INTERFACE
 
!On alloue les tableaux dynamiques correspondant au maillage
ALLOCATE(MX(NI),MY(NJ),MX_RAFF(NI_RAFF),MY_RAFF(NJ_RAFF))
 
!On alloue les tableaux dynamiques correspondant aux valeur du champ
ALLOCATE(TAB(NJ,NI),TAB_RAFF(NJ_RAFF,NI_RAFF))
TAB=0.0
 
!On remplit TAB de valeurs bidons pour l'exemple
DO I=1,size(TAB,dim=1)
MX(I)=I
end do
 
DO J=1,size(TAB,dim=2)
MY(J)=J
END DO
 
DO I=1,size(TAB,dim=1)
DO J=1,size(TAB,dim=2)
TAB(I,J)=I+J
END DO
end do
 
DO I=2,size(TAB_RAFF,dim=1)
MX_RAFF(I)=I/2.0
end do
DO J=2,size(TAB_RAFF,dim=2)
MY_RAFF(J)=J/2.0
END DO
 
print'(11f5.1)',TAB
!On veut maintenant avoir les vlaeurs de TAB_RAFF en fonction de celles de TAB en !utilisant la routine INTERPOL_BARI
DO I=1,size(TAB_RAFF,dim=1)
DO J=1,size(TAB_RAFF,dim=2)
CALL INTERPOL_BARI((/MX_RAFF(I),MY_RAFF(J)/),TAB,MX,MY,TAB_RAFF(I,J))
END DO
END DO
print*,''
print'(22f5.1)',TAB_RAFF
END PROGRAM Raff
 
 
 
 
 
 
 
 
subroutine INTERPOL_BARI(X,CHAMP,MAILLX,MAILLY,INTERPOL)
!==========INTERPOL_BARI===============
!SUBROUTINE permettant l'interpolation lineaire en un point donne
!afin de calculer la valeur d'un vecteur en ce meme point.
!Entree : X, vecteur position, MAILL : maillage, CHAMP : les valeur du champ
!dont on souhaite connaitre la valeur en X.
!Sortie : INTERPOL: la valeur interpolee du champ a la position X
!Marlan, 2011
implicit none
	real(kind=8),dimension(:),intent(in) :: X, MAILLX, MAILLY
	real(kind=8),dimension(:,:),intent(in) :: CHAMP
	real(kind=8),intent(out)::INTERPOL
	real(kind=8),dimension(size(MAILLX)):: TEMPMAILLX
	real(kind=8),dimension(size(MAILLY)):: TEMPMAILLY
	real(kind=8) :: HI,HJ,INTERPOL1,INTERPOL2
	integer,dimension(2,2):: POSXY
 
	TEMPMAILLX=MAILLX
	POSXY(1,1)=minloc(abs(X(1)-TEMPMAILLX),dim=1)
	if(minval(abs(X(1)-TEMPMAILLX))/=0.0) TEMPMAILLX(POSXY(1,1))=TEMPMAILLX(POSXY(1,1))+10.0
	POSXY(2,1)=minloc(abs(X(1)-TEMPMAILLX),dim=1)
	TEMPMAILLY=MAILLY
	POSXY(1,2)=minloc(abs(X(2)-TEMPMAILLY),dim=1)
	if(minval(abs(X(2)-TEMPMAILLY))/=0.0) TEMPMAILLY(POSXY(1,2))=TEMPMAILLY(POSXY(1,2))+10.0
	POSXY(2,2)=minloc(abs(X(2)-TEMPMAILLY),dim=1)
 
	if(minval(abs(X(1)-TEMPMAILLX))/=0.0) HI=-(MAILLX(POSXY(1,1))-X(1))/(MAILLX(POSXY(2,1))-MAILLX(POSXY(1,1)))
	if(minval(abs(X(2)-TEMPMAILLY))/=0.0) HJ=-(MAILLY(POSXY(1,2))-X(2))/(MAILLY(POSXY(2,2))-MAILLY(POSXY(1,2)))
 
 
	INTERPOL1=(1-HI)*CHAMP(POSXY(1,1),POSXY(1,2))+HI*CHAMP(POSXY(2,1),POSXY(1,2))
 
	INTERPOL2=(1-HI)*CHAMP(POSXY(1,1),POSXY(2,2))+HI*CHAMP(POSXY(2,1),POSXY(2,2))
 
	INTERPOL=(1-HJ)*INTERPOL1+HJ*INTERPOL2
end subroutine INTERPOL_BARI

ATTENTION : Les valeurs de TAB doivent être rangées par : colonnes pour la direction 'X', ligne pour la direction 'Y' pour les interpolations.

En espérant que ca t'aurras débloqué,

**abs11** · 11/01/2012, 14h25

Merci beaucoup vous m'avez sauvez la vie !!!!!!!!!
Juste une autre petite question svp
j'ai voulu représenter mon champs (TAB) en fonction des coordonnées de la grille
(Mx,My) donc j'ai enregistré les valeurs dans un fichier classique
do I=1,NI
do J=1,NJ
write(2,*) Mx(I),My(J),TAB(I,J)
end do
end do
mais la dernière valeurs de de NJ n'est pas balayé!!
merci encore de votre aide

**Marlan** · 11/01/2012, 14h32

Essayes avec TAB(J,I) ça devrait mieux marcher (principe de ranger les éléments suivant Y dans les lignes de TAB, les éléments suivant X dans les indices variable de colonne). Je précise que ce n'est là qu'une façon de coder ce genre de problème, pas un absolu. (je le préfère pour ma pars car il correspond aux représentation graphiques que l'on peut se faire de ce genre de problème)

**abs11** · 11/01/2012, 14h51

j'ai obtenu le bon profil
en faisant
do I=1,NI
do J=1,NJ
write(2,*)MY(I),MX(J),TAB(I,J)
end do
end do

**Marlan** · 11/01/2012, 15h00

Heu... tu es sûr que ce ne serait pas plutôt :

Code :

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

1
2
3
4
5
do I=1,NJ
do J=1,NI
write(2,*)MY(I),MX(J),TAB(I,J)
end do
end do

Car MY a NJ éléments dans l'exemples que je t'ai donné plus haut et MX NI éléments...

Enfin bon, content d'avoir pu t'aider. N'oublies pas de cliquer sur le 'Résolu' si tu n'as plus de soucis.

@ Peluche,

Marlan

**abs11** · 11/01/2012, 15h18

je pense c 'est du au fait qu'en allouant le tableau TAB(NJ,NI)
puis on a rempli les MX par :
do I=1,size(TAB,dim=1)
et les MY par:
do j=1,size(TAB,dim=2)

dans tous le cas mercii de votre aide c'était très important pour moi!!!

**abs11** · 12/01/2012, 12h09

Bonjour
J'ai essayé d'insérer ce programme dans mon programme initial pour faire l'interpolation d'un champs déja calculé par le programme initial
J'ai voulu alors mettre mon ancien programme sous forme d'une subroutine
donc j'aurai:
1-Déclaration des variables
2-Construction de la grille MX,MY
3-arrivant à l'instruction ou je dois remplir mon champs TAB j'appelle ma subroutine (ancien programme) qui calcule le champs de temperature
T(I,J)=TAB(I,J)

Mais dans ce cas ma subroutine n'a pas été lu par le programme
et la notion d'interface que j'ai pas pu encore maîtriser ne ma pas laisser déclarer aucune nouvelle variable!!
Merci de votre aide

**Marlan** · 13/01/2012, 11h11

Est ce que tu pourrais nous montrer ton code?

A moins que le soucis ne soit résolu?

@ plus,

Marlan

**abs11** · 13/01/2012, 11h29

Bonjour
Merci beaucoup de votre aide,j'ai pu résoudre ce souci mais j'ai rencontré un autre problème

En prenant NI différent de NJ:
(NI=50,NJ=100) un message d'erreur s'affiche!! " no matching symbolic information found"
Voici le code avec l'interpolation qui marche très bien pour NI=NJ
mais pas pour NI différent de NJ.
Et puisque pour mon cas, j'utilise un canal rectangulaire, donc je dois avoir NI beaucoup plus supérieure que NJ.

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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
	PROGRAM Raff
	IMPLICIT NONE
	parameter n=100,m=100,NI=100,NJ=100, NI_RAFF=200, NJ_RAFF=200
	REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: TAB, TAB_RAFF !Un tableau et 
	!un tableau raffiné correspondant aux valeurs de champs 2D sur le maillage défini !ci-dessous
	REAL(KIND=8), ALLOCATABLE, DIMENSION(:) :: MX, MY, MX_RAFF, 
     &      MY_RAFF 
	real f1(0:n,0:m),f2(0:n,0:m),f3(0:n,0:m),f4(0:n,0:m)
	real rho(0:n,0:m),feq
	!Points de grille et points de grille raffinés
!	INTEGER :: NI=10, NJ=10, NI_RAFF=20, NJ_RAFF=20 !Nombre de points de 
	!grille associés au maillage et à la version raffinée
	INTEGER :: I,J,K !Les très classiques indices de boucles
 
 
	INTERFACE
	subroutine INTERPOL_BARI(X,CHAMP,MAILLX,MAILLY,INTERPOL)
	!==========INTERPOL_BARI===============
	!SUBROUTINE permettant l'interpolation lineaire en un point donne
	!afin de calculer la valeur d'un vecteur en ce meme point.
	!Entree : X, vecteur position, MAILL : maillage, CHAMP : les valeur du champ
	!dont on souhaite connaitre la valeur en X.
	!Sortie : INTERPOL: la valeur interpolee du champ a la position X
	!Marlan, 2011
	implicit none
	real(kind=8),dimension(:),intent(in) :: X, MAILLX, MAILLY
	real(kind=8),dimension(:,:),intent(in) :: CHAMP
	real(kind=8),intent(out)::INTERPOL
	END SUBROUTINE INTERPOL_BARI
	END INTERFACE
	open (2,file='champ.dat')
	open (3,file='champinter.dat')
         open(4,file='ro.dat')
	!On alloue les tableaux dynamiques correspondant au maillage
	ALLOCATE(MX(NI),MY(NJ),MX_RAFF(NI_RAFF),MY_RAFF(NJ_RAFF))
 
	!On alloue les tableaux dynamiques correspondant aux valeur du champ
	ALLOCATE(TAB(NI,NJ),TAB_RAFF(NI_RAFF,NJ_RAFF))
	TAB=0.0
 
	!On remplit TAB de valeurs bidons pour l'exemple
	DO I=1,size(TAB,dim=1)
	MX(I)=I
	print*,I,MX(I)
	end do
 
	pause
	DO J=1,size(TAB,dim=2)
	MY(J)=J
	print*,J,MY(J)
 
	END DO
	pause
 
	call dab(n,m,f1,f2,f3,f4,rho,feq)
      DO  I=1,NI
	DO J=1,NJ 
 
	TAB(I,J)=rho(I,J)
 
	end do
 
	END DO
 
 
 
 
	continue
      DO I=1,NI
	DO J=1,NJ
100	format(1X,3f12.6,1X)
	write(2,100) MY(I),MX(J),TAB(I,J)
	end do
	end do
	DO I=1,size(TAB_RAFF,dim=1)
	MX_RAFF(I)=I/2.0
	end do
	DO J=1,size(TAB_RAFF,dim=2)
	MY_RAFF(J)=J/2.0
	END DO
 
 
	!On veut maintenant avoir les vlaeurs de TAB_RAFF en fonction de celles de TAB en !utilisant la routine INTERPOL_BARI
	DO I=1,size(TAB_RAFF,dim=1)
	DO J=1,size(TAB_RAFF,dim=2)
	CALL INTERPOL_BARI((/MX_RAFF(I),MY_RAFF(J)/),TAB,MX,MY,
     &      TAB_RAFF(I,J))
	END DO
	END DO
 
	 DO I=1,NI_RAFF
	DO J=1,NJ_RAFF
	write(3,100) MX_RAFF(I),MY_RAFF(J),TAB_RAFF(I,J)
	end do
	end do
	END PROGRAM Raff
 
 
 
 
 
c***********************************************************************************************
 
 
 
	subroutine dab(n,m,f1,f2,f3,f4,rho,feq)
 
	!parameter (n=100,m=200)
	real f1(0:n,0:m),f2(0:n,0:m),f3(0:n,0:m),f4(0:n,0:m)
	real rho(0:n,0:m),feq
	!integer i,j,kk
 
	dx=1.
	dy=dx
	dt=1.
 
	csq=dx*dx/(dt*dt)
	alpha=0.25
	omega=1.0/(2.*alpha/(dt*csq)+0.5)
	mstep=400
	do j=0,m
	do i=0,n
	rho(i,j)=0.0
	end do
	end do
 
	do j=0,m
	do i=0,n
	f1(i,j)=0.25*rho(i,j)
	f2(i,j)=0.25*rho(i,j)
	f3(i,j)=0.25*rho(i,j)
	f4(i,j)=0.25*rho(i,j)
	end do
	end do
 
 
	do kk=1, mstep
	do j=0,m
	do i=0,n
	feq=0.25*rho(i,j)
	f1(i,j)=omega*feq+(1.-omega)*f1(i,j)
	f2(i,j)=omega*feq+(1.-omega)*f2(i,j)
	f3(i,j)=omega*feq+(1.-omega)*f3(i,j)
	f4(i,j)=omega*feq+(1.-omega)*f4(i,j)
	end do
	end do
 
 
	do j=0,m
	do i=1,n
	f1(n-i,j)=f1(n-i-1,j)
	f2(i-1,j)=f2(i,j)
	end do
	end do
	do i=0,n
	do j=1,m
	f3(i,m-j)=f3(i,m-j-1)
	f4(i,j-1)=f4(i,j)
	end do
	end do
 
 
 
 
	do j=0,m
	f1(0,j)=0.5-f2(0,j)
	f3(0,j)=0.5-f4(0,j)
	f1(n,j)=0.0
	f2(n,j)=0.0
	f3(n,j)=0.0
	f4(n,j)=0.0
	end do
 
	do i=0,n
	f1(i,m)=0.0
	f2(i,m)=0.0
	f3(i,m)=0.0
	f4(i,m)=0.0
	f1(i,0)=f1(i,1)
	f2(i,0)=f1(i,1)
	f3(i,0)=f1(i,1)
	f4(i,0)=f1(i,1)
 
	end do
	do j=0,m
	do i=0,n
	rho(i,j)=f1(i,j)+f2(i,j)+f3(i,j)+f4(i,j)
	end do
	end do
	end do
	do j=0,m
	do i=0,n
	write(4,*)i,j,rho(i,j)
	end do
	end do
 
	return
	end
 
c*********************************************************
	subroutine INTERPOL_BARI(X,CHAMP,MAILLX,MAILLY,INTERPOL)
	!==========INTERPOL_BARI===============
	!SUBROUTINE permettant l'interpolation lineaire en un point donne
	!afin de calculer la valeur d'un vecteur en ce meme point.
	!Entree : X, vecteur position, MAILL : maillage, CHAMP : les valeur du champ
	!dont on souhaite connaitre la valeur en X.
	!Sortie : INTERPOL: la valeur interpolee du champ a la position X
	!Marlan, 2011
	implicit none
	real(kind=8),dimension(:),intent(in) :: X, MAILLX, MAILLY
	real(kind=8),dimension(:,:),intent(in) :: CHAMP
	real(kind=8),intent(out)::INTERPOL
	real(kind=8),dimension(size(MAILLX)):: TEMPMAILLX
	real(kind=8),dimension(size(MAILLY)):: TEMPMAILLY
	real(kind=8) :: HI,HJ,INTERPOL1,INTERPOL2
	integer,dimension(2,2):: POSXY
 
	TEMPMAILLX=MAILLX
	POSXY(1,1)=minloc(abs(X(1)-TEMPMAILLX),dim=1)
	if(minval(abs(X(1)-TEMPMAILLX))/=0.0) TEMPMAILLX(POSXY(1,1))=
     &	  TEMPMAILLX(POSXY(1,1))+10.0
	POSXY(2,1)=minloc(abs(X(1)-TEMPMAILLX),dim=1)
	TEMPMAILLY=MAILLY
	POSXY(1,2)=minloc(abs(X(2)-TEMPMAILLY),dim=1)
	if(minval(abs(X(2)-TEMPMAILLY))/=0.0) TEMPMAILLY(POSXY(1,2))=
     &      TEMPMAILLY(POSXY(1,2))+10.0
	POSXY(2,2)=minloc(abs(X(2)-TEMPMAILLY),dim=1)
 
	if(minval(abs(X(1)-TEMPMAILLX))/=0.0) HI=
     &      -(MAILLX(POSXY(1,1))-X(1))/
     &      (MAILLX(POSXY(2,1))-MAILLX(POSXY(1,1)))
	if(minval(abs(X(2)-TEMPMAILLY))/=0.0) HJ=
     &      -(MAILLY(POSXY(1,2))-X(2))/
     &      (MAILLY(POSXY(2,2))-MAILLY(POSXY(1,2)))
 
 
	INTERPOL1=(1-HI)*CHAMP(POSXY(1,1),POSXY(1,2))+HI*
     &      CHAMP(POSXY(2,1),POSXY(1,2))
 
	INTERPOL2=(1-HI)*CHAMP(POSXY(1,1),POSXY(2,2))+HI*
     &      CHAMP(POSXY(2,1),POSXY(2,2))
 
	INTERPOL=(1-HJ)*INTERPOL1+HJ*INTERPOL2
	end subroutine INTERPOL_BARI

**Marlan** · 13/01/2012, 13h59

Je n'ai pas d'erreurs pour ma pars. Cela étant, tu as inversé les indices donc les interpolations étaient fausses.

J'ai mis en rouge les modifications. Après, je n'ai pas eu l'erreur que tu as décrit... peut être que cela provient (mais j'en doute, et je ne saurai pas le dire) d'une différence de compilateur ou de langage (j'utilise ifort et fortran90 pour ma pars). Par contre, fais bien attention à déclarer les variables réelles de la même façon. Dans le cas des codes écris précédemment, elles sont toutes en double précision (kind=8). Le type réel de base est en simple précision (kind=4), valeur par défaut lorsque tu déclares un réel. Autre solution pour ne pas s'embêter à tout changer : dans les options de compilation préciser une option -r8 qui forcera tous les réels à être de double précision. (idem avec -r4)

Cela devrait mieux marcher maintenant :

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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
PROGRAM Raff
	IMPLICIT NONE
	parameter n=5,m=10,NI=5,NJ=10, NI_RAFF=10, NJ_RAFF=20
	REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: TAB, TAB_RAFF !Un tableau et 
	!un tableau raffiné correspondant aux valeurs de champs 2D sur le maillage défini !ci-dessous
	REAL(KIND=8), ALLOCATABLE, DIMENSION(:) :: MX, MY, MX_RAFF, MY_RAFF 
	real(kind=8) f1(0:n,0:m),f2(0:n,0:m),f3(0:n,0:m),f4(0:n,0:m)
	real(kind=8) rho(0:n,0:m),feq
	!Points de grille et points de grille raffinés
!	INTEGER :: NI=10, NJ=10, NI_RAFF=20, NJ_RAFF=20 !Nombre de points de 
	!grille associés au maillage et à la version raffinée
	INTEGER :: I,J,K !Les très classiques indices de boucles
	

	INTERFACE
	subroutine INTERPOL_BARI(X,CHAMP,MAILLX,MAILLY,INTERPOL)
	!==========INTERPOL_BARI===============
	!SUBROUTINE permettant l'interpolation lineaire en un point donne
	!afin de calculer la valeur d'un vecteur en ce meme point.
	!Entree : X, vecteur position, MAILL : maillage, CHAMP : les valeur du champ
	!dont on souhaite connaitre la valeur en X.
	!Sortie : INTERPOL: la valeur interpolee du champ a la position X
	!Marlan, 2011
	implicit none
	real(kind=8),dimension(:),intent(in) :: X, MAILLX, MAILLY
	real(kind=8),dimension(:,:),intent(in) :: CHAMP
	real(kind=8),intent(out)::INTERPOL
	END SUBROUTINE INTERPOL_BARI
	END INTERFACE
	open (2,file='champ.dat')
	open (3,file='champinter.dat')
         open(4,file='ro.dat')
	!On alloue les tableaux dynamiques correspondant au maillage
	ALLOCATE(MX(NI),MY(NJ),MX_RAFF(NI_RAFF),MY_RAFF(NJ_RAFF))

	!On alloue les tableaux dynamiques correspondant aux valeur du champ
	ALLOCATE(TAB(NJ,NI),TAB_RAFF(NJ_RAFF,NI_RAFF))
	TAB=0.0

	!On remplit TAB de valeurs bidons pour l'exemple
	DO I=1,size(TAB,dim=1)
	MX(I)=I
	!~ print*,I,MX(I)
	end do
	
	!~ pause
	DO J=1,size(TAB,dim=2)
	MY(J)=J
	!~ print*,J,MY(J)
	
	END DO
	!~ pause

	call dab(n,m,f1,f2,f3,f4,rho,feq)
      DO  I=1,NI
	DO J=1,NJ 
	
	TAB(J,I)=rho(I,J)
	
	end do

	END DO



	
	!~ continue
      DO I=1,NI
	DO J=1,NJ
100	format(1X,3f12.6,1X)
	write(2,100) MY(I),MX(J),TAB(I,J)
	end do
	end do
	DO I=1,size(TAB_RAFF,dim=1)
	MX_RAFF(I)=I/2.0
	end do
	DO J=1,size(TAB_RAFF,dim=2)
	MY_RAFF(J)=J/2.0
	END DO

	
	!On veut maintenant avoir les vlaeurs de TAB_RAFF en fonction de celles de TAB en !utilisant la routine INTERPOL_BARI
	DO I=1,size(TAB_RAFF,dim=1)
	DO J=1,size(TAB_RAFF,dim=2)
	CALL INTERPOL_BARI((/MX_RAFF(I),MY_RAFF(J)/),TAB,MX,MY,TAB_RAFF(I,J))
	END DO
	END DO
		
	 DO I=1,NI_RAFF
	DO J=1,NJ_RAFF
	write(3,100) MX_RAFF(I),MY_RAFF(J),TAB_RAFF(I,J)
	end do
	end do
	END PROGRAM Raff



	subroutine dab(n,m,f1,f2,f3,f4,rho,feq)

	!parameter (n=100,m=200)
	real(kind=8) f1(0:n,0:m),f2(0:n,0:m),f3(0:n,0:m),f4(0:n,0:m)
	real(kind=8) rho(0:n,0:m),feq
	!integer i,j,kk

	dx=1.
	dy=dx
	dt=1.

	csq=dx*dx/(dt*dt)
	alpha=0.25
	omega=1.0/(2.*alpha/(dt*csq)+0.5)
	mstep=400
	do j=0,m
	do i=0,n
	rho(i,j)=0.0
	end do
	end do

	do j=0,m
	do i=0,n
	f1(i,j)=0.25*rho(i,j)
	f2(i,j)=0.25*rho(i,j)
	f3(i,j)=0.25*rho(i,j)
	f4(i,j)=0.25*rho(i,j)
	end do
	end do


	do kk=1, mstep
	do j=0,m
	do i=0,n
	feq=0.25*rho(i,j)
	f1(i,j)=omega*feq+(1.-omega)*f1(i,j)
	f2(i,j)=omega*feq+(1.-omega)*f2(i,j)
	f3(i,j)=omega*feq+(1.-omega)*f3(i,j)
	f4(i,j)=omega*feq+(1.-omega)*f4(i,j)
	end do
	end do


	do j=0,m
	do i=1,n
	f1(n-i,j)=f1(n-i-1,j)
	f2(i-1,j)=f2(i,j)
	end do
	end do
	do i=0,n
	do j=1,m
	f3(i,m-j)=f3(i,m-j-1)
	f4(i,j-1)=f4(i,j)
	end do
	end do




	do j=0,m
	f1(0,j)=0.5-f2(0,j)
	f3(0,j)=0.5-f4(0,j)
	f1(n,j)=0.0
	f2(n,j)=0.0
	f3(n,j)=0.0
	f4(n,j)=0.0
	end do

	do i=0,n
	f1(i,m)=0.0
	f2(i,m)=0.0
	f3(i,m)=0.0
	f4(i,m)=0.0
	f1(i,0)=f1(i,1)
	f2(i,0)=f1(i,1)
	f3(i,0)=f1(i,1)
	f4(i,0)=f1(i,1)

	end do
	do j=0,m
	do i=0,n
	rho(i,j)=f1(i,j)+f2(i,j)+f3(i,j)+f4(i,j)
	end do
	end do
	end do
	do j=0,m
	do i=0,n
	write(4,*)i,j,rho(i,j)
	end do
	end do

	return
	end

	subroutine INTERPOL_BARI(X,CHAMP,MAILLX,MAILLY,INTERPOL)
	!==========INTERPOL_BARI===============
	!SUBROUTINE permettant l'interpolation lineaire en un point donne
	!afin de calculer la valeur d'un vecteur en ce meme point.
	!Entree : X, vecteur position, MAILL : maillage, CHAMP : les valeur du champ
	!dont on souhaite connaitre la valeur en X.
	!Sortie : INTERPOL: la valeur interpolee du champ a la position X
	!Marlan, 2011
	implicit none
	real(kind=8),dimension(:),intent(in) :: X, MAILLX, MAILLY
	real(kind=8),dimension(:,:),intent(in) :: CHAMP
	real(kind=8),intent(out)::INTERPOL
	real(kind=8),dimension(size(MAILLX)):: TEMPMAILLX
	real(kind=8),dimension(size(MAILLY)):: TEMPMAILLY
	real(kind=8) :: HI,HJ,INTERPOL1,INTERPOL2
	integer,dimension(2,2):: POSXY

	TEMPMAILLX=MAILLX
	POSXY(1,1)=minloc(abs(X(1)-TEMPMAILLX),dim=1)
	if(minval(abs(X(1)-TEMPMAILLX))/=0.0) TEMPMAILLX(POSXY(1,1))=TEMPMAILLX(POSXY(1,1))+10.0
	POSXY(2,1)=minloc(abs(X(1)-TEMPMAILLX),dim=1)
	TEMPMAILLY=MAILLY
	POSXY(1,2)=minloc(abs(X(2)-TEMPMAILLY),dim=1)
	if(minval(abs(X(2)-TEMPMAILLY))/=0.0) TEMPMAILLY(POSXY(1,2))=TEMPMAILLY(POSXY(1,2))+10.0
	POSXY(2,2)=minloc(abs(X(2)-TEMPMAILLY),dim=1)

	if(minval(abs(X(1)-TEMPMAILLX))/=0.0) HI=-(MAILLX(POSXY(1,1))-X(1))/(MAILLX(POSXY(2,1))-MAILLX(POSXY(1,1)))
	if(minval(abs(X(2)-TEMPMAILLY))/=0.0) HJ=-(MAILLY(POSXY(1,2))-X(2))/(MAILLY(POSXY(2,2))-MAILLY(POSXY(1,2)))
	

	INTERPOL1=(1-HI)*CHAMP(POSXY(1,1),POSXY(1,2))+HI*CHAMP(POSXY(2,1),POSXY(1,2))

	INTERPOL2=(1-HI)*CHAMP(POSXY(1,1),POSXY(2,2))+HI*CHAMP(POSXY(2,1),POSXY(2,2))

	INTERPOL=(1-HJ)*INTERPOL1+HJ*INTERPOL2
	end subroutine INTERPOL_BARI

**abs11** · 13/01/2012, 14h27

SVP est ce que vous pouvez essayez avec des nombres plus grands tel que NI=100 et NJ=50?
ca marche de pour des faibles grandeurs mais pas pour les grands nombres!

Interpolation linéaire et grille

Fortran

Discussions similaires

Partager

Partager