PROGRAM SET4NSTEQ1NEW
C------------- VERTICAL SCAN / ITYPE OPTION ----\\----------------------
C----------- HORIZONTAL DIRECTION: LEFT TO RIGHT -----------------------
C=======================================================================
C DATA GENERATING PROGRAM FOR NSETQ8XX.FOR
C PROJECT: FREE CONVECTION PROJECT NAME: FREECV3
C DOMAIN: SQUARE DOMAIN BOUNDED BY 0'C(LEFT) AND 8'C(RIGHT)
C ELEMENT: FOUR-NODED ISOPARAMETRIC ELEMENT
C DOMAIN DISCRETIZATION: UNEVEN ELEMENTS WITH VERTICAL SCAN
C REYNOLDS NUMBER: SPECIFIED IN PARAMETER AS REYNLD
C VISCOSITY RATIO [VRATIO]: 2ND VISCOSITY / 1ST VISCOSITY
C EIJI FUKUMORI DECEMBER 29, 1993
C FEB. 02, 2013
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=4, MXE=1200, MXN=1300, MXB=200, NF=3 )
PARAMETER ( VRATIO=1.D9, VISCO=0.018D0, FLMDA=VISCO*VRATIO,
* DT=0.1D0,DTMAX=10.D0,TMAX=7200.D0,MAXITE=8,ERMAX=0.0005D0,
* C1=1.2D0, CONDUCT=0.001339D0,TREF=4.D0,ITEFIX=4,
* T1 = 0.D0, T2 = 8.D0, TLX = 3.D0, TLY = 4.D0,
* NEY = 20, NEX = 20, DX = TLX / NEX, DY = TLY / NEY,
* NDX=NEX+1, NDY=NEY+1, TIME = 0.D0, ITYPE=2 )
PARAMETER ( INTEPT=2 )
C=======================================================================
DIMENSION NODEX(MXE,ND), XCOORD(MXN), YCOORD(MXN),
* IBNDFX(MXB),IBNDFY(MXB),IBNDT(MXB),BVX(MXB),BVY(MXB),BVT(MXB),
* U(MXN), V(MXN),TEMP(MXN), IBNDS(MXB), BVS(MXB), NODES(NDY)
CHARACTER FNAME(NF)*19,EXTEN(NF)*4,PROJECT*15,EXFILE*3
LOGICAL YES
C=======================================================================
DATA PROJECT / 'FREECONVECTION1' /
DATA EXTEN / '.DAT', '.BIN', '.STM' /
C=======================================================================
FNAME(1)=PROJECT//EXTEN(1)
FNAME(2)=PROJECT//EXTEN(2)
FNAME(3)=PROJECT//EXTEN(3)
OPEN ( UNIT=2, FILE='NSDATA4.FLN', STATUS='UNKNOWN')
WRITE (2,'(A15)') PROJECT
WRITE (2,'(A19)') FNAME(1)
WRITE (2,'(A19)') FNAME(2)
WRITE (2,'(A19)') FNAME(3)
CLOSE (2)
WRITE (*,*) PROJECT
WRITE (*,*) FNAME(1)
WRITE (*,*) FNAME(2)
WRITE (*,*) FNAME(3)
C=======================================================================
C C1: ADDED VISCOSITY COMPUTATION PARAMETER
C C1=0. FOR NO ADDED VISCOSITY, C1=1. FOR UPWIND (OPTIMUM) ADDED
C C1 OPTION IS AVALABLE IN NSEQ8DD.FOR SOLVER
C NEY: NEMBER OF VERTICAL ELEMENTS (NUMBER OF NODES: NEY+1)
C NEX: NEMBER OF HORIZONTAL ELEMENTS (NUMBER OF NODES: NEX+1)
C MAXITE: MAXIMUM NUMBER OF ITERATIONS AT EACH TIME INCREMENT
C ITERATION PROCEDURE IS REQUIRED BECAUSE OF NON-LINEAR CONVECTIVE TERM.
C ERMAX: MAXIMUM DIFFERENCE BETWEEN GUESSED {U&V} AND COMPUTED {U&V}
C HEIGHT: HEIGHT OF DRIVEN CAVITY FLOW DOMAIN
C=======================================================================
WRITE(*,*)' PROJECT NAME = ', PROJECT
WRITE(*,*)' VISCOSITY RATIO (2ND/1ST) =', VRATIO
C=======================================================================
C PREPARATION FOR ELEMENT CREATION
C-------ITYPE=1: 1 3 5 7 ............... 8 6 4 2 ( ZIGZAG: PART 1 )
C-------ITYPE=2: ....... 7 5 3 1 2 4 6 ......... ( ZIGZAG: PART 2 )
C-------ITYPE=3: NDX NDX-1 ........... 5 4 3 2 1 ( DESCENDING ORDER )
C-------ITYPE=4: 1 2 3 4 5 6 7 8 ........... NDX ( INCRESING ORDER )
WRITE(*,*)' ITYPE =',ITYPE
CALL NUMBERIN ( NDY , ITYPE , NODES )
WRITE(*,*)' NODES IN Y COORDINATE'
WRITE(*,*) NODES
C=======================================================================
C ELEMENT CREATION
NE = 0
DO I = 1 , NEY
DO J = 1 , NEX
NE = NE + 1
IF ( NE .GT. MXE ) STOP 'NE > MXE'
NODEX(NE,1) = NDY*(J-1) + NODES(I)
NODEX(NE,2) = NDY* J + NODES(I)
NODEX(NE,3) = NDY* J + NODES(I+1)
NODEX(NE,4) = NDY*(J-1) + NODES(I+1)
END DO
END DO
C=======================================================================
C NODAL COORDINATE CREATION
NNODE = 0
DO I = 1 , NDY
DO J = 1 , NDX
NNODE = NNODE + 1
IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
NODE = NDY*(J-1) + NODES(I)
XCOORD(NODE) = TLX * F( DX*(J-1)/TLX )
YCOORD(NODE) = TLY * G( DY*(I-1)/TLY )
END DO
END DO
C=======================================================================
C BOUNDARY CONDITIONS
C--------- MOMENTUM EQUATIONS AND HEAT EQUATION
NBFX = 0
NBFY = 0
NBT = 0
C--------- FACE OF +Y AND -Y
U0 = 0.
V0 = 0.
DO J = 2 , NEX
C--------- -Y
NBFX = NBFX + 1
NBFY = NBFY + 1
IBNDFX(NBFX) = NODES(1) + (J-1)*NDY
IBNDFY(NBFY) = NODES(1) + (J-1)*NDY
BVX(NBFX) = U0
BVY(NBFY) = V0
C--------- +Y
NBFX = NBFX + 1
NBFY = NBFY + 1
IBNDFX(NBFX) = NDY*(J-1) + NODES(NDY)
IBNDFY(NBFY) = NDY*(J-1) + NODES(NDY)
BVX(NBFX) = U0
BVY(NBFY) = V0
END DO
C--------- FACE OF +X AND -X
C -X
DO I = 1 , NDY
NBFX = NBFX + 1
NBFY = NBFY + 1
NBT = NBT + 1
IBNDFX(NBFX) = NODES(I)
IBNDFY(NBFY) = NODES(I)
IBNDT (NBT ) = NODES(I)
BVX(NBFX) = U0
BVY(NBFY) = V0
BVT(NBT ) = T1
C +X
NBFX = NBFX + 1
NBFY = NBFY + 1
NBT = NBT + 1
IBNDFX(NBFX) = NDY*NEX + NODES(I)
IBNDFY(NBFY) = NDY*NEX + NODES(I)
IBNDT (NBT ) = NDY*NEX + NODES(I)
BVX(NBFX) = U0
BVY(NBFY) = V0
BVT(NBT ) = T2
END DO
C=======================================================================
C INITIAL CONDITIONS
U0 = 0.
V0 = 0.
DO I = 1 , NNODE
U(I) = U0
V(I) = V0
TEMP(I) = TREF
END DO
C=======================================================================
C B. C. FOR STREAM FUNCTION
NBS = NBFX
DO I = 1 , NBS
IBNDS(I) = IBNDFX(I)
BVS(I) = BVX(I)
END DO
C=======================================================================
C DATA FILE INQUIRY
EXFILE = 'NEW'
INQUIRE ( FILE = FNAME(1), EXIST = YES )
IF ( YES ) EXFILE='OLD'
C=======================================================================
C MAKING DATA FILES
C---------- 'PROJECT'.DAT
OPEN ( 1, FILE=FNAME(1), STATUS = EXFILE )
WRITE(1,*) VISCO, FLMDA, TREF, CONDUCT, C1
WRITE(1,*) TMAX, DTMAX
WRITE(1,*) MAXITE, ERMAX, ITEFIX
C---------- ELEMENT
WRITE(1,*) NE
DO I = 1 , NE
WRITE (1,*) I, (NODEX(I,J), J = 1 , ND )
END DO
C---------- NODE
WRITE(1,*) NNODE
DO I = 1 , NNODE
WRITE(1,*) I, XCOORD(I), YCOORD(I)
END DO
C---------- BOUNDARY CONDITIONS
WRITE(1,*) NBFX
DO I = 1 , NBFX
WRITE (1,*) IBNDFX(I), BVX(I)
END DO
WRITE(1,*) NBFY
DO I = 1 , NBFY
WRITE (1,*) IBNDFY(I), BVY(I)
END DO
WRITE(1,*) NBT
DO I = 1 , NBT
WRITE (1,*) IBNDT(I), BVT(I)
END DO
CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.BIN
OPEN ( 1, FILE=FNAME(2), STATUS = 'UNKNOWN' ,FORM='UNFORMATTED' )
WRITE(1) TIME, DT
WRITE(1) ( U(I) , I = 1 , NNODE )
WRITE(1) ( V(I) , I = 1 , NNODE )
WRITE(1) ( TEMP(I) , I = 1 , NNODE )
CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.STM
OPEN ( 1, FILE=FNAME(3), STATUS = 'UNKNOWN' )
WRITE(1,*) NBS
WRITE (1,*) ( IBNDS(I) , BVS(I) , I = 1 , NBS )
CLOSE (1)
C=======================================================================
CALL BANDWID ( ND, MXE, NE, NODEX , NBW )
C------ CREATION OF PARAMETER FILE TO BE USED IN INCLUDE STATEMENT
OPEN ( 1, FILE='PARAM.DAT', STATUS='UNKNOWN' )
WRITE (1,*) ' PARAMETER ( NF=',NF,' )'
WRITE (1,*) ' PARAMETER ( ND=',ND,' )'
WRITE (1,*) ' PARAMETER ( INTEPT=',INTEPT,')'
WRITE (1,*) ' PARAMETER ( MXE=',NE,')'
WRITE (1,*) ' PARAMETER ( MXN=',NNODE,')'
WRITE (1,*) ' PARAMETER ( MXW=',NBW,')'
NB = MAX0 (NBFX, NBFY, NBS, NBT )
WRITE (1,*) ' PARAMETER ( MXB=',NB,' )'
CLOSE (1)
C=======================================================================
STOP "NORMAL TERMINATION"
END
C
C
SUBROUTINE BANDWID ( ND, MXE, NE, NODEX , NBW )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND)
NBW = 0
DO I = 1 , NE
DO J = 1 , ND-1
DO K = J+1 , ND
NBW = MAX0(NBW,IABS(NODEX(I,J)-NODEX(I,K)))
END DO
END DO
END DO
NBW = NBW + 1
WRITE(*,*) 'HALF BANDWIDTH + 1 =', NBW
RETURN
END
C
C
FUNCTION F(XX)
IMPLICIT REAL*8 ( A-H , O-Z )
DATA SLOPE / 0.3D0 /
Z(X,S) = ( X * X * (3.-S+(S-2.) * X ) )
X = XX * 2.
IF ( X .LE. 1. ) THEN
F = 0.5 * Z(X,SLOPE)
ELSE
X = 2. - X
F = 0.5 * ( 1.-Z(X,SLOPE)) + 0.5
ENDIF
RETURN
END
C
C
FUNCTION G(X)
IMPLICIT REAL*8 ( A-H , O-Z )
DATA X1 / 0.2D0 /
H(X) = ( X*X*(3.-2.*X) )
IF ( X .LE. X1 ) G = H(X)
IF ( X .GE. (1.-X1) ) G = H(X)
IF ( X.GT.X1 .AND. X.LT.(1.-X1) ) THEN
SLOPE = ( 1.- 2.*H(X1))/(1.-2.*X1)
B = H(X1) - SLOPE*X1
G = SLOPE * X + B
ENDIF
RETURN
END
C
C
SUBROUTINE NUMBERIN ( NDX , ITYPE , NODES )
DIMENSION NODES(NDX)
NEX = NDX - 1
NEGATIVE = -1
J = NEGATIVE
IF ( ITYPE .EQ. 1 ) THEN
ICND = 1
NODES(ICND) = 1
DO I = NEX , 1 , -1
J = J*NEGATIVE
ICND = ICND + J*I
NODES(ICND) = NDX - I + 1
END DO
ENDIF
C
IF ( ITYPE .EQ. 2 ) THEN
ICND = NDX/2 + ( NDX - NDX/2*2 )
NODES(ICND) = 1
DO I = 1 , NEX
J = J*NEGATIVE
ICND = ICND + J*I
NODES(ICND) = I + 1
END DO
ENDIF
C
IF ( ITYPE .EQ. 3 ) THEN
DO I = 1 , NDX
NODES(I) = NDX - I + 1
END DO
ENDIF
C
IF ( ITYPE .EQ. 4 ) THEN
DO I = 1 , NDX
NODES(I) = I
END DO
ENDIF
RETURN
END