戻る

C ******************************************************************** C
C C
C C V D F 2 D-FEM C
C C
C Two-Dimensional Transient Convection-Diffusion C
C C
C Equation Solver by Modified Gakerkin Method C
C C
C ( V.4.0 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
PROGRAM MAIN
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, NBT
C
DIMENSION X ( 1000), Y ( 1000), DMT( 1000,5),
1 PMT( 1000,5), TDP( 1000,10), BT ( 1000,10),
2 USM( 1000,1000), DEL( 1000,6), PEL( 1000,6),
3 BE ( 1000,3), CE ( 1000,3), EVC( 1000),
4 RVC( 1000), AT ( 1000), Q ( 1000),
5 C ( 1000), CO ( 1000), NDE( 1000,4),
6 NBT( 1000)
C
COMMON /ABC / AR( 1000)
COMMON /BDC / CB( 20), NTAB( 20), NBC
COMMON /CNST/ CFDX( 5), CFDY( 5), DT, LMT, U, V, RK
C
OPEN ( UNIT=7, FILE='CVDF2D_OUT.DATA' )
OPEN ( UNIT=8, FILE='CVDF2D_ANL.DATA' )
C
CALL TTLE
C
900 CONTINUE
C
CALL SBIP ( RL, RR, NX, NNP, NEL, ICR )
C
CALL STND ( X, Y, RL, RR, NX, NDE, NBT, NNP, NEL, MB )
C
TME = 0.D0
LOP = 0
C
CALL INIT ( X, Y, NX, NDE, C, CO, Q, RL, RR, NBT, NNP,
1 NEL, MB, MW, ICR )
C
C ----- Construction of Matrix ( D), ( P), ( DP) ---
C
CALL MTCN ( X, Y, DMT, PMT, TDP, DEL, PEL, AT, BE, CE,
1 NDE, NNP, NEL, MB, MW )
C
999 TME = TME + DT
LOP = LOP + 1
C
C ----- Construction of Roght-hand side vector ---
C
CALL VCON ( PMT, CO, EVC, RVC, Q, AT, NDE, NNP, NEL, MB )
C
C ----- Solution ---
C
CALL USLV ( TDP, USM, RVC, C, BT, NBT, NNP, MW, MB, LOP )
C
CALL POUT ( LOP, TME, C, CO, NBT, NNP, RXT, NPT )
C
CALL CANL ( LOP, NX, TME, NNP )
C
CALL ECHK ( NX, C, NNP, LOP )
C
CALL FNAL ( AR, C, LOP, TME, RXT, RR, NPT, NX, NNP,
1 *999, *900 )
C
STOP
END
C **********************************************************************
C
SUBROUTINE FNAL ( AR, C, LOP, TME, RXT, RR, NPT, NX, NNP, *, * )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CNST/ CFDX( 5), CFDY( 5), DT, LMT, U, V, RK
C
DIMENSION AR( NNP), C( NNP)
C
IF ( LOP.LE.LMT) RETURN 1
C
WRITE(7,2000) LOP, TME, RXT, NPT, RR
2000 FORMAT(I3,2F7.3,I3,F7.3)
WRITE(7,2100) ( C( I), I=1, NNP )
2100 FORMAT(10F7.3)
WRITE(8,2100) ( AR( I), I=1, NX )
C
WRITE(6,2200)
2200 FORMAT(/,' Computation to be continued ? ( 1: Yes 2: No )')
** READ(5,*) ICP
ICP = 2
C
IF ( ICP.EQ.1) RETURN 2
C
RETURN
END
C *********************************************************************
C
SUBROUTINE SBIP ( RL, RR, NX, NNP, NEL, ICR )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CNST/ CFDX( 5), CFDY( 5), DT, LMT, U, V, RK
C
C ----- Attention : Computational Area: ( RL ) X ( RL* RR ) ---
C
RL = 6.D0
RR = 0.1D0
C -------------------------------------------------------------
C
WRITE (6,2000)
2000 FORMAT(/,' NX = ? ')
** READ (5,*) NX
NX = 21
C ************
C
WRITE (6,2100)
2100 FORMAT(' * Kx = ? ')
** READ (5,*) CFDX( 1)
CFDX( 1) = 0.1D0
C
WRITE (6,2200)
2200 FORMAT(' * Ky = ? ')
** READ (5,*) CFDY( 1)
CFDY( 1) = 0.1D0
C
WRITE(6,2300)
2300 FORMAT(' Time Step = ? ')
** READ(5,*) DT
C
DT = 0.2D0
C ***************
C
WRITE(6,2400)
2400 FORMAT(' * U = ? ')
** READ(5,*) U
U = 0.3D0
C
WRITE(6,2500)
2500 FORMAT(' * V = ? ')
** READ(5,*) V
V = 0.D0
C
C ----- Attention ---
C
LMT = 7.D0/ DT
C ===================
C
** WRITE(6,*) ' * NX=', NX,' U=', U,' V=', V,' CFDX(1)=', CFDX(1),
** 1 ' CFDY(1)=', CFDY(1), ' DT=', DT,' LMT=', LMT
C
C ----- Without Correction ( 0) / With Correction ( 1) ---
C
ICR = 1
C
NNP = NX* 3
NEL = 4*( NX-1)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CANL ( LOP, NX, TME, NNP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /ABC / AR( 1000)
COMMON /BDC / CB( 20), NTAB( 20), NBC
COMMON /CNST/ CFDX( 5), CFDY( 5), DT, LMT, U, V, RK
C
RL = 6.D0/ DBLE( NX-1)
C
C ------------------------------------------------
DO 100 I = 1, NX
C
AR( I) = 0.D0
C
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, NX
C
X = RL* DFLOAT( I-1)
C
CF1 = 0.5D0* ( X - U* TME )/ DSQRT( RK* TME )
CF1 = ERFC( CF1)
CF2 = DEXP( X*U/ RK )
C
CF3 = 0.5D0* ( X + U* TME )/ DSQRT( RK* TME )
CF3 = ERFC( CF3)
C
AR( I) = 0.5D0* ( CF1 + CF2* CF3 )* CB( 1)
C
IF ( AR( I).LT.0.00001D0) GO TO 300
C
200 CONTINUE
C ------------------------------------------------
C
300 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,100)
WRITE(6,101)'C ************************************************ C'
WRITE(6,101)'C C'
WRITE(6,101)'C CVDF2D - FEM C'
WRITE(6,101)'C C'
WRITE(6,101)'C ( V.4.0 ) C'
WRITE(6,101)'C C'
WRITE(6,101)'C 2-D Convection-Diffusion Equation Solver C'
WRITE(6,101)'C C'
WRITE(6,101)'C by F.E.M. C'
WRITE(6,101)'C C'
WRITE(6,101)'C ( MGM & GM ) C'
WRITE(6,101)'C C'
WRITE(6,101)'C Copyright 2011 : Yasuhiro MATSUDA C'
WRITE(6,101)'C C'
WRITE(6,101)'C ************************************************ C'
C
100 FORMAT(/)
101 FORMAT(12X,A55)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ECHK ( NX, C, NNP, LOP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION C( NNP)
C
COMMON /ABC / AR( 1000)
COMMON /CNST/ CFDX( 5), CFDY( 5), DT, LMT, U, V, RK
C
C ----- Total Mean Error ---
C
SUM = 0.D0
AMX = 0.D0
K = 1
C
IF ( LOP.EQ.1) TSM = 0.D0
C
C ------------------------------------------------
DO 100 I = 1, NX
C
IF ( AR( I).GT.AMX ) AMX = AR( I)
C
100 CONTINUE
C ------------------------------------------------
C
DO 200 I = 5, NNP, 3 ! Except value at left boundary node
C *****
C
IF ( C( I).GT.0.95D0 .AND. C(I).LT.0.05D0 )
C ***********************************************
1 GO TO 200 ! Attention
C
SUM = SUM + ( DABS( AR( K) - C( I)))/ AMX
K = K + 1
C
200 CONTINUE
C ------------------------------------------------
C
** WRITE(6,*) ' K, AMX=', K, AMX
C
SUM = SUM/ DFLOAT( K-1) ! Per Node
TSM = TSM + SUM
C
C --------------------------------------
IF ( LOP.EQ.LMT) THEN
TSM = TSM/ DFLOAT( LMT)
C
* WRITE(6,2000) TSM* 100.D0, SUM* 100.D0
* 2000 FORMAT(/,' *** Temp. Av. Error =',F8.5,'(%)',' * Av. Error at Last
* 1 Loop =',F8.5,'(%)')
C
WRITE(6,2000) SUM* 100.D0
2000 FORMAT(/,' ** Av. Error at Last Loop =',F7.3,'(%)')
END IF
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
FUNCTION ERFC ( A )
C
IMPLICIT REAL*8(A-H,O-Z)
C
L = 100
DERF = 0.D0
DH = DBLE( A / DFLOAT( L))
C
C ----------------------------------------------------------------
DO 100 I = 1, 100
C
DU = DH* DBLE( I)
DV = DH* DBLE( I-1)
C
DERF = DERF + ( DEXP( - DU* DU) + DEXP( - DV* DV))/ 2.D0* DH
C
100 CONTINUE
C ----------------------------------------------------------------
C
ERFC = DERF* 2.D0 / DSQRT( 4.D0* ATAN( 1.D0))
ERFC = 1.D0 - ERFC
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GDPM ( TDP, DPMT, NNP, MB, MW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- Construction of 'TDP'-Matrix from 'DPMT' ----
C
DIMENSION TDP( NNP, MW), DPMT( NNP, MB)
C
C --------------------------------------
C
CALL ZARY ( NNP, MW, TDP )
C
C --------------------------------------
DO 100 I = 1, MB
C
C --------------------------------------
DO 120 J = 1, I
C
TDP( I,J) = DPMT( J, I-J+1 )
C
120 CONTINUE
C --------------------------------------
C
JS = I + 1
JE = MB + I - 1
C --------------------------------------
DO 140 J = JS, JE
C
TDP( I,J) = DPMT( I, J-I+1 )
C
140 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
MBP1 = MB + 1
MBM1 = MB - 1
C
C -----------------------------------------------
DO 200 I = MBP1, NNP
C
C --------------------------------------
DO 220 J = 1, MBM1
C
TDP( I,J) = DPMT( I+J-MB, MB-J+1 )
C
220 CONTINUE
C --------------------------------------
DO 240 J = MB, MW
C
TDP( I,J) = DPMT( I, J-MB+1 )
C
240 CONTINUE
C --------------------------------------
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE GLMAT ( A, B, NDE, NNP, NEL, MB )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- From Element Matrix to Global Matrix ---
C
INTEGER*2 NDE
C
DIMENSION A( NNP, MB), B( NEL, 6), NDE( 1000, 4), LLL( 3, 3)
C
C ---------------------------------
C
CALL ZARY ( NNP, MB, A )
C
C ---------------------------------
C
LLL( 1,1) = 1
LLL( 1,2) = 2
LLL( 1,3) = 3
C
LLL( 2,1) = 2
LLL( 2,2) = 4
LLL( 2,3) = 5
C
LLL( 3,1) = 3
LLL( 3,2) = 5
LLL( 3,3) = 6
C
C ----------------------------------------------------------
DO 100 NE = 1, NEL
DO 100 I = 1, 3
C ----------------------------------------------------------
C
IG = NDE ( NE,I)
C ----------------------------------------------------------
DO 100 J = I, 3
C
JG = NDE ( NE,J)
IF ( IG - JG) 120, 140, 140
120 A ( IG, JG-IG+1 ) = A ( IG, JG-IG+1 ) + B ( NE, LLL( I,J))
GO TO 100
140 A ( JG, IG-JG+1 ) = A ( JG, IG-JG+1 ) + B ( NE, LLL( I,J))
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GSOL ( A, B, N, MB, MW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C To Solve Matrix A( N, MW) by Gaussian Elimination Method
C
C * * * *
C * * * * *
C * * * * * *
C * * * * * * *
C * * * * * * *
C * * * * * * * Unsymmetric Banded Matrix
C N * ---MW---- *
C * * * * * * * A ( N, MW )
C * * * * * * *
C * * * * * * *
C * * * * * *
C * * * * *
C * * * *
C --MB---
C
DIMENSION A( N, MW), B( N)
C
N1 = N - 1
C
C ------------------------------------------------
DO 100 K = 1, N1
C
IF ( K.LE.( N-MB)) THEN
KMX = MB - 1
ELSE
KMX = MB - ( K - ( N-MB))
END IF
C
K1 = K + 1
IF ( K.LE.MB ) THEN
AKK = 1.D0/ A( K,K)
ELSE
AKK = 1.D0/ A( K,MB)
END IF
C
BK = B( K)
KMX1 = K + MB - 1
C
IF ( KMX1.GT.N ) KMX1 = N
C
C --------------------------------------
DO 110 I = K1, KMX1
C
IF ( I.GT.MB ) THEN
KK = K - ( I-MB)
ELSE
KK = K
END IF
C
AIK = A( I,KK)* AKK
B( I) = B( I) - AIK* BK
C
IF ( I.LE.MB ) THEN
C ----------------------------
DO 120 J = K, KMX + K
C
A( I,J) = A( I,J) - AIK * A( K,J)
120 CONTINUE
C ----------------------------
END IF
C
IF ( I.GT.MB ) THEN
C
C ----------------------------
DO 130 J = KK, KMX + KK
C
IF ( K.LT.MB ) THEN
A( I,J) = A( I,J) - AIK * A( K,K+J-KK)
ELSE
A( I,J) = A( I,J) - AIK * A( K,MB+J-KK)
END IF
C
130 CONTINUE
C ----------------------------
C
END IF
C
110 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
C ----- Backward Substition ---
C
K = N
B( K) = B( K) / A( K,MB)
C
500 K = K - 1
IF ( K.LE.0 ) GO TO 999
K1 = K + 1
SUM = 0.D0
C
C ------------------------------------------------
IF ( K.GE.MB ) THEN
MM = MIN0 ( N-K+MB, MW )
C
C ----------------------------
DO 520 J = MB + 1, MM
C
K2 = K1 + J - MB - 1
SUM = SUM + A( K, J)* B( K2)
C
520 CONTINUE
C ----------------------------
C
B( K) = ( B( K) - SUM ) / A( K, MB)
C
ELSE
C ----------------------------
DO 540 J = K1, MB + K - 1
C
SUM = SUM + A( K,J)* B( J)
540 CONTINUE
C ----------------------------
B( K) = ( B( K) - SUM )/ A( K, K)
C
END IF
C ------------------------------------------------
C
GOTO 500
C
999 CONTINUE
RETURN
END
C **********************************************************************
C
SUBROUTINE INIT ( X, Y, NX, NDE, C, CO, Q, RL, RR, NBT,
1 NNP, NEL, MB, MW, ICR )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, NBT
C
DIMENSION X ( NNP), Y ( NNP), NDE( 1000,4), C( NNP),
1 CO( NNP), NBT( NNP), Q ( NEL)
C
COMMON /BDC / CB( 20), NTAB( 20), NBC
COMMON /CNST/ CFDX( 5), CFDY( 5), DT, LMT, U, V, RK
C
MB = 0
C --------------------------------------
DO 100 NE = 1, NEL
DO 100 I = 1, 3
C --------------------------------------
C
II = NDE ( NE, I)
C
C --------------------------------------
DO 100 J = 1, 3
C
JJ = NDE ( NE, J)
IJ = JJ - II + 1
IF ( IJ.GT.MB ) MB = IJ
C
100 CONTINUE
C --------------------------------------
C
WRITE(6,2000) MB
2000 FORMAT(/,' * Size of Band Matrix : MB =',I3)
C
CALL ZVEC ( NEL, Q )
C
CALL ZVEC ( NNP, CO )
C
CALL ZVEC ( NNP, C )
C
CALL ZVI2 ( NNP, NBT )
C
NBC = 3
C
C --------------------------------------
DO 200 I = 1, NBC
C
NTAB( I) = I
CB ( I) = 1.D0
C
200 CONTINUE
C --------------------------------------
C
WRITE(6,2100) ( I, NTAB( I), CB( I), I = 1, NBC )
2100 FORMAT(/,3(' NO. BD-NDE ''C'' ')/( 3( 2I6,F8.3,2X)))
C
C --------------------------------------
DO 300 I = 1, NBC
C
C ( NTAB( I)) = CB( I)
CO ( NTAB( I)) = CB( I)
NBT( NTAB( I)) = 9999
C
300 CONTINUE
C --------------------------------------
C
C ----- Corrected Diffusivity --- ( Attention : Depending mesh shapes )
C
RK = CFDX( 1)
C
IF ( ICR.EQ.0) GO TO 500
C
C ----- Correction for Diffusivity ---------------
C
C ----- Computational mesh size ---
C
RRX = RL/ DFLOAT( NX-1)
RRY = ( RL* RR)* 0.5D0
C ---------------------------------
C
** WRITE(6,*) ' * RRX,RRY', RRX, RRY
C
BX = ( U* DT)/ RRX
BY = ( V* DT)/ RRY
C
RX = ( CFDX( 1)* DT)/( RRX* RRX)
RY = ( CFDY( 1)* DT)/( RRY* RRY)
C
C ----- Attention : Depends on a mesh shape ------
C
*** CFDX( 1) = CFDX( 1)
*** 1 * ( 1.D0 + ( BX + BY)** 2/( 2.D0*( RX + RY)) )
C
*** CFDY( 1) = CFDY( 1)
*** 1 * ( 1.D0 + ( BX + BY)** 2/( 2.D0*( RX + RY)) )
C
C ----------------------------------------------------
C
CFDX( 1) = CFDX( 1)
1 * ( 1.D0 + ( BX + BY)** 2/( 2.D0*( RX + RY))
2 - ( BX + BY)/( 6.D0*( RX + RY) ) )
C
CFDY( 1) = CFDY( 1)
1 * ( 1.D0 + ( BX + BY)** 2/( 2.D0*( RX + RY))
2 - ( BX + BY)/( 6.D0*( RX + RY) ) )
C
C ------------------------------------------------
C
WRITE(6,2600) U, V, DT, RRX, RRY
2600 FORMAT(/,' ** U =',F6.3,' V =',F6.3,' DT =',F6.3,' RRX =',
1 F7.3,' RRY =',F7.3)
C
WRITE(6,2500) BX, BY, RX, RY
2500 FORMAT(' ** Bx =',F9.5,' By =',F9.5,' Rx =',F9.5,
1 ' Ry =',F9.5)
C
WRITE(6,2200) CFDX( 1), CFDY( 1)
2200 FORMAT(/,' ** New Diffusion : CFDX =',F9.5,' CFDY =',F9.5,
1 ' ***')
C ------------------------------------------------
C
500 CONTINUE
WRITE(6,2300) CFDX( 1), CFDY( 1), DT, LMT, U
2300 FORMAT(/,' * CFDX(1) =',F7.3,' * CFDY(1) =',F7.3,' * DT=',
1 F7.3,' * LMT =',I4,' * U =',F6.3,/)
C
WRITE(7,2400) NX, RRX, U, CFDX( 1), CFDY( 1), DT, ICR, LMT+1
2400 FORMAT(I3, 5F7.3, 2I3)
C
C -------------------
C
MW = 2* MB - 1
C -------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE MTCN ( X, Y, DMT, PMT, TDP, DEL, PEL, AT,
1 BE, CE, NDE, NNP, NEL, MB, MW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- Construction of Matrix 'DMT', 'PMT' and 'TDP' ---
C
INTEGER*2 NDE
C
DIMENSION X ( NNP), Y ( NNP), DMT( NNP,MB),
1 PMT( NNP,MB), TDP( NNP,MW), DEL( NEL,6),
2 PEL( NEL,6), AT ( NEL), BE ( NEL,3),
3 CE ( NEL,3), NDE( 1000,4), XL ( 3),
4 YL ( 3), A ( 3), B ( 3),
5 C ( 3), LLL( 3,3)
C
COMMON /CNST/ CFDX( 5), CFDY( 5), DT, LMT, U, V, RK
C
KYR = 0
C
C ----------------------------------------------------------
DO 100 NE = 1, NEL
C
C ------------------------------------------------
DO 110 K = 1, 3
C
XL( K) = X ( NDE( NE,K))
YL( K) = Y ( NDE( NE,K))
C
110 CONTINUE
C ------------------------------------------------
C
GX = ( XL( 1) + XL( 2) + XL( 3))/ 3.D0
GY = ( YL( 1) + YL( 2) + YL( 3))/ 3.D0
C
C ------------------------------------------------
DO 120 J = 1, 3
C
XL( J) = XL( J) - GX
YL( J) = YL( J) - GY
C
120 CONTINUE
C -----------------------------------------------
C
A( 1) = XL( 2)*YL( 3) - XL( 3)*YL( 2)
A( 2) = XL( 3)*YL( 1) - XL( 1)*YL( 3)
A( 3) = XL( 1)*YL( 2) - XL( 2)*YL( 1)
C
B( 1) = YL( 2) - YL( 3)
B( 2) = YL( 3) - YL( 1)
B( 3) = YL( 1) - YL( 2)
C
C( 1) = XL( 3) - XL( 2)
C( 2) = XL( 1) - XL( 3)
C( 3) = XL( 2) - XL( 1)
C
C ------------------------------------------------
DO 130 L = 1, 3
C
BE( NE,L) = B(L)
CE( NE,L) = C(L)
C
130 CONTINUE
C ------------------------------------------------
C
AT( NE) = ( A( 1) + A( 2) + A( 3))/ 2.D0
C
IF ( AT( NE).LE.0.D0 ) WRITE(6,2000) NE
2000 FORMAT(/,' *** Error ( AREA.LE.0. )--- NE =',I4 )
C
IF ( AT( NE).LE.0.D0 ) KYR = 1
C
A4 = 4.0D0* AT( NE)
C
CDXM = CFDX( NDE( NE,4))/ A4
CDYM = CFDY( NDE( NE,4))/ A4
C
LLL( 1,1) = 1
LLL( 1,2) = 2
LLL( 1,3) = 3
C
LLL( 2,1) = 2
LLL( 2,2) = 4
LLL( 2,3) = 5
C
LLL( 3,1) = 3
LLL( 3,2) = 5
LLL( 3,3) = 6
C
C ------------------------------------------------
DO 140 I = 1, 3
DO 140 J = 1, 3
C
DEL( NE, LLL( I,J)) = CDXM* B( I)* B( J) + CDYM* C( I)* C( J)
C
140 CONTINUE
C ------------------------------------------------
C
PEL( NE,1) = AT( NE)/ 6.D0
PEL( NE,2) = AT( NE)/ 12.D0
PEL( NE,3) = AT( NE)/ 12.D0
PEL( NE,4) = AT( NE)/ 6.D0
PEL( NE,5) = AT( NE)/ 12.D0
PEL( NE,6) = AT( NE)/ 6.D0
C
100 CONTINUE
C
IF ( KYR.EQ.1) STOP
C
C ------------------------------------------------
C
CALL GLMAT ( DMT, DEL, NDE, NNP, NEL, MB )
C
C ------------------------------------------------
C
CALL GLMAT ( PMT, PEL, NDE, NNP, NEL, MB )
C
C ------------------------------------------------
DO 200 I = 1, NNP
DO 200 J = 1, MB
C
DMT( I,J) = PMT( I,J) + DT* DMT( I,J)
C
200 CONTINUE
C ------------------------------------------------
C
CALL GDPM ( TDP, DMT, NNP, MB, MW )
C
C ------------------------------------------------
C
CALL MDTD ( TDP, BE, NDE, NNP, NEL, MB, MW )
C
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MDTD ( TDP, BE, NDE, NNP, NEL, MB, MW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- Modification of 'TDP' considering 'U* DC/DX' ---
C
INTEGER*2 NDE
C
DIMENSION TDP( NNP, MW), BE ( NEL, 3), NDE ( 1000, 4)
C
COMMON /CNST/ CFDX( 5), CFDY( 5), DT, LMT, U, V, RK
C
C ----------------------------------------------------------
DO 100 NE = 1, NEL
C
C -----------------------------------------------
DO 150 I = 1, 3
C
IG = NDE ( NE, I)
C
C ------------------------------------------------
DO 150 J = 1, 3
C
JG = NDE ( NE, J)
C
IF ( IG.LE.MB) TDP( IG,JG) = TDP( IG,JG)
1 + DT * U * BE( NE,J) / 6.D0
IF ( IG.LE.MB) GO TO 150
JT = JG - IG + MB
C
TDP( IG,JT) = TDP( IG,JT) + DT * U * BE( NE,J) / 6.D0
C
150 CONTINUE
C ------------------------------------------------
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE POUT ( LOP, TME, C, CO, NBC, NNP, RXT, N )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C---- Stationarity Check & Print ---
C
INTEGER*2 NBC
C
DIMENSION C( NNP), CO( NNP), NBC( NNP)
C
RXT = 0.D0
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
IF ( NBC( I).EQ.9999 ) GO TO 100
IF ( C( I).LT.0.000001D0 ) GO TO 100
C
DRC = ( C( I) - CO( I))/ C( I)
IF ( DABS( DRC).LT.RXT ) GO TO 100
RXT = DABS( DRC)
N = I
C
100 CONTINUE
C ------------------------------------------------
C
RXT = RXT* 100.D0
C
WRITE(6,2000) LOP, TME, RXT, N
2000 FORMAT(' * LOP=',I5,2X,'( T =',F7.3,')* Var_C =',F7.2,
1 '(%) at NDE =',I4)
C
C ----------------------------
DO 200 I = 1, NNP
C
CO( I) = C( I)
200 CONTINUE
C ----------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE PRUM ( TDP, A, RVC, C, NBC, M, MW, MB, LOP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NBC
C
DIMENSION TDP( M, MW), A ( M, M), RVC( M), C( M), NBC( M)
C
C ------------------------------------------------
DO 100 I = 1, M
C
IF ( NBC( I).EQ.9999 ) GO TO 100
C( I) = RVC( I)
C
100 CONTINUE
C ------------------------------------------------
C
C ----- Fixed Value at Boundary ---
C
IF ( LOP.NE.1 ) GO TO 90
C
C ------------------------------------------------
DO 200 I = 1, M
C
IF ( NBC( I).NE.9999 ) GO TO 200
IF ( I - MB ) 220, 240, 240
C
220 CONTINUE
C ------------------------------------------------
C
TDP( I,I) = TDP( I,I)* 1.D+10
C( I) = C( I)* TDP( I,I)
GOTO 200
C
240 TDP( I,MB) = TDP( I,MB)* 1.D+10
C
C( I) = C( I)* TDP( I,MB)
C
200 CONTINUE
C
GO TO 900
C
90 CONTINUE
C
C ------------------------------------------------
DO 300 I = 1, M
C
IF ( NBC( I).NE.9999 ) GO TO 300
IF ( I - MB ) 320, 340, 340
320 C( I) = C( I)* TDP( I,I)
GO TO 300
C
340 C( I) = C( I)* TDP( I, MB)
C
300 CONTINUE
C ------------------------------------------------
C
900 CONTINUE
C
C ------------------------------------------------
DO 400 I = 1, MB
C
JE = MB + I - 1
C
C ------------------------------------------------
DO 400 J = 1, JE
C
A( I,J) = TDP( I, J)
C
400 CONTINUE
C ------------------------------------------------
DO 500 I = MB + 1, M
C
K = I - MB + 1
JE = MIN0( K + MW - 1, M)
C
C ------------------------------------------------
DO 500 J = K, JE
C
A( I,J) = TDP( I, J-K+1)
500 CONTINUE
C -------------------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE STND ( X, Y, RL, RR, NX, NDE, NBT, NNP, NEL, MB )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE, NBT
C
DIMENSION X( NNP), Y( NNP), NDE( 1000, 4)
C
C ----- Attention : Computational Area: ( RL ) X ( RL* RR ) ---
C
X( 1) = 0.D0
X( 2) = 0.D0
X( 3) = 0.D0
C
C ---------------------------------
DO 100 I = 4, NNP, 3
C
X( I) = X( I-1) + RL/ DBLE( NX-1)
X( I+1) = X( I)
X( I+2) = X( I+1)
C
100 CONTINUE
C ---------------------------------
DO 200 I = 1, NNP, 3
C
Y( I) = RL* RR/ 2.D0
Y( I+1) = 0.D0
Y( I+2) = - RL* RR/ 2.D0
C
200 CONTINUE
C ---------------------------------
C
C ----- Element No. ---
C
NDE( 1,1) = 1
NDE( 1,2) = 2
NDE( 1,3) = 5
C
NDE( 2,1) = 2
NDE( 2,2) = 3
NDE( 2,3) = 5
C
NDE( 3,1) = 1
NDE( 3,2) = 5
NDE( 3,3) = 4
C
NDE( 4,1) = 3
NDE( 4,2) = 6
NDE( 4,3) = 5
C
C ------------------------------------------------
DO 300 I = 5, NEL
DO 300 J = 1, 3
C
NDE( I,J) = NDE( ( I-4),J) + 3
C
300 CONTINUE
C ------------------------------------------------
DO 400 I = 1, NEL
C
NDE( I,4) = 1
C =======================
C
400 CONTINUE
C ------------------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE USLV ( TDP, A, RVC, C, BT, NBC, M, MW, MB, LOP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- To Solve Simultaneous Linear Equations -----
C
INTEGER*2 NBC
C
DIMENSION TDP( M, MW), A( M, M), RVC( M), C( M), NBC( M),
1 BT ( M, MW)
C
C --------------------------------------------------------
C
CALL PRUM ( TDP, A, RVC, C, NBC, M, MW, MB, LOP )
C
C --------------------------------------------------------
DO 100 I = 1, M
DO 100 J = 1, MW
C
BT( I,J) = TDP( I,J)
100 CONTINUE
C --------------------------------------
C
CALL GSOL ( BT, C, M, MB, MW )
C
C --------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE VCON ( A, CO, EVC, RVC, Q, AT, NDE, NNP, NEL, MB )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- 'RVC'- Vector ---
C
INTEGER*2 NDE
C
DIMENSION A( NNP,MB), CO( NNP), EVC( NNP), RVC( NNP),
1 Q( NEL), AT( NEL), NDE( 1000,4)
C
COMMON /CNST/ CFDX( 5), CFDY( 5), DT, LMT, U, V, RK
C
C ---------------------------------
C
CALL ZVEC ( NNP, EVC )
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
IF ( I.EQ.1 ) GO TO 120
JM = I - 1
C
IF ( I.LE.MB ) THEN
JS = 1
GO TO 110
END IF
C
JS = I - MB + 1
110 CONTINUE
C
C --------------------------------------
DO 200 J = JS, JM
C
EVC( I) = EVC( I) + A(J,I-J+1)* CO( J)
200 CONTINUE
C --------------------------------------
C
120 CONTINUE
C
C --------------------------------------
DO 300 J = 1, MB
C
EVC( I) = EVC( I) + A( I,J)* CO(J+I-1)
300 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
CALL ZVEC ( NNP,RVC )
C
C ------------------------------------------------
DO 400 NE = 1, NEL
DO 400 K = 1, 3
C
I = NDE( NE,K)
RVC( I) = RVC( I) - Q( NE)* AT( NE)/ 3.D0
C
400 CONTINUE
C ------------------------------------------------
DO 500 I = 1, NNP
C
RVC( I) = EVC( I) - DT* RVC( I)
500 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ZARY ( M, N, A )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION A( M, N)
C
DO 100 I = 1, M
DO 100 J = 1, N
C
A( I,J) = 0.D0
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ZVEC ( M, B )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION B( M)
C
DO 100 I = 1, M
C
B( I) = 0.D0
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ZVI2 ( M, NC )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NC
C
DIMENSION NC( M)
C
DO 100 I = 1, M
C
NC( I) = 0
100 CONTINUE
C
RETURN
END
C ********************************************************************