Numerical resolution of coupled partial differential equations governing natural flow of fluid in a rectangular cavity

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. :slight_smile:

For a Fortran beginner, I bet that’s tremendous help already.