PROGRAM SET9NSEQ12
C=======================================================================
C DATA GENERATING PROGRAM FOR NSEQ8XX.FOR
C PROJECT: DRIVEN CAVITY FLOW PROJECT NAME: SEE BELOW DATA PROJECT
C DOMAIN: SQUARE ( 1 X 1 )
C ELEMENT: 9-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 JAN. 29, 2013
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=9, MXE=12000, MXN=13000, MXB=1610, NF=3 )
PARAMETER ( REYNLD=1000.D0, VRATIO=1.D7, HEIGHT=1.D0, UAVE=1.D0,
* VISCO=HEIGHT*UAVE/REYNLD, FLMDA=VISCO*VRATIO,
* DT = 0.025D0, TMAX=200.D0, MAXITE=10, ERMAX=0.001D0, C1=0.D0,
* TLX = HEIGHT, TLY = HEIGHT, NEX = 20, NEY = 20,
* TIME = 0.D0, DTMAX=0.05D0, ITEFIX=3 )
PARAMETER ( INTEPT=3 )
C=======================================================================
DIMENSION NODEX(MXE,ND), XCOORD(MXN), YCOORD(MXN), U(MXN), V(MXN),
* IBNDFX(MXB),IBNDFY(MXB),BVX(MXB),BVY(MXB),IBNDS(MXB),BVS(MXB)
CHARACTER FNAME(NF)*19,EXTEN(NF)*4,PROJECT*15,EXFILE*3
LOGICAL YES
C=======================================================================
DATA PROJECT / 'DCF9REYNLD-1000' /
DATA EXTEN / '.DAT','.BIN','.STM' /
C=======================================================================
F(X) = X * X * ( 3.- 2.* X )
G(X) = X * ( (0.1-1.)* X + (2.- 0.1) )
C=======================================================================
FNAME(1)=PROJECT//EXTEN(1)
FNAME(2)=PROJECT//EXTEN(2)
FNAME(3)=PROJECT//EXTEN(3)
OPEN ( UNIT=2, FILE='NSDATA9.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 (*,*)' DRIVEN CAVITY FLOW: REYNOLDS NUMBER OF ',REYNLD
WRITE (*,*)' VISCOSITY RATIO OF ',VRATIO
C=======================================================================
NDX = NEX*2 + 1
NDY = NEY*2 + 1
DX = TLX / (NEX*2)
DY = TLY / (NEY*2)
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) = (J-1)*2 + 1 + (I-1)*(2*NDX)
NODEX(NE,2) = NODEX(NE,1) + 1
NODEX(NE,3) = NODEX(NE,2) + 1
NODEX(NE,4) = NODEX(NE,3) + NDX
NODEX(NE,5) = NODEX(NE,4) + NDX
NODEX(NE,6) = NODEX(NE,5) - 1
NODEX(NE,7) = NODEX(NE,6) - 1
NODEX(NE,8) = NODEX(NE,1) + NDX
NODEX(NE,9) = NODEX(NE,8) + 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'
XCOORD(NNODE) = TLX * F( DX*(J-1)/TLX )
YCOORD(NNODE) = TLY * G( DY*(I-1)/TLY )
END DO
END DO
C=======================================================================
C BOUNDARY CONDITIONS
C--------- MOMENTUM EQUATIONS
C--------- FACE OF -X
NBFX = 0
NBFY = 0
U0 = 0.
V0 = 0.
DO I = 2 , NDY-2
NBFX = NBFX + 1
NBFY = NBFY + 1
IBNDFX(NBFX) = NDX*(I-1)+1
IBNDFY(NBFY) = NDX*(I-1)+1
BVX(NBFX) = U0
BVY(NBFY) = V0
END DO
C--------- FACE OF -Y
U0 = 0.
V0 = 0.
DO I = 1 , NDX
NBFX = NBFX + 1
IBNDFX(NBFX) = I
NBFY = NBFY + 1
IBNDFY(NBFY) = I
BVX(NBFX) = U0
BVY(NBFY) = V0
END DO
C--------- FACE OF +Y
U0 = UAVE
V0 = 0.
DO I = 1 , NDX
NBFX = NBFX + 1
IBNDFX(NBFX) = NDX*(NDY-1) + I
NBFY = NBFY + 1
IBNDFY(NBFY) = NDX*(NDY-1) + I
BVX(NBFX) = U0
BVY(NBFY) = V0
END DO
C--------- FACE OF +X
U0 = 0.
V0 = 0.
DO I = 2 , NDY-2
NBFX = NBFX + 1
IBNDFX(NBFX) = NDX*I
NBFY = NBFY + 1
IBNDFY(NBFY) = NDX*I
BVX(NBFX) = U0
BVY(NBFY) = V0
END DO
C=======================================================================
C INITIAL CONDITIONS
U0 = 0.
V0 = 0.
DO I = 1 , NNODE
U(I) = U0
V(I) = V0
END DO
C=======================================================================
C B. C. FOR STREAM FUNCTION
DDY = YCOORD(NNODE) - YCOORD(NNODE-2*NDX)
Q = DDY * UAVE / 2.
NBS = NBFX
DO I = 1 , NBS
IBNDS(I) = IBNDFX(I)
BVS(I) = BVX(I)
IF ( BVX(I) .EQ. UAVE ) BVS(I) = Q
END DO
C=======================================================================
C MAKING DATA FILES
OPEN ( 1, FILE=FNAME(1), STATUS = 'UNKNOWN')
C---------- HEADER
WRITE(1,*) VISCO, FLMDA
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---------- NORDAL INFORMATION
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
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 )
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=======================================================================
C---------- ELEMENT DRAWING
OPEN ( 1, FILE='ELEMENT9.DAT', STATUS = 'UNKNOWN')
DO I = 1, NE
DO J = 1, 3
WRITE (1,*) XCOORD(NODEX(I,J)), YCOORD(NODEX(I,J))
END DO
WRITE (1,*)
DO J = 3, 5
WRITE (1,*) XCOORD(NODEX(I,J)), YCOORD(NODEX(I,J))
END DO
WRITE (1,*)
DO J = 5, 7
WRITE (1,*) XCOORD(NODEX(I,J)), YCOORD(NODEX(I,J))
END DO
WRITE (1,*)
DO J = 7, 8
WRITE (1,*) XCOORD(NODEX(I,J)), YCOORD(NODEX(I,J))
END DO
WRITE (1,*) XCOORD(NODEX(I,1)), YCOORD(NODEX(I,1))
WRITE (1,*)
WRITE (1,*) XCOORD(NODEX(I,9)), YCOORD(NODEX(I,9))
WRITE (1,*)
END DO
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 )
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