pstag.f90

This program incorporates the FV method for solving the Navier-Stokes equations using 2D, Cartesian

!C############################################################
PROGRAM COMET
!C############################################################
!C     This program incorporates the FV method for solving the
!C     Navier-Stokes equations using 2D, Cartesian grids and the
!C     staggered arrangement of variables. Variables are stored
!C     as 2D arrays. SIMPLE method is used for pressure
!C     calculation. UDS and CDS are implemented for the 
!C     discretization of convective terms, CDS is used for the
!C     diffusive terms. The boundary conditions are set for the
!C     lid-driven cavity flow. Only steady flows are considered.
!C
!C                    M. Peric, IfS, Hamburg, 1996
!C#############################################################
	PARAMETER (NX=101,NY=101,NXY=NX*NY)
	COMMON /IDAT/ NI,NJ,NIM,NJM,NIMM,NJMM,IU,IV,IPP,			&
	&       MAXIT,IMON,JMON,IPREF,JPREF,NSWP(3)
	COMMON /RDAT/ RESOR(3),SOR(3),URF(3),SLARGE,SORMAX,ULID,	&
	&       SMALL,GREAT,G(3),ALFA,VIS,DEN,X(NX),Y(NY),			&
	&       XC(NX),YC(NY),R(NY),FX(NX),FY(NY)
	COMMON /COEF/ AE(NX,NY),AW(NX,NY),AN(NX,NY),AS(NX,NY),	&
	&       AP(NX,NY),SU(NX,NY),DU(NX,NY),DV(NX,NY),PP(NX,NY)
	COMMON /VAR/ U(NX,NY),V(NX,NY),P(NX,NY),CX(NX,NY),CY(NX,NY) 
	COMMON /LOGIC/ LCAL(3),LREAD,LWRITE,LTEST,LAKSI
	LOGICAL LCAL,LREAD,LWRITE,LTEST,LAKSI
	CHARACTER*10 FILIN,FILOUT,FILRES,FILGR,FILPL
!C 
!C.....READ FILE NAMES, OPEN FILES 
!C 

	FILIN = 'pstag.inp'
	FILOUT = 'pstag.out'
	FILGR = 'grid.out'
	FILRES = 'pstag.res'
	FILPL = 'pstag.plot'
	
	1 FORMAT(A10)
	OPEN(UNIT=5,FILE=FILIN)
	OPEN(UNIT=1,FILE=FILGR)
	OPEN(UNIT=6,FILE=FILOUT)
	OPEN(UNIT=2,FILE=FILPL,FORM='UNFORMATTED')
	OPEN(UNIT=3,FILE=FILRES,FORM='UNFORMATTED')
	REWIND 1
	REWIND 3
	REWIND 5
	REWIND 6
	
!C
!C.....FORMAT SPECIFICATIONS
!C
	60 FORMAT(20X,'TIME STEP NO.',I3,' ,  TIME =',E10.4,/,20X,36('*'),/) 
	66 FORMAT(2X,'ITER   I-----------ABSOLUTE RESIDUAL SOURCE SUMS----', &
	&      '-------I    I----FIELD VALUES AT MONITORING LOCATION (',	&
	&      I3,',',I3,')---I',/,2X,'NO. ',6X,'UMOM',6X,'VMOM',6X,		&
	&      'MASS',30X,'U',9X,'V',9X,'P',9X,'T',/)
	76 FORMAT(2X,I4,3X,1P3E10.3,24X,1P3E10.3)
	
!C
!C.....INITIALIZATION, GRID DEFINITION, BOUNDARY CONDITIONS
!C
	CALL INIT 
!C
	ITER=0
