PROGRAM SETHL4
C========= WAVE AROUND CYLINDER ========
C=======================================================================
C DATA GENERATING PROGRAM FOR HELM4Q.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------------- DIAMETER OF CYLINDER = 1. -------------------------------
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=4, MXE=20000, MXN=30000, MXB=3000 )
PARAMETER ( WAVEH=1.,ALPHA=10.D0,
* ISIDE = 5, IFRONT = 10, IBACK = 20, NEX = IFRONT + 8 + IBACK,
* NEY = ISIDE + 8 + ISIDE, DX = 0.4, DY = 0.4,
* NDX = NEX + 1, NDY = NEY + 1 )
C======================================================================
C
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
* IBNDF(MXB),BV(MXB), ITYPE(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 PROJECT*12,EXFILE*3
LOGICAL YES
C=======================================================================
DATA PROJECT / 'HEL4DATA.QQQ' /
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=======================================================================
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=======================================================================
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======================================================================
C--------- BOUNDARY CONDITIONS ---------
C............ MOMENTUM EQUATIONS .......
C FACE OF -X
NBF = 0
DO J = 1 , NDY
NBF = NBF + 1
IBNDF(NBF) = J
BV(NBF) = WAVEH
ITYPE(NBF) = 1
END DO
C------- FACE OF +X
V0 = 0.
DO 29 I = 1 , NDY
NBF = NBF + 1
IBNDF(NBF) = NEX*NDY + I
BV(NBF) = 0.
ITYPE(NBF) = 1
29 CONTINUE
C======================================================================
CALL OPTIM4(NODEX,IREL,INOD,NEW,JNT,NJNT,XCOORD,YCOORD,XTEMP,
* YTEMP,MXE,MXN,NE,NNODE,6 )
DO I = 1 , NBF
IBNDF(I) = NEW(IBNDF(I))
END DO
C=======================================================================
C DATA FILE INQUIRY
EXFILE = 'NEW'
INQUIRE ( FILE = PROJECT, EXIST = YES )
IF ( YES ) EXFILE='OLD'
C=======================================================================
C MAKING DATA FILES
OPEN ( 1, FILE=PROJECT, STATUS = EXFILE )
C------ CONTROL PARAMETERS
WRITE(1,'(F10.5)') ALPHA
WRITE(1,'(I5)' ) NE
C------ ELEMENTS
DO I = 1 , NE
WRITE (1,'(5I5)') I, (NODEX(I,J), J = 1 , ND )
END DO
C------ NODAL COORDINATES
WRITE(1,'(I5)' ) NNODE
WRITE(1,120) ( I,XCOORD(I), YCOORD(I) , I = 1, NNODE )
120 FORMAT ( I5, 2F10.6 )
C------ BOUNDARY CONDITIONS
WRITE(1,'(I5)' ) NBF
WRITE (1,125) ( IBNDF(I), ITYPE(I),BV(I) , I = 1 , NBF )
125 FORMAT ( 2I5, F10.5 )
CLOSE (1)
STOP
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