戻る

C
PROGRAM MAIN
C ******************************************************************** C
C C
C Error Analysis Software C
C C
C for 2D - Convectin & Diffusion Equation C
C ＝＝ C
C ( for FEMs & FDMs ) C
C C
C ( V.4.0 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CST/ PAI
COMMON /MTV/ MT( 50), TV( 14)
COMMON /NUM/ MTS, NT, NRX, NBX, NKX, NRY, NBY, NKY
COMMON /ABC/ RVX(10), RVY(10), BVX(10), BVY(10), KMX(14), KMY(14)
C
COMMON /OPM/ OPTE( 20,50), OPTT( 20,50), OPTR( 20,50),
1 OPTG( 20,50), OPET( 50), OPT ( 50),
2 OPR ( 50), OPGZ( 50)
C
OPEN ( 7, FILE='CVDF2D_EA.RESULT' )
C
CALL TTLE
C
CALL PCON
C
C -------------------------------------------------------------------
DO 100 NMH = 1, MTS
C
CALL PTTL ( NMH)
C
KT = 0
NTT = NT
IF ( NMH.NE.1) NTT = 1
C
C ----------------------------------------------------------
DO 200 J = 1, NTT
C
T = TV( J)
KT = KT + 1
KR = 0
C
NBXY = NBX* NBY
NKXY = NKX* NKY
C
TSGA = 0.D0
TSGC = 0.D0
C
C ------------------------------------------------
DO 300 IR = 1, NRX
C
RX = RVX( IR)
RY = RVY( IR)
C
SGA = 0.D0
SGC = 0.D0
RMA = 0.D0
C
CALL INIT
C
C --------------------------------------
DO 400 IBX = 1, NBX
C
BX = BVX( IBX)
C
C --------------------------------------
DO 400 IBY = 1, NBY
C
BY = BVY( IBY)
BSA = 0.D0
BSC = 0.D0
RMX = 0.D0
C
C ---------------------------
DO 500 IKX = NKX, 1, -1
C
KX = KMX( IKX)
VAX = 2.D0* PAI/ DFLOAT( KX)
C
C ---------------------------
DO 500 IKY = NKY, 1, -1
C
KY = KMY( IKY)
VAY = 2.D0* PAI/ DFLOAT( KY)
C
CALL SMTD ( NMH, P1, P2, P3, P4, RX, RY, BX, BY, VAX,
1 VAY, T )
C
CALL CMP1 ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY,
1 A, C, DB0, BA, BC, TS, BSA, BSC, RMX, SGA,
2 SGC, RMA )
C
500 CONTINUE
C ---------------------------
C
CALL CMP2 ( NKXY, BSA, BSC, TBS )
C
400 CONTINUE
C --------------------------------------
C
CALL CMP3 ( NBXY, NKXY, SGA, SGC, TSGA, TSGC, TVR )
C
IF ( NMH.EQ.1) CALL COPT1 ( KT, KR, T, RX, RMA, TVR )
C
IF ( NMH.NE.1) CALL OUTP1 ( RX, RY, RMA, SGA, SGC, TVR )
C
300 CONTINUE
C -----------------------------------------------
C
CALL CMP4 ( NRX, NBXY, NKXY, DAB3, TSGA, TSGC, SUMT )
C
200 CONTINUE
C ----------------------------------------------------------
C
IF ( NMH.EQ.1) THEN
C
CALL COPT2 ( KT, KR )
C
CALL OUTPT ( KR ) ! Opt. THETA
C
END IF
C
100 CONTINUE
C --------------------------------------------------------------------
C
WRITE(6,2000)
2000 FORMAT(/10X,' *** END *** ')
C
CLOSE( 7)
C
STOP
END
C ***********************************************************************
C
BLOCK DATA
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
COMMON /NUM/ MTS, NT, NRX, NBX, NKX, NRY, NBY, NKY
COMMON /ABC/ RVX(10), RVY(10), BVX(10), BVY(10), KMX(14), KMY(14)
C
DATA TV / 0.D0, 0.1D0, 0.2D0, 0.3D0, 0.4D0, 0.5D0, 0.6D0, 0.7D0,
1 0.8D0, 0.9D0, 1.D0, 1.1D0, 1.2D0, 1.3D0 /
C
DATA NKX, NKY / 14, 14 /
DATA KMX / -3, -4, -6, -8, - 10, - 20, -40, 3, 4, 6, 8, 10,20,40/
DATA KMY / -3, -4, -6, -8, - 10, - 20, -40, 3, 4, 6, 8, 10,20,40/
C
C ----- Check ---
C
DATA MTS, NT / 12, 14 / ! No. of Methods / No. of THETA
C
DATA NRX, NRY / 6, 6 /
DATA RVX / 0.01D0, 0.02D0, 0.04D0, 0.1D0, 0.2D0, 0.3D0, 0.D0,
1 0.D0, 0.D0, 0.D0 /
DATA RVY / 0.01D0, 0.02D0, 0.04D0, 0.1D0, 0.2D0, 0.3D0, 0.D0,
1 0.D0, 0.D0, 0.D0 /
C
DATA NBX, NBY / 4, 4 /
DATA BVX / 0.05D0, 0.1D0, 0.15D0, 0.2D0, 0.D0, 0.D0, 0.D0, 0.D0,
1 0.D0, 0.D0 /
DATA BVY / 0.05D0, 0.1D0, 0.15D0, 0.2D0, 0.D0, 0.D0, 0.D0, 0.D0,
1 0.D0, 0.D0 /
C
END
C ***********************************************************************
C
SUBROUTINE CMP1 ( P1, P2, P3, P4, RX, RY, BX, BY,
1 VAX, VAY, AA, CC, DB0, BA, BC, TS,
2 BSA, BSC, RMX, SGA, SGC, RMA )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON / CST / PAI
C
AE = ( P1* P3 + P2* P4 )/( P3* P3 + P4* P4 )
BE = ( P2* P3 - P1* P4 )/( P3* P3 + P4* P4 )
AA = DSQRT( AE* AE + BE* BE )
C
IF ( DABS( BE).LT.0.00001D0.AND.AE.GE.0.D0 ) TA = 0.D0
TA = DATAN2( BE, AE )
C
CC = ( TA + VAX* BX + VAY* BY)
IF ( DABS( CC).GT.PAI ) CC = 2.D0* PAI - DABS( CC)
C
CC = CC + 1.D0
C
C ------------------------------
C
VAL = VAX* BX + VAY* BY
C
C ------------------------------
C
AP = DEXP( - VAX* VAX* RX - VAY* VAY* RY)
AE = AP* DCOS( VAX* BX + VAY* BY) - AE
BE = - AP* DSIN( VAX* BX + VAY* BY) - BE
DB0 = DSQRT( AE* AE + BE* BE )
C
IF ( AA.GE.RMX) RMX = AA
C
BSA = BSA + ( AA - AP )* ( AA - AP )
BSC = BSC + ( CC - 1.D0)* ( CC - 1.D0)
C
C -----------------------------------------
C
BA = DSQRT(( AA - AP )* ( AA - AP ))
C
C -----------------------------------------
BC = DSQRT(( CC - 1.D0)* ( CC - 1.D0))
TS = DSQRT( BA* BA + BC* BC )
C
IF ( AA.GE.RMA ) RMA = AA
C
SGA = SGA + ( AA - AP )* ( AA - AP )
SGC = SGC + ( CC - 1.D0)* ( CC - 1.D0)
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE CMP2 ( NKXY, BSA, BSC, TBS )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
BSA = DSQRT( BSA/ DFLOAT( NKXY))
BSC = DSQRT( BSC/ DFLOAT( NKXY))
C
TBS = DSQRT( BSA* BSA + BSC * BSC )
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE CMP3 ( NBXY, NKXY, SGA, SGC, TSGA, TSGC, TVR )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
TSGA = TSGA + SGA
SGA = DSQRT( SGA/( DFLOAT( NBXY)* DFLOAT( NKXY)))
C
TSGC = TSGC + SGC
SGC = DSQRT( SGC/( DFLOAT( NBXY)* DFLOAT( NKXY)))
C
TVR = DSQRT( SGA* SGA + SGC* SGC )
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE CMP4 ( NRX, NBXY, NKXY, DAB3, TSGA, TSGC, SUMT )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DAB3 = DSQRT( DAB3/( DFLOAT(NRX)* DFLOAT(NBXY)* DFLOAT(NKXY)))
C
TSGA = DSQRT( TSGA/( DFLOAT(NRX)* DFLOAT(NBXY)* DFLOAT(NKXY)))
TSGC = DSQRT( TSGC/( DFLOAT(NRX)* DFLOAT(NBXY)* DFLOAT(NKXY)))
C
SUMT = DSQRT( TSGA* TSGA + TSGC* TSGC )
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE COPT1 ( KT, KR, T, RX, RMA, TVR )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /OPM/ OPTE( 20,50), OPTT( 20,50), OPTR( 20,50),
1 OPTG( 20,50), OPET( 50), OPT ( 50),
2 OPR ( 50), OPGZ( 50)
C
KR = KR + 1
C
OPTE( KT,KR ) = TVR
OPTT( KT,KR ) = T
OPTR( KT,KR ) = RX
OPTG( KT,KR ) = RMA
C
RETURN
END
C **********************************************************************
C
SUBROUTINE COPT2 ( KT, KR )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /OPM/ OPTE( 20,50), OPTT( 20,50), OPTR( 20,50),
1 OPTG( 20,50), OPET( 50), OPT ( 50),
2 OPR ( 50), OPGZ( 50)
C
C ------------------------------------------------
DO 100 J = 1, KR
C
L = 0
C
C --------------------------------------
DO 200 I = 1, KT
C
IF ( OPTG( I,J).LE.1.D0) L = 1
200 CONTINUE
C --------------------------------------
DO 300 I = 1, KT
C
IF ( L.EQ.0) GO TO 10
IF ( OPTG( I,J).LE.1.D0 .AND. OPTE( I,J).LT.OPET( J)) GO TO 20
GO TO 300
C
10 IF ( OPTG( I,J).LT.OPGZ( J)) GO TO 20
GO TO 300
C
20 CONTINUE
OPET( J) = OPTE( I,J)
OPR ( J) = OPTR( I,J)
OPT ( J) = OPTT( I,J)
OPGZ( J) = OPTG( I,J)
C
300 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CRNS ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ******************************
C
F1( X) = 1.D0 - DCOS( X)
F3( X) = DSIN( X)
C
C ******************************
C
S1 = 2.D0
S2 = - 2.D0* ( RX* F1( VAX) + RY* F1( VAY))
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
S6 = - ( BX* F3( VAX) + BY* F3( VAY))
C
S7 = 2.D0
S8 = 2.D0* ( RX* F1( VAX) + RY* F1( VAY))
S9 = 0.D0
S10= 0.D0
S11= 0.D0
S12= BX* F3( VAX) + BY* F3( VAY)
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EFU3 ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
C
C ******************************
C
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C
C ******************************
C
RXF = 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
RYF = 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
C
G = 1.D00
C
IF (( RX+RY).GE.0.001) THEN
F = 1.D0 + ( BX+BY)* ( BX+BY)/( 2.D0* ( RX+RY))
RXF = RX* F
RYF = RY* F
C
G = 1.D00
END IF
C
S1 = 1.D0
S2 = - 2.D0* ( RXF* F1( VAX) + RYF* F1( VAY))
C
S3 = G* ( BX* ( 4.D0* F2( VAX) - F2( 2.D0* VAX) - 3.D0)
1 + BY* ( - F2( 2.D0* VAY) - 3.D0 ))/ 6.D0
C
S4 = 0.D0
S5 = 0.D0

