Hello Dear all, I create this topic to have a help about Numerical Simulations of coupled partial differential equations governing natural convection in a rectangular cavity filled with fluid and heated from below. I need one exemple of Fortran code using finite volume method or Lattice Boltzmann method or Differential Quadrature method or any other method. In the following picture, I presented the Equations (1) to (9) that I want to solve.
Hi Justin, is this a homework problem? Please tag it as homework if so.
What does your current solution look like? What did you try and where did you get stuck? If you provide as much detail as possible, people here will be able to help with exactly the problem that you have.
This is a problem that I have been instructed to solve in order to familiarize myself with the numerical solution of similar equations. I did the discretization by the differential quadrature method but I have difficulties writing a Fortran code.
Did you write, compile, and run any Fortran code before? If not, Learn Fortran - Fortran Programming Language should get you going with the basics. Give it a try and please provide information about a specific thing that you get stuck with.
The code I wrote for the sample case when Ha = 0 is the following, but when I compile it by taking:
- N=21
- M=21
- NTT=1
- H=0.001
- EPS1=0.000001
- EPS2=0.000001
- Pr=0.71
- Ra=1000
- ALAY=0.66
the error is very large and I can’t get the exact results. It would be difficult to explain the discretization method and the various parameters used as well as the acronyms, which is why I am asking if anyone has any other sample code of any method to provide me with detail.
Program pdqnc
implicit real*8(A-H,O-Z)
dimension AX(50,50),BX(50,50),XX(50),AY(50,50),BY(50,50),YY(50)
dimension U(50,50),V(50,50),PD(401,401),P(50,50)
dimension T(50,50),W(50,50),DT(50,50),DW(50,50)
! Imput number of mesh points, maximum iteration number, time step
! size, convergence criteria, prandtl number and rayleigh number
write(6,*) 'N,M,NTT,H,EPS1,EPS2,PR,RA'
read(5,*) N,M,NTT,H,EPS1,EPS2,PR,RA
! Mesh generation and set-up of initial flow field
call GRID(N,M,XX,YY,U,V,W,T,DT,DW)
! Computation of weighting coefficient of the first and second order
! derivative in the x and y direction.
call WCOEF(1,N,XX,AX,BX)
call WCOEF(1,M,YY,AY,BY)
! Form the matrix of equation system 6.21. The obtained matrix is
! then decomposed by LU decomposition (see equation 5.27).
call PDCO(N,M,AX,BX,AY,BY,PD)
NT=1
! Start iteration
88 continue
! Implement boundary conditions for vorticity and temperature
call BOCO(N,M,U,V,W,T,AX,AY)
!4-Stage Runge-Kutta method to solve resultant ordinary diffrencial
! equationqs for vortivity and temperature
call RUNGU(N,M,H,RA,PR,AX,BX,AY,BY,U,V,T,W,DT,DW)
! Compute stream function
call COREP(N,M,W,P,PD)
!Update stream function on the boundary
call BPCO(N,M,AX,AY,P)
!Compute velocity U,V
call CORUV(N,M,AX,AY,U,V,P)
!Find out maximum residuals for equations 6.22 and 6.23
ERR1=DT(1,1)**2
ERR2=DW(1,1)**2
do 91 I=1,N
do 91 J=1,M
ERR1=DT(I,J)**2+ERR1
ERR2=DW(I,J)**2+ERR2
91 continue
NM=N*M
ERR1=sqrt(ERR1/NM)
ERR2=sqrt(ERR2/NM)
write(6,113) NT,ERR1,ERR2
113 format(1X,I10,5X,' ERR1= ',E10.5,5X,' ERR2= ',E10.5)
! Check conergence criteria
if (ERR1.le.EPS1.and.ERR2.le.EPS2) go to 92
if (NT.ge.NTT) go to 92
NT=NT+1
go to 88
92 continue
! Compute nusselt numbers
call NUSEC(N,M,XX,YY,AX,U,T)
! Interpolate numerical results from coarse mesh to fine mesh for
! plotting and other purpose
call OUTPLT(N,M,XX,YY,U,V,W,T,P)
stop
end program pdqnc
! This subroutine is to compute the PDQ weighting coefficients of
! the first and second order derivatives when the coordinates of
! mesh points are input
subroutine WCOEF(N0,NN,XX,AX,BX)
implicit real*8(A-H,O-Z)
dimension AX(50,50),BX(50,50),XX(50)
dimension AIX(50)
NC=NN-N0+1
! Compute the weighting coefficients of the first order derivatives
do 72 I=N0,NN
II=I-N0+1
AX1=1.
do 73 J=N0,NN
if (J.eq.I) go to 73
AX1=AX1*(XX(I)-XX(J))
73 continue
72 AIX(II)=AX1
do 74 I=1,NC
II=I+N0-1
do 74 J=1,NC
JJ=J+N0-1
if (J.eq.I) go to 74
AX(I,J)=AIX(I)/((XX(II)-XX(JJ))*AIX(J))
74 continue
do 77 I=1,NC
AXI=0.
do 78 J=1,NC
if (J.eq.I) go to 78
AXI=AXI+AX(I,J)
78 continue
AX(I,I)=-AXI
77 continue
! Compute the weighting of the second order derivate
do 300 I=1,NC
II=I+N0-1
do 300 J=1,NC
JJ=J+N0-1
if (J.eq.I) go to 300
BX(I,J)=2.*AX(I,J)*(AX(I,I)-1./(XX(II)-XX(JJ)))
300 continue
do 400 I=1,NC
BXI=0.
do 410 J=1,NC
if(J.eq.I) go to 410
BXI=BXI+BX(I,J)
410 continue
BX(I,I)=-BXI
400 continue
return
end
! This subroutine is to generate the coordinates of mesh points
subroutine GRID(N,M,XX,YY,U,V,W,T,DT,DW)
implicit real*8(A-H,O-Z)
dimension XX(50),YY(50),U(50,50),V(50,50)
dimension T(50,50),W(50,50),DT(50,50),DW(50,50)
write(6,*) 'ALAX,ALAY'
read(5,*) ALAX,ALAY
PI=4.*ATAN(1.)
! Lx=1.
! Ly=1.
do 71 I=1,N
ALFA=PI*FLOAT(I-1)/FLOAT(N-1)
! ALFA1=(cos(PI/(2.*N))-cos((((2.*I)-1.)*PI)/(2.*N)))
! ALFA2=(cos(PI/(2.*N))-cos((((2.*N)-1.)*PI)/(2.*N)))
71 XX(I)=0.5*(1.-cos(ALFA))
! 71 XX(I)=(ALFA1/ALFA2)*Lx
X11=XX(1)
! Stretch the mesh towards the boundary, ALAX.lt.1
do 76 I=1,N
76 XX(I)=(1.-ALAX)*(3.*PI*XX(I)**2-2.*XX(I)**3)/PI**2+ALAX*XX(I)
do 81 J=1,M
ALFA=PI*FLOAT(J-1)/FLOAT(M-1)
! ALFA1=(cos(PI/(2.*M))-cos((((2.*J)-1.)*PI)/(2.*M)))
! ALFA2=(cos(PI/(2.*M))-cos((((2.*M)-1.)*PI)/(2.*M)))
81 YY(J)=0.5*(1.-cos(ALFA))
! 81 YY(J)=(ALFA1/ALFA2)*Ly
! Stretch the mesh towards the boundary, ALAY.lt.1
do 86 J=1,M
86 YY(J)=(1.-ALAY)*(3.*PI*YY(J)**2-2.*YY(J)**3)/PI**2+ALAY*YY(J)
YY(MH)=0.5
! Set up initial flow field
do 18 J=1,M
do 10 I=1,N
U(I,J)=0.
V(I,J)=0.
T(I,J)=0.
W(I,J)=0.
DT(I,J)=0.
DW(I,J)=0.
10 continue
T(1,J)=1.0
18 continue
return
end
! Implement boudary conditions for vorticity and temperature
subroutine BOCO(N,M,U,V,W,T,AX,AY)
implicit real*8(A-H,O-Z)
dimension AX(50,50),AY(50,50),U(50,50),V(50,50)
dimension T(50,50),W(50,50)
! Update vorticity on the Bottom and top boundaries
do 72 I=1,N
SUM1=0.
SUM2=0.
do 10 J=1,M
SUM1=SUM1+AY(1,J)*U(I,J)
10 SUM2=SUM2+AY(M,J)*U(I,J)
W(I,1)=SUM1
72 W(I,M)=SUM2
! Update vorticity on the left and right boundaries
do 73 J=1,M
SUM3=0.
SUM4=0.
do 20 I=1,N
SUM3=SUM3+AX(1,I)*V(I,J)
20 SUM4=SUM4+AX(N,I)*V(I,J)
W(1,J)=-SUM3
W(N,J)=-SUM4
! Temperature condition at left and right boundaries
T(1,J)=1.
T(N,J)=0.
73 continue
return
end
! 4-Runge-Kutta method to solve ordinary differential equations
! for vorticity and temperature. The scheme is shown in equation 5.15
subroutine RUNGU(N,M,H,RA,PR,AX,BX,AY,BY,U,V,T,W,DT,DW)
implicit real*8(A-H,O-Z)
dimension AX(50,50),BX(50,50),AY(50,50),BY(50,50)
dimension U(50,50),V(50,50),T(50,50),W(50,50)
dimension CT(50,50),CW(50,50),DT(50,50),DW(50,50)
! Stage 1
call TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,T,W,DW)
call TRFT(N,M,AX,BX,AY,BY,U,V,T,DT)
do 3 I=1,N
do 3 J=1,M
CW(I,J)=W(I,J)+0.25*H*DW(I,J)
3 CT(I,J)=T(I,J)+0.25*H*DT(I,J)
! Stage 2
call TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,CT,CW,DW)
call TRFT(N,M,AX,BX,AY,BY,U,V,CT,DT)
do 5 I=1,N
do 5 J=1,M
CW(I,J)=W(I,J)+H*DW(I,J)/3.
5 CT(I,J)=T(I,J)+H*DT(I,J)/3.
! Stage 3
call TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,CT,CW,DW)
call TRFT(N,M,AX,BX,AY,BY,U,V,CT,DT)
do 7 I=1,N
do 7 J=1,M
CW(I,J)=W(I,J)+H*DW(I,J)/2.
7 CT(I,J)=T(I,J)+H*DT(I,J)/2.
! Stage 4
call TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,CT,CW,DW)
call TRFT(N,M,AX,BX,AY,BY,U,V,CT,DT)
do 9 I=2,N-1
do 9 J=2,M-1
W(I,J)=W(I,J)+H*DW(I,J)
9 T(I,J)=T(I,J)+H*DT(I,J)
call TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,T,W,DW)
call TRFT(N,M,AX,BX,AY,BY,U,V,T,DT)
return
end
! Compute the residual of vorticity equation 6.22
subroutine TRFW(N,M,RA,PR,AX,BX,AY,BY,U,V,CT,CW,DW)
implicit real*8(A-H,O-Z)
dimension AX(50,50),BX(50,50),AY(50,50),BY(50,50)
dimension U(50,50),V(50,50),CT(50,50),CW(50,50),DW(50,50)
do 10 I=2,N-1
do 10 J=2,M-1
DWI=0.
DTI=0.
DDWI=0.
do 20 K=1,N
DDWI=DDWI+BX(I,K)*CW(K,J)
DTI=DTI+AX(I,K)*CT(K,J)
20 DWI=DWI+AX(I,K)*CW(K,J)
DWJ=0.
DDWJ=0.
do 30 K=1,M
DWJ=DWJ+AY(J,K)*CW(I,K)
30 DDWJ=DDWJ+BY(J,K)*CW(I,K)
DW(I,J)=(DDWI+DDWJ)*PR-U(I,J)*DWI-V(I,J)*DWJ-DTI*PR*RA
10 continue
return
end
! Compute the residual of temperature equation 6.23
subroutine TRFT(N,M,AX,BX,AY,BY,U,V,CT,DT)
implicit real*8(A-H,O-Z)
dimension AX(50,50),BX(50,50),AY(50,50),BY(50,50)
dimension U(50,50),V(50,50),CT(50,50),DT(50,50)
! Update the temperature on the top and bottom boundaries using
! equations 6.42b and 6.42c
do 70 I=2,N-1
C2=0.
C3=0.
do 80 J=2,M-1
C2=C2+(AY(1,J)*AY(M,M)-AY(M,J)*AY(1,M))*CT(I,J)
C3=C3+(AY(1,J)*AY(M,1)-AY(M,J)*AY(1,1))*CT(I,J)
80 continue
CT(I,1)=-C2/(AY(1,1)*AY(M,M)-AY(M,1)*AY(1,M))
70 CT(I,M)=-C3/(AY(1,M)*AY(M,1)-AY(M,M)*AY(1,1))
! Compute residuals of temperature equation 6.23 for all interior
! points
do 10 I=2,N-1
do 10 J=2,M-1
DTI=0.
DDTI=0.
do 20 K=1,N
DDTI=DDTI+BX(I,K)*CT(K,J)
20 DTI=DTI+AX(I,K)*CT(K,J)
DTJ=0.
DDTJ=0.
do 30 K=1,M
DTJ= DTJ+AY(J,K)*CT(I,K)
30 DDTJ=DDTJ+BY(J,K)*CT(I,K)
DT(I,J)=DDTI+DDTJ-U(I,J)*DTI-V(I,J)*DTJ
10 continue
return
end
! Form the matrix of equation system 6.21. the two layer approach
! is used to implement the boundary condition for stream function.
! The obtained matrix is decomposed by LU decomposition
subroutine PDCO(N,M,AX,BX,AY,BY,PD)
implicit real*8(A-H,O-Z)
dimension AX(50,50),BX(50,50),AY(50,50),BY(50,50)
dimension PD(401,401),PL(401,401),PU(401,401)
NN=N-4
MM=M-4
NM=NN*MM
N1=N-1
M1=M-1
! AXN and AYM are used in equation 6.30 and 6.31 for two-layer
! approach
AXN=AX(1,2)*AX(N,N1)-AX(N,2)*AX(1,N1)
AYM=AY(1,2)*AY(M,M1)-AY(M,2)*AY(1,M1)
do 10 I=3,N-2
do 10 J=3,M-2
IJ=MM*(I-3)+J-2
do 20 K=1,NM
20 PD(IJ,K)=0.
do 30 K1=3,N-2
KJ=MM*(K1-3)+J-2
! Equation 6.30 is directly substituted into equation 6.21
AXK1=AX(1,K1)*AX(N,N1)-AX(N,K1)*AX(1,N1)
AXKN=AX(N,K1)*AX(1,2)-AX(1,K1)*AX(N,2)
30 PD(IJ,KJ)=BX(I,K1)-(AXK1*BX(I,2)+AXKN*BX(I,N1))/AXN
do 40 K2=3,M-2
IK=MM*(I-3)+K2-2
! Equation 6.31 is directly substituted into equation 6.21
AYK1=AY(1,K2)*AY(M,M1)-AY(M,K2)*AY(1,M1)
AYKM=AY(M,K2)*AY(1,2)-AY(1,K2)*AY(M,2)
40 PD(IJ,IK)=PD(IJ,IK)+BY(J,K2)-(AYK1*BY(J,2)+AYKM*BY(J,M1))/AYM
10 continue
!LU decomposition using equation 5.27
do 50 I=1,NM
do 50 J=1,NM
if(J.ge.I) go to 51
J1=J-1
SUM1=0.
if(J1.eq.0) go to 57
do 55 JJ=1,J1
55 SUM1=SUM1+PL(I,JJ)*PU(JJ,J)
57 PL(I,J)=(PD(I,J)-SUM1)/PU(J,J)
go to 50
51 if(I.eq.1) go to 52
I1=I-1
SUM2=0.
do 53 II=1,I1
53 SUM2=SUM2+PL(I,II)*PU(II,J)
PU(I,J)=PD(I,J)-SUM2
go to 50
52 PU(I,J)=PD(I,J)
50 continue
! The tringular matrices [PL] and [PU] are stored in [PD]
do 60 I=1,NM
do 60 J=1,NM
if(J.ge.I) go to 61
PD(I,J)=PL(I,J)
go to 60
61 PD(I,J)=PU(I,J)
60 continue
return
end
! The algebric equation system for stream function is solved by LU
! decomposition
subroutine COREP(N,M,W,P,PD)
implicit real*8(A-H,O-Z)
dimension PD(401,401),P(50,50),W(50,50)
dimension B(401),BI(401),BJ(401)
NN=N-4
MM=M-4
NM=NN*MM
do 10 I=3,N-2
do 10 J=3,M-2
IJ=MM*(I-3)+J-2
B(IJ)=W(I,J)
10 continue
! Equation 5.26a is used to get the solution at the first step
do 40 I=1,NM
if(I.eq.1) go to 41
I1=I-1
SUM1=0.
do 42 J=1,I1
42 SUM1=SUM1+PD(I,J)*BI(J)
BI(I)=B(I)-SUM1
go to 40
41 BI(1)=B(1)
40 continue
! Equation 5.26b is used to get the solution at the second step
do 50 I=1,NM
IV=NM-I+1
if(I.eq.1) GO TO 51
I2=I-1
SUM2=0.
do 52 J=1,I2
IVJ=NM-J+1
52 SUM2=SUM2+PD(IV,IVJ)*BJ(IVJ)
BJ(IV)=(BI(IV)-SUM2)/PD(IV,IV)
go to 50
51 BJ(IV)=BI(IV)/PD(IV,IV)
50 continue
do 60 I=3,N-2
do 60 J=3,M-2
IJ=MM*(I-3)+J-2
60 P(I,J)=BJ(IJ)
return
end
! Implement boundary condition for stream function. The two-layer
! approach is used
subroutine BPCO(N,M,AX,AY,P)
implicit real*8(A-H,O-Z)
dimension AX(50,50),AY(50,50),P(50,50)
N1=N-1
M1=M-1
! AXN and AYM are used in equations 6.30 and 6.31 for two-layer
!approach
AXN=AX(1,2)*AX(N,N1)-AX(N,2)*AX(1,N1)
AYM=AY(1,2)*AY(M,M1)-AY(M,2)*AY(1,M1)
! Equation 6.31 is used to update the stream function
do 20 I=3,N-2
SUM3=0.
SUM4=0.
do 21 J=3,M-2
SUM3=SUM3-(AY(1,J)*AY(M,M1)-AY(M,J)*AY(1,M1))*P(I,J)
21 SUM4=SUM4-(AY(M,J)*AY(1,2)-AY(1,J)*AY(M,2))*P(I,J)
P(I,2)=SUM3/AYM
20 P(I,M1)=SUM4/AYM
! Equation 6.30 is used to update the stream function
do 10 J=2,M-1
SUM1=0.
SUM2=0.
do 11 I=3,N-2
SUM1=SUM1-(AX(1,I)*AX(N,N1)-AX(N,I)*AX(1,N1))*P(I,J)
11 SUM2=SUM2-(AX(N,I)*AX(1,2)-AX(1,I)*AX(N,2))*P(I,J)
P(2,J)=SUM1/AXN
10 P(N1,J)=SUM2/AXN
do 30 J=1,M
P(1,J)=0.
30 P(N,J)=0.
do 40 I=2,N1
P(I,1)=0.
40 P(I,M)=0.
return
end
! Compute velocity U,V from stream fonction
subroutine CORUV(N,M,AX,AY,U,V,P)
implicit real*8(A-H,O-Z)
dimension AX(50,50),AY(50,50),U(50,50),V(50,50),P(50,50)
do 10 I=2,N-1
do 10 J=2,M-1
DPI=0.
! Compute V velocity for interior points
DO 20 k=1,N
20 DPI=DPI+AX(I,K)*P(K,J)
V(I,J)=-DPI
DPJ=0.
! Compute U velocity for interior points
do 30 k=1,M
30 DPJ=DPJ+AY(J,K)*P(I,K)
10 U(I,J)=DPJ
do 40 I=1,N
U(I,1)=0.
V(I,1)=0.
U(I,M)=0.
40 V(I,M)=0.
do 50 J=2,M-1
U(1,J)=0.
V(1,J)=0.
U(N,J)=0.
50 V(N,J)=0.
return
end
! Compute the Nusselt numbers
subroutine NUSEC(N,M,XX,YY,AX,U,T)
implicit real*8(A-H,O-Z)
dimension XX(50),YY(50),X2(201),Y2(201),AX(50,50),U(50,50)
dimension T(50,50),QQ(50,50),QI(50),QT(50,201),CNU(50),CNUT(201)
! Compute the local heat flux in a horizontal direction Q(X,Y)
do 10 I=1,N
do 10 J=1,M
DTI=0.
do 11 K=1,N
11 DTI=DTI+AX(I,K)*T(K,J)
10 QQ(I,J)=U(I,J)*T(I,J)-DTI
! 10 QQ(I,J)=DTI
! Set up a fine mesh
NN=101
MM=101
TDX=(XX(N)-XX(1))/FLOAT(NN-1)
TDY=(YY(M)-YY(1))/FLOAT(MM-1)
do 20 I=1,NN
20 X2(I)=FLOAT(I-1)*TDX
do 30 J=1,MM
30 Y2(J)=FLOAT(J-1)*TDY
! Interpolate Q(X,Y) to the fine mesh in the Y direction
do 40 I=1,N
do 41 J=1,M
QI(J)=QQ(I,J)
41 continue
do 42 J=1,MM
YJ=Y2(J)
call CZ(M,YJ,YY,QI,QOUT)
QT(I,J)=QOUT
42 continue
40 continue
! Compute average Nusselt number on the vertical line
do 188 I=1,N
SUM=0.0
do 198 J=2,MM
SUM=SUM+0.5*(QT(I,J)+QT(I,J-1))*TDY
198 continue
CNU(I)=SUM
188 continue
do 52 I=1,NN
XI=X2(I)
call CZ(N,XI,XX,CNU,CNUOUT)
CNUT(I)=CNUOUT
52 continue
! Compute average Nusselt number throughout the cavity
CNUAVG=0.0
do 918 I=2,NN
CNUAVG=CNUAVG+0.5*(CNUT(I)+CNUT(I-1))*TDX
918 continue
! Find out the maximum and minimum local Nusselt numbers at X=0
! together with their positions
AMIN=QT(1,1)
AMAX=QT(1,1)
do 518 J=1,MM
if(QT(1,J).le.AMAX) go to 888
AMAX=QT(1,J)
YMAX=Y2(J)
888 if (QT(1,J).ge.AMIN) go to 518
AMIN=QT(1,J)
YMAX=Y2(J)
518 continue
NH=(N+1)/2
! OUTPUT THE Nusselt Numbers
write(6,70) CNU(1)
70 FORMAT(2X,'AVERAGE NUSSELT NUMBER AT X=0',5X,F15.4)
write(6,72) CNU(NH)
72 FORMAT(2X,'AVERAGE NUSSELT NUMBER AT X=0.5',3X,F15.4)
write(6,74) CNUAVG
74 FORMAT(2X,'AVERAGE NUSSELT NUMBER OVER CAVITY',F15.4)
write(6,76) AMAX,YMAX
76 FORMAT(2X,'NUMAX AND POSITION AT X=0',5X,F10.4,5X,F10.4)
write(6,78) AMIN,YMIN
78 FORMAT(2X,'NUMIN AND POSITION AT X=0',5X,F10.4,5X,F10.4)
return
end
! Lagrange interpolation
subroutine CZ(N,XX,X,F,FV)
implicit real*8(A-H,O-Z)
dimension X(50),F(50)
FV=0.
do 10 I=1,N
AD=1.
do 20 J=1,N
if(J.eq.I) go to 21
AD1=(XX-X(J))/(X(I)-X(J))
go to 22
21 AD1=1.
22 AD=AD*AD1
20 continue
FV=FV+AD*F(I)
10 continue
return
end
! Interpolation the numerical results from a coarse mesh to a fine
! mesh
subroutine OUTPLT(N,M,X1,Y1,U,V,W,G,P)
implicit real*8(A-H,O-Z)
dimension X1(50),Y1(50),X2(101),Y2(101),P(50,50),W(50,50)
dimension U(50,50),V(50,50),G(50,50),UT(101,101),VT(101,101)
dimension GT(101,101),PT(101,101),WT(101,101),UI(51),UJ(51)
dimension VI(51),VJ(51),WI(51),WJ(51),PI(51),PJ(51),GI(50),GJ(50)
X2(1)=X1(1)
Y2(1)=Y1(1)
write(6,*) 'Read the values of NN, MM for interpolation'
read(5,*) NN,MM
TDX=(X1(N)-X1(1))/FLOAT(NN-1)
TDY=(Y1(M)-Y1(1))/FLOAT(MM-1)
! Generate fine mesh
do 20 I=2,NN
20 X2(I)=X2(I-1)+TDX
do 30 J=2,MM
30 Y2(J)=Y2(J-1)+TDY
! Interpolate in the X direction
do 40 J=1,M
do 41 I=1,N
UI(I)=U(I,J)
VI(I)=V(I,J)
WI(I)=W(I,J)
GI(I)=G(I,J)
PI(I)=P(I,J)
41 continue
do 42 I=1,NN
XI=X2(I)
call CZ(N,XI,X1,UI,UOUT)
call CZ(N,XI,X1,VI,VOUT)
call CZ(N,XI,X1,WI,WOUT)
call CZ(N,XI,X1,GI,GOUT)
call CZ(N,XI,X1,PI,POUT)
UT(I,J)=UOUT
VT(I,J)=VOUT
WT(I,J)=WOUT
GT(I,J)=GOUT
PT(I,J)=POUT
42 continue
40 continue
! Interpolate in the Y direction
do 50 I=1,NN
do 51 J=1,M
UJ(J)=UT(I,J)
VJ(J)=VT(I,J)
WJ(J)=WT(I,J)
GJ(J)=GT(I,J)
PJ(J)=PT(I,J)
51 continue
do 52 J=1,MM
YJ=Y2(J)
call CZ(M,YJ,Y1,UJ,UOUT)
call CZ(M,YJ,Y1,VJ,VOUT)
call CZ(M,YJ,Y1,WJ,WOUT)
call CZ(M,YJ,Y1,GJ,GOUT)
call CZ(M,YJ,Y1,PJ,POUT)
UT(I,J)=UOUT
VT(I,J)=VOUT
WT(I,J)=WOUT
GT(I,J)=GOUT
PT(I,J)=POUT
52 continue
50 continue
! Output data for contour plotting by tecplot (at software)
open(7,file='NSTO.dat')
! Write(7,*) 'TITLE= " " '
! write(7,*) 'VARIABLE=X,Y,U,V,T,P,W'
! write(7,*) 'ZONE I=',NN,'J=',MM,',F=POINT'
do 850 J=1,MM
do 850 I=1,NN
write(7,950) X2(I),Y2(J),UT(I,J),VT(I,J),GT(I,J),PT(I,J),WT(I,J)
850 continue
950 FORMAT(F7.4,1X,F7.4,1X,E11.5,2X,E11.5,2X,E11.5,2X,E11.5,2X,E11.5)
close(7)
return
end
My generic advice is to print out intermediate results of your program to figure out where things are going wrong.
It is better to declare variables explicitly and use IMPLICIT NONE instead of writing
implicit real*8(A-H,O-Z)
The code
dimension AX(50,50),BX(50,50),XX(50),AY(50,50),BY(50,50),YY(50)
dimension U(50,50),V(50,50),PD(401,401),P(50,50)
dimension T(50,50),W(50,50),DT(50,50),DW(50,50)
should be rewritten by defining a PARAMETER, say ngrid equal to 50 and using it in the declarations. Is the 401 related to the 50 in some way? Does the code work if the discretization is made finer by increasing the value of ngrid?
For your code
write(6,*) 'N,M,NTT,H,EPS1,EPS2,PR,RA'
read(5,*) N,M,NTT,H,EPS1,EPS2,PR,RA
can you check that the values are read properly by printing them out? Reading many parameters in a single statement is error-prone.
1 Like
Thanks for your comments. I will use it to review the code and come back. But I can’t check that the values are read properly by printing them out.
The code don’t work if the discretization is made finer by increasing the value of ngrid.
You MUST use explicit typing if you want assistance from me. There is an undefined variable reference in the code you posted. Luckily I have a tool that helps me find it. It is so easy to make a typo and use a new name (whose value is undefined) when you meant to use an existing name with a defined value. With implicit typing, how are you going to find it?
1 Like
I agree with @themos about adding IMPLICIT NONE. Compilng the code with gfortran -Wall, I get
202 | YY(MH)=0.5
|
Warning: 'mh' is used uninitialized [-Wuninitialized]
If a user sees many compiler warnings that can be disregarded, she may miss the important ones. Before the message above, gfortran emits numerous messages such as
797 | dimension GT(101,101),PT(101,101),WT(101,101),UI(51),UJ(51)
| 1
Warning: Array 'gt' at (1) is larger than limit set by '-fmax-stack-var-size=', moved from stack to static storage. This makes the procedure unsafe when called recursively, or concurrently from multiple threads. Consider using '-frecursive', or increase the '-fmax-stack-var-size=' limit, or change the code to use an ALLOCATABLE array. [-Wsurprising]
xcavity.f90:797:32:
Compiling without the -Wall option, there is no warning about uninitialized mh and the only messages are ones such as
62 | do 91 J=1,M
| 1
Warning: Fortran 2018 deleted feature: Shared DO termination label 91 at (1)
xcavity.f90:118:18:And the first advice (“use implicit none”) is free. 
For a Fortran beginner, I bet that’s tremendous help already.