!C
!C.....READ RESULTS OF PREVIOUS RUN IF CONTINUATION
!C
	IF(LREAD) READ(3) ITER,TIME,NI,NJ,NIM,NJM,NIJ,(X(I),I=1,NI),	&
	&          (Y(J),J=1,NJ),((CX(I,J),J=1,NJ),I=1,NI),			&
	&          ((CY(I,J),J=1,NJ),I=1,NI),((U(I,J),J=1,NJ),I=1,NI),	&
	&          ((V(I,J),J=1,NJ),I=1,NI),((P(I,J),J=1,NJ),I=1,NI)  
	REWIND 3
!C
!C.....INITIAL PRINTOUT
!C
	CALL OUT1 
	IF (LTEST) THEN 
		IF (LCAL(IU)) then
			CALL PRINT(U,'U VEL.') 
		end if
		IF (LCAL(IV)) then 
			CALL PRINT(V,'V VEL.') 
		end if
		IF (LCAL(IPP)) then 
			CALL PRINT(P,'PRESS.')
		end if
	END IF
	WRITE(6,66) IMON,JMON 
!C 
!C.....OUTER ITERATIONS 
!C 
	ITER = 1

	SOURCE = 1
	DO while (ITER .LT. MAXIT .AND. SOURCE.GT.SORMAX)
		IF (LCAL(IU)) then
			CALL CALCU 
		end if
		IF (LCAL(IV)) then 
			CALL CALCV 
		end if
		IF (LCAL(IPP)) then 
			CALL CALCP
		end if
		WRITE(6,76) ITER,RESOR(IU),RESOR(IV),RESOR(IPP),U(IMON,JMON),V(IMON,JMON),P(IMON,JMON) 
		SOURCE=MAX(RESOR(IU),RESOR(IV),RESOR(IPP))
		if (SOURCE.GT.SLARGE) then
			!C
			!C==============================================================
			!C.....OUTER ITERATIONS DIVERGING
			!C==============================================================
			!C
			PRINT *,'     *** TERMINATED - OUTER ITERATIONS DIVERGING ***'
			STOP
		end if
		
		ITER = ITER+1
		print *,ITER
	END DO
!C
!C.....OUTER ITERATIONS CONVERGED -> SAVE AND PRINT
!C
	
	NIJ=NI*NJ
	IF(LCAL(IU)) CALL PRINT(U,'U VEL.') 
	IF(LCAL(IV)) CALL PRINT(V,'V VEL.') 
	IF(LCAL(IPP)) CALL PRINT(P,'PRESS.')
!C
	IF(LWRITE) WRITE(3) ITER,TIME,NI,NJ,NIM,NJM,NIJ,(X(I),I=1,NI),	&
	&         (Y(J),J=1,NJ),((CX(I,J),J=1,NJ),I=1,NI),				&
	&         ((CY(I,J),J=1,NJ),I=1,NI),((U(I,J),J=1,NJ),I=1,NI),		&
	&         ((V(I,J),J=1,NJ),I=1,NI),((P(I,J),J=1,NJ),I=1,NI)		
	REWIND 3
!C
!C==============================================================
!C.....PREPARE AND SAVE DATA FOR PLOT PROGRAM (U -> AE, V -> AW)
!C==============================================================
!C     Interpolate velocities for scalar CV centers, as needed in
!C     the plot program
!C
	itim=1.
	time=0.

	AE(1,1)=U(1,1)
	AW(1,1)=V(1,1)
	AE(NI,1)=U(NIM,1)
	AW(NI,1)=V(NI,1)
!C
	DO J=2,NJM
		AE(1,J)=U(1,J)
		AW(1,J)=0.5*(V(1,J)+V(1,J-1))
		AE(NI,J)=U(NIM,J)
		AW(NI,J)=0.5*(V(NI,J)+V(NI,J-1))
	END DO
!C
	AE(1,NJ)=U(1,NJ)
	AW(1,NJ)=V(1,NJM)
	AE(NI,NJ)=U(NIM,NJ)
	AW(NI,NJ)=V(NI,NJM)
