PROGRAM SET3D27NSEQ1_NEWFORMAT
C=======================================================================
C DATA GENERATING PROGRAM FOR 3D-NSEQ8XX.FOR
C PROJECT: DRIVEN CAVITY FLOW PROJECT NAME: DCFR999
C DOMAIN: SQUARE ( 1 X 1 X 1 )
C ELEMENT: 3D HEXA 27-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 JULY 11, 2010
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=27, MXE=22000, MXN=53000, MXB=20000, NF=9 )
C=======================================================================
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),ZCOORD(MXN),
* U(MXN),V(MXN),W(MXN),
* IBNDFX(MXB),IBNDFY(MXB),IBNDFZ(MXB),BVX(MXB),BVY(MXB),BVZ(MXB)
CHARACTER FNAME(NF)*11,EXTEN(NF)*4,PROJECT*7,EXFILE*3
LOGICAL YES
C=======================================================================
DATA PROJECT / 'DCFR999' /
DATA EXTEN / '.JNK', '.ELE', '.NOD', '.BND', '.INI', '.BIN',
* '.XXX', '.MAX', '.STM' /
C=======================================================================
F(X) = X * X * ( 3.D0- 2.D0* X )
G(X) = X * ( (0.1D0-1.D0)* X + (2.D0- 0.1D0) )
C=======================================================================
C----------> FLOW DOMAIN PARARETERS
TLX = 1.D0
TLY = 1.D0
TLZ = 1.D0
C----------> ELEMENT INFORMATION
NEX = 10
NEY = 10
NEZ = 10
C----------> ELEMENT INFORMATION CONTINUED
DX = TLX/(2*NEX)
DY = TLY/(2*NEY)
DZ = TLZ/(2*NEZ)
NDX = 2*NEX + 1
NDY = 2*NEY + 1
NDZ = 2*NEZ + 1
C=======================================================================
C----------> FLUID PARAMETERS
REYNLD = 1000.D0
VRATIO = 1.D7
UAVE = 1.D0
VISCO = TLZ*UAVE/REYNLD
FLMDA = VISCO*VRATIO
C=======================================================================
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 TLZ: HEIGHT OF DRIVEN CAVITY FLOW DOMAIN
C=======================================================================
CALL FILEMG ( NF, FNAME, EXTEN, PROJECT )
WRITE (*,*)'DRIVEN CAVITY FLOW'
WRITE (*,*)'REYNOLDS NUMBER OF ',REYNLD
WRITE (*,*)'VISCOSITY RATIO OF ',VRATIO
WRITE (*,*)'VISCOSITY ',VISCO
WRITE (*,*)'2ND VISCOSITY ',FLMDA
C=======================================================================
DT = 0.005D0
TMAX = 200.D0
MAXITE = 10
TIME = 0.D0
ERMAX = 0.001D0
C=======================================================================
C NODAL COORDINATE CREATION
NNODE = 0
DO I = 1 , NDX
DO K = 1 , NDZ
DO J = 1 , NDY
NNODE = NNODE + 1
IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
XCOORD(NNODE) = TLX * F( DX*(I-1)/TLX )
YCOORD(NNODE) = TLY * F( DY*(J-1)/TLY )
ZCOORD(NNODE) = TLZ * G( DZ*(K-1)/TLZ )
IF ( K .EQ. NDZ ) ZCOORD(NNODE) = TLZ
END DO
END DO
END DO
WRITE (*,*) 'NNODE =', NNODE
C=======================================================================
C ELEMENT CREATION
NE = 0
DO I = 1 , NEX
DO J = 1 , NEZ
DO K = 1 , NEY
NE = NE + 1
IF ( NE .GT. MXE ) STOP 'NE > MXE'
C------ BOTTOM
NODEX(NE,1) = (I-1)*NDY*NDZ*2 + (J-1)*NDY*2 + (K-1)*2 + 1
NODEX(NE,2) = NODEX(NE,1) + NDY*NDZ
NODEX(NE,3) = NODEX(NE,2) + NDY*NDZ
NODEX(NE,4) = NODEX(NE,3) + 1
NODEX(NE,5) = NODEX(NE,4) + 1
NODEX(NE,8) = NODEX(NE,1) + 1
NODEX(NE,7) = NODEX(NE,8) + 1
NODEX(NE,9) = NODEX(NE,2) + 1
NODEX(NE,6) = NODEX(NE,9) + 1
C------ MID
NODEX(NE,10) = NODEX(NE,1) + NDY
NODEX(NE,17) = NODEX(NE,10) + 1
NODEX(NE,16) = NODEX(NE,17) + 1
NODEX(NE,11) = NODEX(NE,10) + NDY*NDZ
NODEX(NE,18) = NODEX(NE,11) + 1
NODEX(NE,15) = NODEX(NE,18) + 1
NODEX(NE,12) = NODEX(NE,11) + NDY*NDZ
NODEX(NE,13) = NODEX(NE,12) + 1
NODEX(NE,14) = NODEX(NE,13) + 1
C------ TOP
NODEX(NE,19) = NODEX(NE,10) + NDY
NODEX(NE,26) = NODEX(NE,19) + 1
NODEX(NE,25) = NODEX(NE,26) + 1
NODEX(NE,20) = NODEX(NE,19) + NDY*NDZ
NODEX(NE,27) = NODEX(NE,20) + 1
NODEX(NE,24) = NODEX(NE,27) + 1
NODEX(NE,21) = NODEX(NE,20) + NDY*NDZ
NODEX(NE,22) = NODEX(NE,21) + 1
NODEX(NE,23) = NODEX(NE,22) + 1
END DO
END DO
END DO
WRITE (*,*) 'NE =', NE
C=======================================================================
CALL BANDWID ( MXE, ND, NE, NODEX, NBW )
C=======================================================================
C BOUNDARY CONDITIONS
C--------- MOMENTUM EQUATIONS
NBFX = 0
NBFY = 0
NBFZ = 0
C--------- FACE OF -X
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE -X
DO K = 2 , NDZ-1
DO J = 1 , NDY
C
NBFX = NBFX + 1
IBNDFX(NBFX) = J + (K-1)*NDY
BVX(NBFX) = 0.D0
C
NBFY = NBFY + 1
IBNDFY(NBFY) = J + (K-1)*NDY
BVY(NBFY) = 0.D0
C
NBFZ = NBFZ + 1
IBNDFZ(NBFZ) = J + (K-1)*NDY
BVZ(NBFZ) = 0.D0
C
END DO
END DO
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE +X
DO K = 2 , NDZ-1
DO J = 1 , NDY
C
NBFX = NBFX + 1
IBNDFX(NBFX) = J + (K-1)*NDY + (NDX-1)*NDY*NDZ
BVX(NBFX) = 0.D0
C
NBFY = NBFY + 1
IBNDFY(NBFY) = J + (K-1)*NDY + (NDX-1)*NDY*NDZ
BVY(NBFY) = 0.D0
C
NBFZ = NBFZ + 1
IBNDFZ(NBFZ) = J + (K-1)*NDY + (NDX-1)*NDY*NDZ
BVZ(NBFZ) = 0.D0
C
END DO
END DO
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE -Y
DO K = 2 , NDZ-1
DO I = 2 , NDX-1
C
NBFX = NBFX + 1
IBNDFX(NBFX) = 1 + (K-1)*NDY + (I-1)*NDY*NDZ
BVX(NBFX) = 0.D0
C
NBFY = NBFY + 1
IBNDFY(NBFY) = 1 + (K-1)*NDY + (I-1)*NDY*NDZ
BVY(NBFY) = 0.D0
C
NBFZ = NBFZ + 1
IBNDFZ(NBFZ) = 1 + (K-1)*NDY + (I-1)*NDY*NDZ
BVZ(NBFZ) = 0.D0
C
END DO
END DO
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE +Y
DO K = 2 , NDZ-1
DO I = 2 , NDX-1
C
NBFX = NBFX + 1
IBNDFX(NBFX) = K*NDY + (I-1)*NDY*NDZ
BVX(NBFX) = 0.D0
C
NBFY = NBFY + 1
IBNDFY(NBFY) = K*NDY + (I-1)*NDY*NDZ
BVY(NBFY) = 0.D0
C
NBFZ = NBFZ + 1
IBNDFZ(NBFZ) = K*NDY + (I-1)*NDY*NDZ
BVZ(NBFZ) = 0.D0
C
END DO
END DO
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE -Z
DO I = 1 , NDX
DO J = 1 , NDY
C
NBFX = NBFX + 1
IBNDFX(NBFX) = J + (I-1)*NDY*NDZ
BVX(NBFX) = 0.D0
C
NBFY = NBFY + 1
IBNDFY(NBFY) = J + (I-1)*NDY*NDZ
BVY(NBFY) = 0.D0
C
NBFZ = NBFZ + 1
IBNDFZ(NBFZ) = J + (I-1)*NDY*NDZ
BVZ(NBFZ) = 0.D0
C
END DO
END DO
C----------------------------------------------
C---------- ALL DIRECTIONS ON FACE +Z
DO I = 1 , NDX
DO J = 1 , NDY
C
NBFX = NBFX + 1
IBNDFX(NBFX) = I*NDY*NDZ - NDY + J
BVX(NBFX) = UAVE
C
NBFY = NBFY + 1
IBNDFY(NBFY) = I*NDY*NDZ - NDY + J
BVY(NBFY) = 0.D0
C
NBFZ = NBFZ + 1
IBNDFZ(NBFZ) = I*NDY*NDZ - NDY + J
BVZ(NBFZ) = 0.D0
C
END DO
END DO
WRITE (*,*) 'NBFX =', NBFX
WRITE (*,*) 'NBFY =', NBFY
WRITE (*,*) 'NBFZ =', NBFZ
C=======================================================================
C INITIAL CONDITIONS
DO I = 1 , NNODE
U(I) = 0.D0
V(I) = 0.D0
W(I) = 0.D0
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'.JNK
OPEN ( 1, FILE=FNAME(1), STATUS = EXFILE )
WRITE(1,*) VISCO, FLMDA
WRITE(1,*) DT, TMAX
WRITE(1,*) MAXITE, ERMAX
CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.ELE
OPEN ( 1, FILE=FNAME(2), STATUS = EXFILE )
WRITE(1,*) NE
DO I = 1 , NE
WRITE (1,*) I, (NODEX(I,J), J = 1 , ND )
END DO
CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.NOD
OPEN ( 1, FILE=FNAME(3), STATUS = EXFILE )
WRITE(1,*) NNODE
DO I = 1, NNODE
WRITE(1,*) I,XCOORD(I), YCOORD(I), ZCOORD(I)
END DO
CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.BND
OPEN ( 1, FILE=FNAME(4), STATUS = EXFILE )
WRITE(1,*) NBFX
DO I = 1 , NBFX
WRITE (1,*) IBNDFX(I), BVX(I)
END DO
C
WRITE(1,*) NBFY
DO I = 1 , NBFY
WRITE (1,*) IBNDFY(I), BVY(I)
END DO
C
WRITE(1,*) NBFZ
DO I = 1 , NBFZ
WRITE (1,*) IBNDFZ(I), BVZ(I)
END DO
CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.BIN
OPEN ( 1, FILE=FNAME(6), STATUS = EXFILE,FORM='UNFORMATTED' )
WRITE(1) TIME
WRITE(1) ( U(I) , I = 1 , NNODE )
WRITE(1) ( V(I) , I = 1 , NNODE )
WRITE(1) ( W(I) , I = 1 , NNODE )
CLOSE (1)
C=======================================================================
C---------- 'PROJECT'.MAX
OPEN ( 1, FILE=FNAME(8), STATUS = EXFILE )
UMIN = 0.
UMAX = 0.
VMIN = 0.
VMAX = 0.
PMIN = 0.
PMAX = 0.
SMIN = 0.
SMAX = 0.
WRITE (1,*) UMIN , UMAX
WRITE (1,*) VMIN , VMAX
WRITE (1,*) PMIN , PMAX
WRITE (1,*) SMIN , SMAX
CLOSE (1)
C=======================================================================
STOP
END
C
C
SUBROUTINE FILEMG ( NF, FNAME, EXTEN, PROJECT )
CHARACTER FNAME(NF)*11, PROJECT*7, EXTEN(NF)*4, FNWOE(7)*1,
* EXFILE*3
LOGICAL YES
C------- RETURN VALUE: FNAME
C
EXFILE ='NEW'
INQUIRE ( FILE='NSDATA.FLN', EXIST=YES )
IF ( YES ) EXFILE='OLD'
OPEN ( UNIT=2, FILE='NSDATA.FLN', STATUS=EXFILE )
C
NCHAR = 7
DO I = 1 , NCHAR
FNWOE(I) = PROJECT (I:I)
END DO
CALL GENFNM ( NCHAR, FNWOE, NF, FNAME, EXTEN )
WRITE (2,'(I5)') NCHAR
WRITE (2,'(7A1)') ( FNWOE(I) , I = 1 , NCHAR )
DO I = 1 , NF
WRITE (*,'(1X,A11)') FNAME(I)
WRITE (2,'( A11)') FNAME(I)
END DO
CLOSE (2)
RETURN
END
C
C
SUBROUTINE GENFNM ( NCHAR, FNWOE, NF, FNAME, EXTEN )
CHARACTER*1 FNWOE(7)
CHARACTER*4 EXTEN(NF)
CHARACTER*11 FNAME(NF)
DO I = 1 , NF
DO J = 1 , NCHAR
FNAME (I) (J:J) = FNWOE(J)
END DO
DO J = NCHAR+1,NCHAR+4
FNAME (I) (J:J) = EXTEN(I) (J-NCHAR:J-NCHAR)
END DO
END DO
RETURN
END
C
C
SUBROUTINE BANDWID ( MXE , ND , NE , NODEX , NBW )
DIMENSION NODEX(MXE,ND)
C------- RETURN VALUE: NBW
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 (*,*) ' HALH BANDWIDTH =', NBW
RETURN
END