S6 = G* ( BX* ( F3( 2.D0* VAX) - 8.D0* F3( VAX))
1 + BY* ( F3( 2.D0* VAY) - 8.D0* F3( VAY)))/ 6.D0
C
P1 = S1 + S2 + S3
P2 = S6
P3 = 1.D0
P4 = 0.D0
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXCHH ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ******************************

F1( X) = 1.D0 - DCOS( X)
F2( X) = 2.D0 + DCOS( X)
F3( X) = DSIN( X)
C
C ******************************
C
RXF = 0.5D0* (( BX+BY)* ( BX+BY) - 1.D0)/ 3.D0
RYF = 0.5D0* (( BX+BY)* ( BX+BY) - 1.D0)/ 3.D0
G = 1.D00
C
IF (( RX+RY).GE.0.001) THEN
F = 1.D0 + 0.5D0 * ( BX+BY)* ( BX+BY)/( RX+RY)
1 - 1.D0/( 3.D0* ( RX+RY))
RXF = RX* F
RYF = RY* F
C
G = 1.D00
END IF
C
S1 = F2( VAX)* F2( VAY)
S2 = - 6.D0* ( RXF* F1( VAX)* F2( VAY) + RYF* F2( VAX)* F1( VAY))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
S6 = - 3.D0* G* ( BX* F3( VAX)* F2( VAY) + BY* F2( VAX)* F3( VAY))
C
S7 = 27.D0
S8 = 0.D0
S9 = 0.D0
S10= 0.D0
S11= 0.D0
S12= 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = 0.D0
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE EXFU ( P1, P2, P3, P4, RX, RY,BX, BY, VAX, VAY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
C
RXF = 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0 - 0.5D0* ( BX+BY)/ 3.D0
RYF = 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0 - 0.5D0* ( BX+BY)/ 3.D0
C
G = 1.D00
C
IF (( RX+RY).GE.0.001) THEN
F = 1.D0 + (( BX+BY)* ( BX+BY))/( 2.D0* ( RX+RY))
RXF = RX* F
RYF = RY* F
C
G = 1.D0
END IF
C
S1 = 1.D0
C
S2 = - 2.D0* ( RXF* ( 1.D0 - DCOS( VAX))
1 + RYF* ( 1.D0 - DCOS( VAY)))
2 - G* ( BX * ( 1.D0 - DCOS( VAX))
3 + BY * ( 1.D0 - DCOS( VAY)))
C
S3 = - G* ( BX* DSIN( VAX) + BY* DSIN( VAY))
C
S4 = 1.D0
S5 = 0.D0
S6 = 0.D0
C
P1 = S1 + S2
P2 = S3
P3 = S4 + S5
P4 = S6
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE EXHH ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C *******************************
C
F1( X) = 1.D0 - DCOS( X)
F2( X) = 2.D0 + DCOS( X)
F3( X) = DSIN( X)
C
C *******************************
C
S1 = F2( VAX)* F2( VAY)
S2 = - 6.D0* ( RX* F1( VAX)* F2( VAY) + RY* F2( VAX)* F1( VAY))
C
S3 = 0.D0
C
S4 = 0.D0
S5 = 0.D0
S6 = - 3.D0* ( BX* F3( VAX)* F2( VAY)+ BY* F2( VAX)* F3( VAY))
C
S7 = 27.D0
S8 = 0.D0
S9 = 0.D0
C
S10= 0.D0
S11= 0.D0
S12= 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE FTCS ( P1, P2, P3, P4, RX, RY,BX, BY, VAX, VAY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
C
RXF = 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
RYF = 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
C
G = 1.D00
C
IF (( RX+RY).GE.0.001) THEN
F = 1.D0 + (( BX+BY)* ( BX+BY))/( 2.D0* ( RX+RY))
RXF = RX* F
RYF = RY* F
C
G = 1.D00
END IF
C
RXF = RX
RYF = RY
C
G = 1.D00
C
S1 = 1.D0
S2 = - 2.D0*( RXF*( 1.D0 -DCOS( VAX))+ RYF*( 1.D0 - DCOS( VAY)))
C
S3 = - G* ( BX* DSIN( VAX) + BY* DSIN( VAY))
S4 = 1.D0
S5 = 0.D0
S6 = 0.D0
C
P1 = S1 + S2
P2 = S3
P3 = S4 + S5
P4 = S6
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FULY ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
C
RXF = - 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
RYF = - 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
C
G = 1.D00
C
IF (( RX+RY).GE.0.001) THEN
C
F = 1.D0 - (( BX+BY)* ( BX+BY))/( 2.D0* ( RX+RY))
RXF = RX* F
RYF = RY* F
C
G = 1.D00
C
END IF
C
S1 = 1.D0
S2 = 0.D0
S3 = 0.D0
S4 = 1.D0
S5 = 2.D0*( RXF*( 1.D0 -DCOS( VAX)) + RYF*( 1.D0 -DCOS( VAY)))
C
S6 = G* ( BX* DSIN( VAX) + BY* DSIN( VAY))
C
P1 = S1 + S2
P2 = S3
P3 = S4 + S5
P4 = S6
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE GHH ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY,
1 T )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C *******************************
C
F1( X) = 1.D0 - DCOS( X)
F2( X) = 2.D0 + DCOS( X)
F3( X) = DSIN( X)
C
C *******************************
C
S1 = F2( VAX)* F2( VAY)
S2 = - 6.D0* ( 1.D0 - T )* ( RX* F1( VAX)* F2( VAY)
1 + RY* F2( VAX)* F1( VAY))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
S6 = - 3.D0* ( 1.D0 - T )* ( BX* F3( VAX)* F2( VAY)
1 + BY* F2( VAX)* F3( VAY))
C
S7 = S1
S8 = 6.D0* T* ( RX* F1( VAX)* F2( VAY) + RY* F2( VAX)* F1( VAY))
C
S9 = 0.D0
S10= 0.D0
S11= 0.D0
C
S12= 3.D0* T* ( BX* F3( VAX)* F2( VAY) + BY* F2( VAX)* F3( VAY))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE IMFU ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
C
C *******************************
C
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C
C *******************************
C
RXF = - 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0 - 0.5D0* ( BX+BY)/ 3.D0
C
RYF = - 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0 - 0.5D0* ( BX+BY)/ 3.D0
C
G = 1.D00
C
IF (( RX+RY).GE.0.001) THEN
F = 1.D0 - 0.5D0* ( BX+BY)* ( BX+BY)/( RX+RY)
1 - 0.5D0* ( BX+BY)/( RX+RY)
RXF = RX* F
RYF = RY* F
C
G = 1.D00
END IF
C
S1 = 1.D0
S2 = 0.D0
S3 = 0.D0
S4 = 0.D0
C
S5 = 0.D0
S6 = 0.D0
S7 = 1.D0
C
S8 = 2.D0* ( RXF* F1( VAX) + RYF* F1( VAY))
S9 = G* ( BX* F1( VAX) + BY* F1( VAY))
S10= 0.D0
S11= 0.D0
S12= G* ( BX* F3( VAX) + BY* F3( VAY))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE INIT
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /OPM/ OPTE( 20,50), OPTT( 20,50), OPTR( 20,50),
1 OPTG( 20,50), OPET( 50), OPT ( 50),
2 OPR ( 50), OPGZ( 50)
C
DO 100 I = 1, 50
C
OPET( I) = 100.D0
OPGZ( I) = 100.D0
C
100 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE KSLA ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
C
C ********************************
C
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C
C ********************************
C
RXF = - 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
RYF = - 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
C
G = 1.D00
C
IF (( RX+RY).GE.0.001) THEN
F = 1.D0 - ( BX+BY)* ( BX+BY)/( 2.D0* ( RX+RY))
RXF = RX* F
RYF = RY* F
C
G = 1.D00
END IF
C
S1 = 1.D0
S2 = 0.D0
S3 = G* ( BX* F1( VAX) + BY* F1( VAY))
C
S4 = 0.D0
S5 = 0.D0
S6 = 0.D0
C
S7 = 1.D0
S8 = 2.D0* ( RXF* F1( VAX) + RYF* F1( VAY))
S9 = G* ( BX * F1( VAX) + BY * F1( VAY))
C
S10 = 0.D0
S11 = 0.D0
S12 = G* ( BX* F3( VAX) + BY* F3( VAY))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE LXWD ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
C
C *******************************
C
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C
C *******************************
C
RXF = 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
1 - 0.5D0* ( BX* BX + BY* BY)/ 3.D0
C
RYF = 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
1 - 0.5D0* ( BX* BX + BY* BY)/ 3.D0
C
G = 1.D00
C
IF (( RX+RY).GE.0.001) THEN
F = 1.D0 + ( BX+BY)* ( BX+BY)/( 2.D0* ( RX+RY))
RXF = RX* F
RYF = RY* F
C
G = 1.D00
END IF
C
S1 = 1.D0
S2 = - 2.D0* ( RXF* F1( VAX) + RYF* F1( VAY))
1 - G* ( BX* BX* F1( VAX) + BY* BY* F1( VAY))
S3 = 0.D0
C
S4 = 0.D0
S5 = 0.D0
S6 = - G* ( BX* F3( VAX) + BY* F3( VAY))
C
S7 = 1.D0
S8 = 0.D0
S9 = 0.D0
C
S10= 0.D0
S11= 0.D0
S12= 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MGHH ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY, T )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C *******************************
C
F1( X) = 1.D0 - DCOS( X)
F2( X) = 2.D0 + DCOS( X)
F3( X) = DSIN( X)
C
C *******************************
C
RXF = 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
RYF = 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
RZF = 0.5D0* ( BX+BY)* ( BX+BY)/ 3.D0
C
G = 1.D0
C
IF (( RX+RY).GE.1.D-10 ) THEN
F = 1.D0 + 0.5D0* ( BX+BY)* ( BX+BY)/( RX+RY)
RXF = RX* F
RYF = RY* F
C
G = 1.D0
END IF
C
S1 = F2( VAX)* F2( VAY)
C
S2 = - 6.D0* ( 1.D0 - T )* ( RXF* F1( VAX)* F2( VAY)
1 + RYF* F2( VAX)* F1( VAY))
C
S6 = - 3.D0* G* ( BX* F3( VAX)* F2( VAY)+ BY* F2( VAX)* F3( VAY))
C
S7 = S1
S8 = 6.D0* T* ( RXF* F1( VAX)* F2( VAY) + RYF* F2( VAX)* F1( VAY))
C
P1 = S1 + S2
P2 = S6
P3 = S7 + S8
P4 = 0.D0
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE OUTP1 ( RX, RY, RMA, SGA, SGC, TVR )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
WRITE(6,2000) RX, RMA, TVR
WRITE(7,2000) RX, RMA, TVR
2000 FORMAT(' R =',F5.2,3X,'Gzai=',F9.6,3X,'Et.Error=',F9.6)
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE OUTPT ( KR )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
COMMON /NUM/ MTS, NT, NRX, NBX, NKX, NRY, NBY, NKY
COMMON /OPM/ OPTE( 20,50), OPTT( 20,50), OPTR( 20,50),
1 OPTG( 20,50), OPET( 50), OPT ( 50),
2 OPR ( 50), OPGZ( 50)
C
WRITE(6,2000) ( TV( I), I = 1, NT )
WRITE(7,2000) ( TV( I), I = 1, NT )
2000 FORMAT(/,' * Check Time Sch. Par.=',11F5.1/24X,11F5.1)
C
WRITE(6,2100)
WRITE(7,2100)
2100 FORMAT(/' ** Optimum Time Scheme Par. **')
C
DO 100 I = 1, KR
C
WRITE(6,2200) OPR( I), OPT( I), OPGZ( I), OPET( I)
WRITE(7,2200) OPR( I), OPT( I), OPGZ( I), OPET( I)
2200 FORMAT(3X,'R =',F5.2,3X,'Opt.Time Sc.Par.=',F5.2,
1 3X,'Gzai Max.=',F9.6,3X,'ET=',F9.6)
C
100 CONTINUE
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE PCON
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CST/ PAI
COMMON /MTV/ MT( 50), TV( 14)
COMMON /NUM/ MTS, NT, NRX, NBX, NKX, NRY, NBY, NKY
COMMON /ABC/ RVX(10), RVY(10), BVX(10), BVY(10), KMX(14), KMY(14)
C
PAI = 3.14159 26535 89793 D0
C
WRITE(7,2000) MTS,NT,NRX,NRY,NBX,NBY,NKX,NKY
WRITE(6,2000) MTS,NT,NRX,NRY,NBX,NBY,NKX,NKY
2000 FORMAT(/ 2X,' * No. of Methods = ',I3,5X,' * No. of Theta =',I3,
1 // 5X,' * No. of Rx = ',I3,' N_Ry=',I3, /
2 5X,' * No. of Bx = ',I3,' N_By=',I3 /
3 5X,' * No. of Kx = ',I3,' N_Ky=',I3 )
C
WRITE(6,2100) ( TV( I),I = 1, NT )
WRITE(7,2100) ( TV( I),I = 1, NT )
2100 FORMAT(/1X,'THETA =',14F5.1,/,8X,14F5.1)
C
WRITE(6,2200) ( RVX( I),I = 1, NRX)
WRITE(6,2300) ( RVY( I),I = 1, NRY)
WRITE(7,2200) ( RVX( I),I = 1, NRX)
WRITE(7,2300) ( RVY( I),I = 1, NRY)
C
2200 FORMAT(3X,' RX = ',15(F6.3,2X))
2300 FORMAT(3X,' RY = ',15(F6.3,2X))
C
WRITE(6,2400) ( BVX( I),I = 1, NBX)
WRITE(6,2500) ( BVY( I),I = 1, NBY)
WRITE(7,2400) ( BVX( I),I = 1, NBX)
WRITE(7,2500) ( BVY( I),I = 1, NBY)
C
2400 FORMAT(/3X,' BX = ',15(F6.3,2X))
2500 FORMAT( 3X,' BY = ',15(F6.3,2X))
C
WRITE(6,2600) ( KMX( I),I = 1, NKX)
WRITE(6,2700) ( KMY( I),I = 1, NKY)
WRITE(7,2600) ( KMX( I),I = 1, NKX)
WRITE(7,2700) ( KMY( I),I = 1, NKY)
C
2600 FORMAT(/2X,' KX = ',16(I3,2X),/,16X,10(I3,2X))
2700 FORMAT( 2X,' KY = ',16(I3,2X),/,16X,10(I3,2X))
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE PTTL ( NMH )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
C
DIMENSION NMT( 12)
CHARACTER*9 NMT
C
C ----- FEMs & FDMs ---
C
DATA NMT/ 'MGHH', 'EXCHH', 'LXWD ', 'SZE ', 'EXFU ',
1 'EFU3', 'FTCS', 'FULY', 'CRNS', 'RBWS ',
2 'KSLA ', 'IMFU ' /
C
WRITE(6,2000) NMH, NMT( NMH )
WRITE(7,2000) NMH, NMT( NMH )
C
2000 FORMAT(/1X,'METHOD = ',I3,2X,'* ',A9,' *')
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RBWS ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C *******************************
C
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C
C *******************************
C
S1 = 2.D0
S2 = - 2.D0* ( RX* F1( VAX) + RY* F1( VAY))
S3 = ( BX* F1( VAX) + BY* F1( VAY))
S4 = 0.D0
S5 = 0.D0
S6 = - ( BX* F3( VAX) + BY* F3( VAY))
C
S7 = 2.D0
S8 = 2.D0* ( RX* F1( VAX) + RY* F1( VAY))
S9 = BX* F1( VAX) + BY* F1( VAY)
S10= 0.D0
S11= 0.D0
S12= BX* F3( VAX) + BY* F3( VAY)
C
P1 = S1 + S2 + S3
P2 = S6
P3 = S7 + S8 + S9
P4 = S12
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE SMTD ( NMH, P1, P2, P3, P4, RX, RY, BX, BY,
1 VAX, VAY, T )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
C
GO TO ( 1, 2, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20 ), NMH
C
C ****** FEMs ***************************************************
C
C
C ----- MGM ---
C
1 CALL MGHH ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY, T )
C
GO TO 99
C
C ----- Explicit ---
C
2 CALL EXCHH ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
GO TO 99
C
C ****** FDMs ***************************************************
C
C ----- Explicit FDMs ---
C
11 CALL LXWD ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
GO TO 99
C
12 CALL SZE ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
GO TO 99
C
13 CALL EXFU ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
GO TO 99
C
14 CALL EFU3 ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
GO TO 99
C
15 CALL FTCS ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
GO TO 99
C
C ----- Implicit FDMs ---
C
16 CALL FULY ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
GO TO 99
C
17 CALL CRNS ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
GO TO 99
C
18 CALL RBWS ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
GO TO 99
C
19 CALL KSLA ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
GO TO 99
C
20 CALL IMFU ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
99 RETURN
END
C ***********************************************************************
C
SUBROUTINE SZE ( P1, P2, P3, P4, RX, RY, BX, BY, VAX, VAY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /MTV/ MT(50), TV(14)
C
C ******************************
C
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C
C ******************************
C
RXF = ( BX* BY)/ 3.D0
RYF = ( BX* BY)/ 3.D0
RZF = ( BX* BY)/ 3.D0
C
G = 1.D0
C
IF (( RX+RY).GE.0.001D0 ) THEN
F = 1.D0 + ( BX* BY)/( RX + RY)
RXF = RX* F
RYF = RY* F
C
G = 1.D0
END IF
C
S1 = 1.D0
S2 = - 2.D0* ( RXF* F1(VAX) + RYF* F1(VAY))
C
S3 = G* ( - ( BX* BX + BY* BY)/ 4.D0 - 3.D0* ( BX + BY)/ 4.D0
1 + BX* F2( VAX) + BY* F2( VAY)
2 + ( BX* BX - BX)* F2( 2.D0* VAX)/ 4.D0
3 + ( BY* BY - BY)* F2( 2.D0* VAY)/ 4.D0 )
C
S4 = 0.D0
S5 = 0.D0
C
S6 = 1.D0* G* ( ( BX* BX - 3.D0* BX)* F3( VAX)/ 2.D0
1 + ( BY* BY - 3.D0* BY)* F3( VAY)/ 2.D0
2 - ( BX* BX - BX)* F3( 2.D0* VAX)/ 4.D0
3 - ( BY* BY - BY)* F3( 2.D0* VAY)/ 4.D0 )
C
S7 = 1.D0
S8 = 0.D0
S9 = 0.D0
C
S10= 0.D0
S11= 0.D0
S12= 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C ************************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,*) ' '
WRITE(6,*) '*****************************************************'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Error Analysis Software for C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C 2-Dimensional Convection- Diffusion Equation C'
WRITE(6,*) 'C === C'
WRITE(6,*) 'C ( V.4.0 ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Copyright 2013 : Yasuhiro MATSUDA C'
WRITE(6,*) 'C C'
WRITE(6,*) '*****************************************************'
C
RETURN
END
C ***********************************************************************