PROGRAM SET4NSEQ14
C========= FLOW AROUND CYLINDER ========
C=======================================================================
C DATA GENERATING PROGRAM FOR NSEQ8XX.FOR
C PROJECT: DRIVEN CAVITY FLOW PROJECT NAME: CYLINDR
C DOMAIN: RECTANGULER WITH A CYLINDER OBJECT IN FLOW FIELD
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 JAN. 31, 2013
C------------- DIAMETER OF CYLINDER = 1. -------------------------------
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=4, MXE=2200, MXN=4400, MXB=616, NF=3 )
PARAMETER ( VRATIO=1.E7, UAVE=1.,REYNLD=250.D0,
* VISCO = UAVE / REYNLD, FLMDA=VISCO*VRATIO,
* DT = 0.05,TMAX=100., MAXITE=10, ERMAX=0.001, C1=0., TIME=0.,
* ISIDE = 2, IFRONT = 10, IBACK = 20, NEX = IFRONT + 8 + IBACK,
* NEY = ISIDE + 8 + ISIDE, DX = 0.4, DY = 0.4,DTMAX=0.05,
* NDX = NEX + 1, NDY = NEY + 1 , ITEFIX=4 )
PARAMETER ( INTEPT=2 )
C======================================================================
C
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
* IBNDFX(MXB),IBNDFY(MXB),BVX(MXB),BVY(MXB),U(MXN),V(MXN),
* IBNDS(MXB),BVS(MXB)
DIMENSION ICR(7,32),ICREL(192),ICRND(192),XC(32),YC(32)
DIMENSION IREL(MXN,10),INOD(MXN),NEW(MXN),JNT(MXN),
* NJNT(MXN),XTEMP(MXN),YTEMP(MXN)
CHARACTER FNAME(NF)*19,EXTEN(NF)*4,PROJECT*15,EXFILE*3
LOGICAL YES
C=======================================================================
DATA PROJECT / 'FLOWARNDCYLINDR' /
DATA EXTEN / '.DAT', '.BIN', '.STM' /
C------------ COORDINATES OF BASE CYLINDER (-1, -1) - (+1, +1) ---------
DATA XC /
*-.707,-.555,-.383,-.195,0.,0.195,0.383,0.555,0.707,0.831,
*0.924,0.981,1.,0.981,0.924,0.831,0.707,0.555,0.383,0.195,
*0.,-.195,-.383,-.555,-.707,-.831,-.924,-.981,-1.00,-.981,
*-.924,-.831 /
DATA YC /
*-.707,-.831,-.924,-.981,-1.,-.981,-.924,-.831,-.707,-.555,
*-.383,-.195,0.,0.195,0.383,0.555,0.707,0.831,0.924,0.981,
*1.,0.981,0.924,0.831,0.707,0.555,0.383,0.195,0.,-0.195,
*-0.383,-0.555 /
C=======================================================================
F(X) = X
G(X) = X
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=======================================================================
WRITE (*,*)' FLOW AROUND CYLINDER'
WRITE (*,*)' VISCOSITY RATIO OF ',VRATIO
C=======================================================================
NE = 0
DO 10 J = 1 , NEX
DO 10 I = 1 , NEY
NE = NE + 1
NODEX(NE,1) = NDY*(J-1) + I
NODEX(NE,2) = NODEX(NE,1) + NDY
NODEX(NE,3) = NODEX(NE,2) + 1
NODEX(NE,4) = NODEX(NE,1) + 1
10 CONTINUE
C
TLX = NEX * DX
TLY = NEY * DY
NNODE = 0
DO 20 I = 1 , NDX
DO 20 J = 1 , NDY
NNODE = NNODE + 1
XCOORD(NNODE) = TLX*F( DX*(I-1)/TLX )
YCOORD(NNODE) = TLY*G( DY*(J-1)/TLY )
20 CONTINUE
C======== CENTER OF CYLINDER ========
I = NDY*(IFRONT+4) + ISIDE + 4 + 1
XCENTER = XCOORD (I)
YCENTER = YCOORD (I)
C======== NODES AND ELEMENTS AROUND CYLINDER =======
K = 0
DO 71 I = 1 , 8
K = K + 1
ICR(7,K) = NDY*(IFRONT+I-1) + ISIDE + 1
71 CONTINUE
DO 72 I = 1 , 8
K = K + 1
ICR(7,K) = NDY*(IFRONT+8) + ISIDE + I
72 CONTINUE
DO 73 J = 1 , 8
I = 8 - J + 1
K = K + 1
ICR(7,K) = NDY*(IFRONT+I) + ISIDE + 8 + 1
73 CONTINUE
DO 74 J = 1 , 8
I = 8 - J + 1
K = K + 1
ICR(7,K) = NDY*IFRONT + ISIDE + 1 + I
74 CONTINUE
K = 0
DO 75 I = 1 , 7
DO 75 J = 1 , 7
K = K + 1
ICRND(K) = NDY*(IFRONT+1+I-1) + ISIDE + 1 + J
75 CONTINUE
DO 76 I = K+1 , 192
NNODE = NNODE + 1
ICRND(I) = NNODE
76 CONTINUE
K = 0
DO 77 I = 1 , 6
DO 77 J = 1 , 32
K = K + 1
ICR(I,J) = ICRND(K)
77 CONTINUE
K = 0
DO 78 I = 1 , 8
DO 78 J = 1 , 8
K = K + 1
ICREL(K) = NEY*(IFRONT+I-1) + ISIDE + J
78 CONTINUE
DO 79 I = K+1 , 192
NE = NE + 1
ICREL(I) = NE
79 CONTINUE
WRITE (*,*)' NUMBER OF ELEMENTS =',NE
WRITE (*,*)' NUMBER OF NODES =',NNODE
K = 0
DO 55 I = 1 , 6
DO 44 J = 1 , (32-1)
K = K + 1
IEL = ICREL(K)
NODEX(IEL,1) = ICR(I+1,J)
NODEX(IEL,2) = ICR(I+1,J+1)
NODEX(IEL,3) = ICR(I,J+1)
NODEX(IEL,4) = ICR(I,J)
44 CONTINUE
K = K + 1
IEL = ICREL(K)
NODEX(IEL,1) = ICR(I+1,32)
NODEX(IEL,2) = ICR(I+1,1)
NODEX(IEL,3) = ICR(I,1)
NODEX(IEL,4) = ICR(I,32)
55 CONTINUE
C------- COORDINATES -----
DO 41 I = 1 , 32
XCOORD(ICR(1,I)) = 1.25*DX*XC(I) + XCENTER
YCOORD(ICR(1,I)) = 1.25*DX*YC(I) + YCENTER
41 CONTINUE
DO 51 J = 1 , 32
RX = (XCOORD(ICR(7,J))-XCOORD(ICR(1,J))) / 6.
RY = (YCOORD(ICR(7,J))-YCOORD(ICR(1,J))) / 6.
DO 51 I = 2 , 6
XCOORD(ICR(I,J)) = XCOORD(ICR(1,J)) + RX*(I-1)
YCOORD(ICR(I,J)) = YCOORD(ICR(1,J)) + RY*(I-1)
51 CONTINUE
C
C--------- BOUNDARY CONDITIONS ---------
C............ MOMENTUM EQUATIONS .......
C FACE OF -X
NBFX = 0
NBFY = 0
U0 = UAVE
V0 = 0.
DO 30 J = 2 , ( NDY-1)
NBFX = NBFX + 1
NBFY = NBFY + 1
IBNDFX(NBFX) = J
IBNDFY(NBFY) = J
BVX(NBFX) = U0
BVY(NBFY) = V0
30 CONTINUE
C FACE OF -Y
U0 = 0.
V0 = 0.
DO 40 I = 1 , NDX
NBFX = NBFX + 1
NBFY = NBFY + 1
IBNDFX(NBFX) = 1 + NDY*(I-1)
IBNDFY(NBFY) = 1 + NDY*(I-1)
BVX(NBFX) = U0
BVY(NBFY) = V0
40 CONTINUE
C FACE OF +Y
U0 = 0.
V0 = 0.
DO 50 I = 1 , NDX
NBFX = NBFX + 1
NBFY = NBFY + 1
IBNDFX(NBFX) = I*NDY
IBNDFY(NBFY) = I*NDY
BVX(NBFX) = U0
BVY(NBFY) = V0
50 CONTINUE
C------- FACE OF +X
CCCCCCCCCCC U0 = UAVE
V0 = 0.
DO 29 I = 2 , NEY
CCCCCCCCCCC NBFX = NBFX + 1
NBFY = NBFY + 1
CCCCCCCCCCC IBNDFX(NBFX) = NEX*NDY + I
IBNDFY(NBFY) = NEX*NDY + I
CCCCCCCCCCC BVX(NBFX) = U0
BVY(NBFY) = V0
29 CONTINUE
DO 28 I = 1 , 32
NBFX = NBFX + 1
NBFY = NBFY + 1
IBNDFX(NBFX) = ICR(1,I)
IBNDFY(NBFY) = ICR(1,I)
BVX(NBFX) = 0.
BVY(NBFY) = 0.
28 CONTINUE
C======================================================================
C-------------- B.C. FOR STREAM FUNCTION -------------
NBS = 0
C-------- FACE OF +Y
Q = 0.
DO I = NDX, 1 , -1
NBS = NBS + 1
IBNDS(NBS) = I*NDY
BVS(NBS) = Q
END DO
C------- FACE OF -X
Q = UAVE*(YCOORD(NEY)-YCOORD(NDY))/2.
IBNDS(NBS) = NEY
BVS(NBS) = Q
DO J = NEY-1 , 2 , -1
NBS = NBS + 1
IBNDS(NBS) = J
DDY = YCOORD(J)-YCOORD(J+1)
Q = Q + DDY * UAVE
BVS(NBS) = Q
END DO
Q = Q + UAVE*(YCOORD(1)-YCOORD(2))/2.
C------- FACE OF -Y
DO I = 1 , NDX
NBS = NBS + 1
IBNDS(NBS) = 1 + NDY*(I-1)
BVS(NBS) = Q
END DO
C------- FACE OF +X -------------
C>>>>>>> D(STREAM)/DN = 0.
C------- ON THE CYLINDER --------
DO 27 I = 1 , 32
NBS = NBS + 1
IBNDS(NBS) = ICR(1,I)
BVS(NBS) = 0.5*Q
27 CONTINUE
C======================================================================
CALL OPTIM4(NODEX,IREL,INOD,NEW,JNT,NJNT,XCOORD,YCOORD,XTEMP,
* YTEMP,MXE,MXN,NE,NNODE,6 )
DO I = 1 , NBFX
IBNDFX(I) = NEW(IBNDFX(I))
END DO
DO I = 1 , NBFY
IBNDFY(I) = NEW(IBNDFY(I))
END DO
DO I = 1 , NBS
IBNDS(I) = NEW(IBNDS(I))
END DO
C======================== INITIAL CONDITIONS ==========================
U0 = UAVE
V0 = 0.
XCC = (IFRONT+4)*DX
YCC = TLY / 2.
DO 60 I = 1 , NNODE
U(I) = U0
V(I) = V0
IF((XCOORD(I)-XCC)**2+(YCOORD(I)-YCC)**2
* .LE. 1.2**2) U(I) = 0.
IF(XCOORD(I).GE. XCC .AND. YCOORD(I) .GE.(YCC-0.6) .AND.
* YCOORD(I).LE.(YCC+0.6) ) U(I) = 0.
60 CONTINUE
C=======================================================================
C MAKING DATA FILES
C---------- 'PROJECT'.DAT
OPEN ( 1, FILE=FNAME(1), STATUS = 'UNKNOWN' )
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---------- NODE
WRITE(1,*) NNODE
DO I = 1 , NNODE
WRITE(1,*) I, XCOORD(I), YCOORD(I)
END DO
C---------- BOUNDARY
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=======================================================================
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
C
C
SUBROUTINE OPTIM4 ( ICON,IREL, INOD, NEW, JNT, NJNT, X, Y,
& XTEMP, YTEMP, MAXEL, MAXNOD, NEL, NNOD, IW )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION ICON(MAXEL,4), IREL(MAXNOD,10), INOD(MAXNOD),
& NEW(MAXNOD),JNT(MAXNOD),NJNT(MAXNOD),X(MAXNOD),Y(MAXNOD),
& XTEMP(MAXNOD), YTEMP(MAXNOD),JTEMP(10), ITEMP(4)
C
CALL SETUP( ICON,IREL,INOD,MAXEL,MAXNOD,NEL,NNOD,IBWOR )
IMP = 0
CALL OPTIM(IREL,INOD,JNT,NJNT,NEW,MAXNOD,NNOD,IBWOR,MAXDIF,IMP)
IBW = MAXDIF
DO 10 I = 1, NNOD
JEND = INOD(I)
DO 20 J = 1, JEND
JR = JEND - J + 1
JTEMP(J) = IREL(I,JR)
20 CONTINUE
DO 10 K = 1, JEND
IREL(I,K) = JTEMP(K)
10 CONTINUE
CALL OPTIM(IREL,INOD,JNT,NJNT,NEW,MAXNOD,NNOD,IBW,MAXDIF,IMP)
IF ( IMP .LE. 0 ) GO TO 90
IBOLD = IBWOR + 1
MAXBWH = MAXDIF + 1
WRITE(IW,200) IBOLD, MAXBWH
DO 30 I = 1, NEL
DO 40 J = 1, 4
ITEMP(J) = ICON(I,J)
40 CONTINUE
DO 30 J = 1, 4
ICON(I,J) = NEW( ITEMP(J) )
30 CONTINUE
DO 50 I = 1, NNOD
XTEMP(I) = X(I)
YTEMP(I) = Y(I)
50 CONTINUE
DO 60 I = 1, NNOD
N = NEW(I)
X(N) = XTEMP(I)
Y(N) = YTEMP(I)
60 CONTINUE
C WRITE(IW,300)
DO 70 I = 1, NNOD
C WRITE(IW,400) I, NEW(I)
70 CONTINUE
GO TO 100
90 IBOLD = IBWOR + 1
WRITE(IW,500) IBOLD
100 RETURN
200 FORMAT ( /,6X,'HALF-BANDWIDTH HAS BEEN REDUCED FROM',I4,' TO ',
& I4,/ )
300 FORMAT ( /,6X,'NEW OPTIMIZED NODE NUMBERING:',//,
& 10X,'OLD NODE NO.',10X,'NEW NODE NO.',/ )
400 FORMAT ( 14X,I4,4X,'----->',4X,I4 )
500 FORMAT ( 6X,'NODE NUMBERING CANNOT BE IMPROVED.',//,
& 6X,'HALF-BANDWIDTH REMAINS AT ',I4 )
END
C
C
SUBROUTINE SETUP ( ICON,IREL,INOD,MAXEL,MAXNOD,NEL,NNOD,IBWOR )
DIMENSION ICON(MAXEL,4), IREL(MAXNOD,10), INOD(MAXNOD), KPERM(4,3)
DATA KPERM / 2, 3, 4, 1, 3, 4, 1, 2, 4, 1, 2, 3 /
IBWOR = 0
DO 100 I = 1, NNOD
KOUNT = 0
DO 200 J = 1, NEL
DO 300 K = 1, 4
NOD1 = ICON(J,K)
NOD2 = ICON(J,KPERM(K,1))
NOD3 = ICON(J,KPERM(K,2))
NOD4 = ICON(J,KPERM(K,3))
IBW = MAX0 ( IABS(NOD1-NOD2), IABS(NOD1-NOD3), IABS(NOD1-NOD4),
& IABS(NOD2-NOD3), IABS(NOD2-NOD4), IABS(NOD3-NOD4) )
IF ( I .NE. NOD1 ) GO TO 300
IF ( KOUNT .EQ. 0 ) GO TO 11
DO 10 L = 1, KOUNT
IF ( IREL(I,L) .EQ. NOD2 ) GO TO 12
10 CONTINUE
11 KOUNT = KOUNT + 1
IREL(I,KOUNT) = NOD2
12 DO 13 L = 1, KOUNT
IF ( IREL(I,L) .EQ. NOD3 ) GO TO 14
13 CONTINUE
KOUNT = KOUNT + 1
IREL(I,KOUNT) = NOD3
14 DO 15 L = 1, KOUNT
IF ( IREL(I,L) .EQ. NOD4 ) GO TO 300
15 CONTINUE
KOUNT = KOUNT + 1
IREL(I,KOUNT) = NOD4
300 CONTINUE
IBWOR = MAX0 ( IBW, IBWOR )
200 CONTINUE
INOD(I) = KOUNT
100 CONTINUE
RETURN
END
C
C
SUBROUTINE OPTIM ( IREL, INOD, JNT, NJNT, NEW, MAXNOD, NNOD, IDIFF
& ,MAXDIF, IMP )
DIMENSION IREL(MAXNOD,10),INOD(MAXNOD),JNT(MAXNOD),NJNT(MAXNOD),
& NEW(MAXNOD)
MAXDIF = IDIFF
DO 60 IX = 1, NNOD
DO 20 J = 1, NNOD
JNT(J) = 0
NJNT(J) = 0
20 CONTINUE
MAX = 0
I = 1
NJNT(1) = IX
JNT(IX) = 1
K = 1
30 JEND = INOD(NJNT(I))
DO 40 JJ = 1, JEND
KX = IREL(NJNT(I),JJ)
IF ( JNT(KX) .GT. 0 ) GO TO 40
K = K + 1
NJNT(K) = KX
JNT(KX) = K
NDIFF = IABS ( I - K )
IF ( NDIFF .GE. MAXDIF ) GO TO 60
MAX = MAX0 ( NDIFF, MAX )
40 CONTINUE
IF ( K .EQ. NNOD ) GO TO 50
I = I + 1
GO TO 30
50 MAXDIF = MAX
DO 70 J = 1, NNOD
NEW(J) = JNT(J)
70 CONTINUE
IMP = IMP + 1
60 CONTINUE
RETURN
END