!C
	DO I=2,NIM
		AE(I,1)=0.5*(U(I,1)+U(I-1,1))
		AW(I,J)=V(I,J)
		AE(I,NJ)=0.5*(U(I,NJ)+U(I-1,NJ))
		AW(I,NJ)=V(I,NJM)
	END DO
!C
	DO I=2,NIM
		DO J=2,NJM
			AE(I,J)=0.5*(U(I,J)+U(I-1,J))
			AW(I,J)=0.5*(V(I,J)+V(I,J-1))
		END DO
	END DO
!C
	WRITE(2) ITIM,TIME,NI,NJ,NIM,NJM,NIJ,							&
	&        ((X(I),J=1,NJ),I=1,NI),((Y(J),J=1,NJ),I=1,NI),		&
	&        ((XC(I),J=1,NJ),I=1,NI),((YC(J),J=1,NJ),I=1,NI),		&
	&        ((CX(I,J),J=1,NJ),I=1,NI),((CY(I,J),J=1,NJ),I=1,NI),	&
	&        ((AE(I,J),J=1,NJ),I=1,NI),((AW(I,J),J=1,NJ),I=1,NI),	&
	&        ((P(I,J),J=1,NJ),I=1,NI),((P(I,J),J=1,NJ),I=1,NI)
	REWIND 4
	
	
	
	
END program
!C
!C
!C##############################################################
SUBROUTINE CALCU
!C##############################################################
!C     This routine assembles and solves U-equation. Along west
!C     and east boundary, there are half-CVs which are not
!C     included in the integration region. CX and CY are the mass 
!C     fluxes across east and north faces of mass CVs (positive
!C     for flow out of CV). The geometrical and other coefficients
!C     are first initialized along west and south boundary; then,
!C     in a loop over CVs, east and north cell faces are treated.
!C==============================================================
	PARAMETER (NX=101,NY=101,NXY=NX*NY)						
	COMMON /IDAT/ NI,NJ,NIM,NJM,NIMM,NJMM,IU,IV,IPP,				&
	&       MAXIT,IMON,JMON,IPREF,JPREF,NSWP(3)
	COMMON /RDAT/ RESOR(3),SOR(3),URF(3),SLARGE,SORMAX,ULID,		&
	&       SMALL,GREAT,G(3),ALFA,VIS,DEN,X(NX),Y(NY),				&
	&       XC(NX),YC(NY),R(NY),FX(NX),FY(NY)
	COMMON /COEF/ AE(NX,NY),AW(NX,NY),AN(NX,NY),AS(NX,NY),		&
	&       AP(NX,NY),SU(NX,NY),DU(NX,NY),DV(NX,NY),PP(NX,NY)
	COMMON /VAR/ U(NX,NY),V(NX,NY),P(NX,NY),CX(NX,NY),CY(NX,NY) 
	COMMON /LOGIC/ LCAL(3),LREAD,LWRITE,LTEST,LAKSI
	DIMENSION DW(NY),CW(NY)
	LOGICAL LCAL,LREAD,LWRITE,LTEST,LAKSI
!C============================================================== 
!C
	URFUR=1./URF(IU)
!C 
!C.....INITIALIZE QUANTITIES ALONG WEST BOUNDARY
!C 
	DXER=1./(X(2)-X(1))
	DO J=2,NJM 
		SW=0.5*(R(J)+R(J-1))*(Y(J)-Y(J-1)) 
		DW(J)=VIS*SW*DXER 
		CW(J)=0.5*(CX(2,J)+CX(1,J)) 
	END DO
!C 
!C.....INITIALIZE QUANTITIES ALONG SOUTH BOUNDARY 
!C 
	DO I=2,NIMM
		DX=XC(I+1)-XC(I)
		DYNR=1./(YC(2)-YC(1))
		SS=DX*R(1)
		DXER=1./(X(I+1)-X(I)) 
		DS=VIS*SS*DYNR
		CS=0.5*(CY(I,1)+CY(I+1,1)) 
!C
!C=========================================================== 
!C.....MAIN LOOP - EAST AND NORTH SIDE 
!C=========================================================== 
!C
		DO J=2,NJM 
!C
!C.....CELL FACE AREA, RECIPROCAL DISTANCE FROM P TO N
!C
			SE=0.5*(R(J)+R(J-1))*(Y(J)-Y(J-1)) 
			SN=DX*R(J)
			DYNR=1./(YC(J+1)-YC(J))
!C
!C.....MASS FLUXES ACROSS U-VEL. CV FACES
!C
			CN=0.5*(CY(I,J)+CY(I+1,J)) 
			CE=0.5*(CX(I,J)+CX(I+1,J)) 
!C
!C.....DIFFUSION COEFFICIENTS
!C
			DN=VIS*SN*DYNR
			DE=VIS*SE*DXER 
!C
!C......PRESSURE SOURCE TERM
!C
			SUP=SE*(P(I,J)-P(I+1,J))
!C 
!C.....EXPLICIT SUMM OF CONVECTIVE FLUXES USING UDS AND CDS 
!C 
			AE1=-MIN(CE,0.) 
			AW1= MAX(CW(J),0.) 
			AN1=-MIN(CN,0.) 
			AS1= MAX(CS,0.) 
!C 
			SUCU=AE1*U(I+1,J)+AW1*U(I-1,J)+AN1*U(I,J+1)+AS1*U(I,J-1)-	&
			&         (AE1+AW1+AN1+AS1)*U(I,J)
			SUCC=0.5*((U(I,J)+U(I-1,J))*CW(J)-(U(I+1,J)+U(I,J))*CE)-	&
			&         CN*(U(I,J+1)*FY(J)+U(I,J)*(1.-FY(J)))+				&
			&         CS*(U(I,J)*FY(J-1)+U(I,J-1)*(1.-FY(J-1)))
!C
!C.....MATRIX COEFFICIENTS (ONLY UDS)
!C
			AE(I,J)=-(AE1+DE)
			AW(I,J)=-(AW1+DW(J)) 
			AN(I,J)=-(AN1+DN) 
			AS(I,J)=-(AS1+DS)
!C
!C.....SOURCE TERM
!C 
			SU(I,J)=SUP+G(IU)*(SUCC-SUCU)
!C
!C.....SAVE COEFFICIENTS OF EAST AND NORTH SIDE SIDE 
!C
			DW(J)=DE
			CW(J)=CE
			CS=CN
			DS=DN
!C
		END DO
	END DO
!C 
!C=============================================================
!C.....FINAL ASSEMBLY (NO BOUNDARY MODIFICATIONS NECESSARY)
!C============================================================= 
!C
	DO I=2,NIMM
		DO J=2,NJM 
			AP(I,J)=-(AE(I,J)+AW(I,J)+AN(I,J)+AS(I,J))*URFUR
			SU(I,J)=SU(I,J)+(1.-URF(IU))*AP(I,J)*U(I,J) 
			DU(I,J)=1./AP(I,J) 
		END DO 
	END DO
!C 
!C.....SOLUTION
!C 
	CALL SIPSOL(U,IU,NIMM,NJM)
!C
	
end subroutine 
!C
!C
!C##############################################################
SUBROUTINE CALCV
!C##############################################################
!C     This routine assembles and solves V-equation. Along south
!C     and north boundary, there are half-CVs which are not
!C     included in the integration region. CX and CY are the mass 
!C     fluxes across east and north faces of mass CVs (positive
!C     for flow out of CV). The geometrical and other coefficients
!C     are first initialized along west and south boundary; then,
!C     in a loop over CVs, east and north cell faces are treated.
!C==============================================================
	PARAMETER (NX=101,NY=101,NXY=NX*NY)
	COMMON /IDAT/ NI,NJ,NIM,NJM,NIMM,NJMM,IU,IV,IPP,				&
	&       MAXIT,IMON,JMON,IPREF,JPREF,NSWP(3)
	COMMON /RDAT/ RESOR(3),SOR(3),URF(3),SLARGE,SORMAX,ULID,		&
	&       SMALL,GREAT,G(3),ALFA,VIS,DEN,X(NX),Y(NY),				&
	&       XC(NX),YC(NY),R(NY),FX(NX),FY(NY)
	COMMON /COEF/ AE(NX,NY),AW(NX,NY),AN(NX,NY),AS(NX,NY),		&
	&       AP(NX,NY),SU(NX,NY),DU(NX,NY),DV(NX,NY),PP(NX,NY)
	COMMON /VAR/ U(NX,NY),V(NX,NY),P(NX,NY),CX(NX,NY),CY(NX,NY) 
	COMMON /LOGIC/ LCAL(3),LREAD,LWRITE,LTEST,LAKSI
	DIMENSION DW(NY),CW(NY)
	LOGICAL LCAL,LREAD,LWRITE,LTEST,LAKSI
!C===============================================================
!C
	URFVR=1./URF(IV)
!C 
!C.....INITIALIZE QUANTITIES ALONG WEST BOUNDARY
!C 
	DXER=1./(XC(2)-XC(1))
	DO J=2,NJMM
		SW=0.5*(YC(J+1)-YC(J))*(R(J+1)+R(J-1))
		DW(J)=VIS*SW*DXER 
		CW(J)=0.5*(CX(1,J)+CX(1,J+1))
	END DO
!C 
!C.....INITIALIZE QUANTITIES ALONG SOUTH BOUNDARY 
!C 
	DO I=2,NIM 
		DX=X(I)-X(I-1) 
		DXER=1./(XC(I+1)-XC(I))
		DYNR=1./(Y(2)-Y(1))
		SS=0.5*(R(2)+R(1))*DX
		DS=VIS*SS*DYNR
		CS=0.5*(CY(I,2)+CY(I,1))
!C 
!C==============================================================
!C.....MAIN LOOP - EAST AND NORTH SIDES
!C============================================================== 
!C
		DO J=2,NJMM
!C
!C.....CELL FACE AREA, RECIPROCAL DISTANCE FROM P TO N
!C
			DY=YC(J+1)-YC(J)
			SE=DY*0.5*(R(J+1)+R(J-1))
			SN=0.5*(R(J+1)+R(J))*DX
			DYNR=1./(Y(J+1)-Y(J)) 
!C
!C.....MASS FLUXES THROUGH EAST AND NORTH FACE OF V-CV
!C
			CN=0.5*(CY(I,J+1)+CY(I,J)) 
			CE=0.5*(CX(I,J+1)+CX(I,J)) 
!C
!C......DIFFUSION COEFFICIENTS
!C
			DN=VIS*SN*DYNR 
			DE=VIS*SE*DXER
!C
!C.....PRESSURE SOURCE TERM
!C
			SUP=DX*R(J)*(P(I,J)-P(I,J+1))
!C 
!C.....EXPLICIT SUM OF CONVECTIVE FLUXES USING UDS AND CDS
!C 
			AE1=-MIN(CE,0.) 
			AW1= MAX(CW(J),0.) 
			AN1=-MIN(CN,0.) 
			AS1= MAX(CS,0.) 
!C 
			SUCU=AE1*V(I+1,J)+AW1*V(I-1,J)+AN1*V(I,J+1)+AS1*V(I,J-1)-	&
			&         (AE1+AW1+AN1+AS1)*V(I,J)
			SUCC=0.5*((V(I,J)+V(I,J-1))*CS-(V(I,J+1)+V(I,J))*CN)-		&
			&         CE*(V(I+1,J)*FX(I)+V(I,J)*(1.-FX(I)))+				&
			&         CW(J)*(V(I,J)*FX(I-1)+V(I-1,J)*(1.-FX(I-1)))
!C
!C.....MATRIX COEFFICIENTS (ONLY UDS)
!C
			AE(I,J)=-(AE1+DE)