戻る

C
PROGRAM MAIN
C ******************************************************************** C
C C
C Error Analysis Software C
C C
C for 3D - Diffusion Equation C
C C
C ( For FEMs ) C
C C
C ( V.3.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, NRY, NRZ, NKX, NKY, NKZ
COMMON /ABC/ RVX( 10), RVY( 10), RVZ( 10), KMX( 14), KMY( 14),
1 KMZ( 14)
COMMON /NAMEM/ NMT( 2)
CHARACTER*6 NMT
C
DIMENSION OPTE( 20, 50), OPTT( 20, 50), OPTR( 20, 50),
1 OPTG( 20, 50), OPET( 50), OPT ( 50),
2 OPR ( 50), OPGZ( 50)
C
C --------------------------------------------
C
OPEN ( 7, FILE='DIF3D_FEM_EA.RESULT' )
C
C --------------------------------------------
C
CALL DINPT
C
C -------------------------------------------------------------------
DO 100 I = 1, MTS
C
WRITE(6,2000) NMT( I)
WRITE(7,2000) NMT( I)
2000 FORMAT(/,1X,' *** ',A6,' ***',/ )
C
KT = 0
NTT = NT
IF ( MT(I).EQ.2) NTT = 1
C
C ---------------------------------------------------------
DO 200 J = 1, NTT
C
T = TV( J)
KT = KT + 1
IF ( MT(I).EQ.2) GO TO 10
WRITE(7,2100) T
2100 FORMAT(/,3X,'** Time Scheme Par. = ',F6.3,' **'/)
C
10 CONTINUE
KR = 0
NKXYZ = NKX* NKY* NKZ
C
CALL CLEAR ( TSGA, TSGC, DAB3, DUMY)
C
C ------------------------------------------------
DO 300 IR = 1, NRX
C
RX = RVX( IR)
RY = RVY( IR)
RZ = RVZ( IR)
C
CALL CLEAR ( SGA, SGC, DAB2, RMX )
C
CALL INIT ( OPET, OPGZ)
C
CALL CLEAR ( BSA, BSC, DAB1, RMXB )
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
C ----------------------------
DO 500 IKZ = NKZ, 1, -1
C
KZ = KMZ( IKZ)
VAZ = 2.D0* PAI/ DFLOAT( KZ)
C
C ------ Schemes ---
C
CALL SMTHD ( P1, P2, P3, P4, RX, RY, RZ, VAX, VAY, VAZ, I, T )
C
C -------------------
C
CALL CMP1 ( P1, P2, P3, P4, RX, RY, RZ, VAX, VAY,
1 VAZ, A, C, DAB0, BA, BC, TS, DAB1, BSA,
2 BSC, RMXB, SGA, SGC, RMX )
C
500 CONTINUE
C -----------------------------
C
CALL CMP2 ( NKXYZ, DAB1, DAB2, BSA, BSC, TBS )
C
CALL CMP3 ( NKXYZ, DAB2, DAB3, SGA, SGC, TSGA, TSGC, TVR )
C
C ----------------------------
C
IF ( MT(I).EQ.1) CALL COPT1 ( KT, KR, T, RX, RMX, TVR,
1 OPTE, OPTT, OPTR, OPTG )
C
IF ( MT(I).EQ.2) CALL OUTP1 ( RX, RY, RZ, RMX, SGA, SGC,
1 TVR, DAB2 )
C
300 CONTINUE
C -----------------------------------------------
C
CALL CALC4 ( NRX, NKXYZ, DAB3, TSGA, TSGC, SUMT )
C
200 CONTINUE
C ----------------------------------------------------------
C
IF ( MT(I).EQ.1 ) THEN
C
CALL COPT2 ( KT, KR, OPTE, OPTT, OPTR, OPTG, OPET,
1 OPR, OPT, OPGZ)
C
CALL OUTP2 ( KR, OPET, OPR, OPT, OPGZ) ! Opt. THEtA
C
END IF
C -------------
C
100 CONTINUE
C --------------------------------------------------------------------
WRITE(6,2200)
WRITE(7,2200)
2200 FORMAT(/,5X,'** END **')
C
CLOSE( 7)
C
STOP
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 3-Dimensional Convection- Diffusion Equation C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( for Finite Element MEthods ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( V.3.0 ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Copyright 2011 : Yasuhiro MATSUDA C'
WRITE(6,*) 'C C'
WRITE(6,*) '*****************************************************'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMP1 ( P1, P2, P3, P4, RX, RY, RZ, VAX,
1 VAY, VAZ, AA, CC, DAB0, BA, BC, TS,
2 DAB1, BSA, BSC, RMXB, SGA, SGC, RMX )
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 )
C
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)
C
IF ( DABS( CC).GT.PAI) CC = 2.D0* PAI - DABS( CC)
CC = CC + 1.D0
C
VAL = 0.D0
C
AP = DEXP( - VAX* VAX* RX - VAY* VAY* RY - VAZ* VAZ* RZ)
C
AE = AP - AE
BE = - BE
C
DAB1 = DAB1 + AE* AE + BE* BE
DAB0 = DSQRT( AE* AE + BE* BE )
C
IF ( AA.GE.RMXB ) RMXB = AA
BSA = BSA + ( AA - AP) * ( AA - AP)
BSC = BSC + ( CC - 1.D0)* ( CC - 1.D0)
C
BA = DSQRT( ( AA - AP) * ( AA - AP ))
BC = DSQRT( ( CC - 1.D0)* ( CC - 1.D0))
TS = DSQRT( BA* BA + BC* BC )
C
IF ( AA.GE.RMX) RMX = AA
C
SGA = SGA + ( AA - AP )* ( AA - AP )
SGC = SGC + ( CC - 1.D0)* ( CC - 1.D0)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMP2 ( NKXYZ, DAB1, DAB2, BSA, BSC, TBS )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DAB2 = DAB2 + DAB1
C
DAB1 = DSQRT( DAB1/ DFLOAT( NKXYZ))
BSA = DSQRT( BSA / DFLOAT( NKXYZ))
BSC = DSQRT( BSC / DFLOAT( NKXYZ))
C
TBS = DSQRT( BSA* BSA + BSC * BSC )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMP3 ( NKXYZ, DAB2, DAB3, SGA, SGC, TSGA, TSGC, TVR )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DAB3 = DAB2 + DAB3
DAB2 = DSQRT( DAB2/( DFLOAT( NKXYZ)))
TSGA = TSGA + SGA
C
SGA = DSQRT( SGA/( DFLOAT( NKXYZ)))
TSGC = TSGC + SGC
C
SGC = DSQRT( SGC/( DFLOAT( NKXYZ)))
C
TVR = DSQRT( SGA* SGA + SGC* SGC )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CALC4 ( NRX, NKXYZ, DAB3, TSGA, TSGC, SUMT )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DAB3 = DSQRT( DAB3/( DFLOAT( NRX)* DFLOAT( NKXYZ)))
C
TSGA = DSQRT( TSGA/( DFLOAT( NRX)* DFLOAT( NKXYZ)))
TSGC = DSQRT( TSGC/( DFLOAT( NRX)* DFLOAT( NKXYZ)))
C
SUMT = DSQRT( TSGA* TSGA + TSGC* TSGC )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CLEAR ( A, C, D, RM )
C
IMPLICIT REAL*8(A-H,O-Z)
C
A = 0.D0
C = 0.D0
D = 0.D0
RM = 0.D0
C
RETURN
END
C **********************************************************************
C
SUBROUTINE COPT1 ( KT, KR, T, RX, RMX, TVR, OPTE, OPTT, OPTR,
1 OPTG )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION OPTE( 20,50), OPTT( 20,50), OPTR( 20,50), OPTG( 20,50)
C
KR = KR + 1
C
OPTE( KT,KR) = TVR
OPTT( KT,KR) = T
OPTR( KT,KR) = RX
OPTG( KT,KR) = RMX
C
RETURN
END
C **********************************************************************
C
SUBROUTINE COPT2 ( KT, KR, OPTE, OPTT, OPTR, OPTG, OPET,
1 OPR, OPT, OPGZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION OPTE( 20,50), OPTT( 20,50), OPTR( 20,50),
1 OPTG( 20,50), OPET( 50), OPR ( 50),
2 OPT ( 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 ----------------------------------------------------------
RETURN
END
C **********************************************************************
C
SUBROUTINE FEM_EXP ( P1, P2, P3, P4, RX, RY, RZ, VAX, VAY, VAZ )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- EXCHH ---
C
C ---------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = 2.D0 + DCOS( X)
F3( X) = DSIN( X)
C ---------------------------------
C
RXF = - 0.5D0 / 3.D0
RYF = - 0.5D0 / 3.D0
RZF = - 0.5D0 / 3.D0
G = 1.D00
C
IF ( ( RX+RY+RZ).GE.0.001) THEN
C
F = 1.D0 - 1.D0/( 2.D0* ( RX + RY + RZ))
C
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
C
END IF
C
S1 = F2( VAX)* F2( VAY)* F2( VAZ)
S2 = - 6.D0* ( RXF* F1( VAX)* F2( VAY)* F2( VAZ)
1 + RYF* F2( VAX)* F1( VAY)* F2( VAZ)
2 + RZF* F2( VAX)* F2( VAY)* F1( VAZ))
C
P1 = S1 + S2
P2 = 0.D0
P3 = 27.D0
P4 = 0.D0
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INIT ( OPET, OPGZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION OPET( 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 FEM_IMP ( P1, P2, P3, P4, RX, RY, RZ, VAX, VAY, VAZ,
1 T )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- MGHH : GHH : GCHH ---
C
C ---------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = 2.D0 + DCOS( X)
F3( X) = DSIN( X)
C ---------------------------------
C
RXF = 0.D0
RYF = 0.D0
RZF = 0.D0
C
G = 1.D0
C
IF ( ( RX + RY + RZ).GE.1.D-10) THEN
F = 1.D0
RXF = RX* F
RYF = RY* F
RZF = RZ* F
C
G = 1.D0
END IF
C
S1 = F2( VAX)* F2( VAY)* F2( VAZ)
C
S2 = - 6.D0* ( 1.D0 - T)
1 * ( RXF* F1( VAX)* F2( VAY)* F2( VAZ)
2 + RYF* F2( VAX)* F1( VAY)* F2( VAZ)
3 + RZF* F2( VAX)* F2( VAY)* F1( VAZ))
C
S7 = S1
S8 = 6.D0* T* ( RXF* F1( VAX)* F2( VAY)* F2( VAZ)
1 + RYF* F2( VAX)* F1( VAY)* F2( VAZ)
2 + RZF* F2( VAX)* F2( VAY)* F1( VAZ))
C
P1 = S1 + S2
P2 = 0.D0
P3 = S7 + S8
P4 = 0.D0
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PCON
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /NAMEM/ NMT( 2)
CHARACTER*6 NMT
COMMON /MTV/ MT( 50), TV( 14)
COMMON /NUM/ MTS, NT, NRX, NRY, NRZ, NKX, NKY, NKZ
COMMON /ABC/ RVX( 10), RVY( 10), RVZ( 10), KMX( 14), KMY( 14),
1 KMZ( 14)
C
WRITE(7,2000) MTS, NMT(1), NMT(2), NT, NRX, NRY, NRZ, NKX, NKY,
1 NKZ
WRITE(6,2000) MTS, NMT(1), NMT(2), NT, NRX, NRY, NRZ, NKX, NKY,
1 NKZ
2000 FORMAT(/3X,' * No. of MEthods = ',I3,' ** ',A6,3X,' ** ',A6,
1 3X,//,5X,' * No. of ThEta =',I3,/,5X,
2 ' * No. of Rx = ',I3,' N_Ry=',I3,' N_Rz=',I3/5X,
3 ' * No. of Kx = ',I3,' N_Ky=',I3,' N_Kz=',I3/ )
C
WRITE(7,2100) ( MT(I),I = 1, MTS )
2100 FORMAT( 3X,' * MEtHOD = ',20I3 )
C
WRITE(6,2200) ( TV( I),I = 1, NT )
WRITE(7,2200) ( TV( I),I = 1, NT )
2200 FORMAT(1X,'THEtA =',11F6.1,/,8X,9F6.1,/,9X,9F6.1,/,9X,9F6.1,/,
1 9X,9F6.1 )
C
WRITE(6,2300) ( RVX( I),I = 1,NRX )
WRITE(6,2400) ( RVY( I),I = 1,NRY)
WRITE(6,2500) ( RVZ( I),I = 1,NRZ)
C
WRITE(7,2300) ( RVX( I),I = 1,NRX )
WRITE(7,2400) ( RVY( I),I = 1,NRY)
WRITE(7,2500) ( RVZ( I),I = 1,NRZ)
2300 FORMAT(/3X,' RX = ',15(F6.3,2X))
2400 FORMAT( 3X,' RY = ',15(F6.3,2X))
2500 FORMAT( 3X,' RZ = ',15(F6.3,2X))
C
WRITE(6,2600) ( KMX( I),I = 1,NKX )
WRITE(6,2700) ( KMY( I),I = 1,NKY)
WRITE(6,2800) ( KMZ( I),I = 1,NKZ)
C
WRITE(7,2600) ( KMX( I),I = 1,NKX )
WRITE(7,2700) ( KMY( I),I = 1,NKY)
WRITE(7,2800) ( KMZ( I),I = 1,NKZ)
2600 FORMAT(/2X,' KX = ',16(I3,2X),/,16X,10(I3,2X))
2700 FORMAT( 2X,' KY = ',16(I3,2X),/,16X,10(I3,2X))
2800 FORMAT( 2X,' KZ = ',16(I3,2X),/,16X,10(I3,2X)// )
C
RETURN
END
C **********************************************************************
C
BLOCK DATA
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON / NAMEM / NMT( 2)
CHARACTER*6 NMT
COMMON /MTV/ MT( 50), TV( 14)
COMMON /NUM/ MTS, NT, NRX, NRY, NRZ, NKX, NKY, NKZ
COMMON /ABC/ RVX( 10), RVY( 10), RVZ( 10), KMX( 14), KMY( 14),
1 KMZ( 14)
C
DATA NMT / 'FEM-Im ','FEM-Ex '/
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, NKZ / 14, 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/
DATA KMZ / -3, -4, -6, -8, - 10, - 20, -40, 3, 4, 6, 8, 10,20,40/
C
DATA MTS, NT / 2, 14/ ! MEthod No. / No. of THEtA
C
DATA NRX, NRY, NRZ / 6, 6, 6/
DATA RVX / 0.01, 0.02, 0.04, 0.1, 0.2, 0.3, 0., 0., 0., 0./
DATA RVY / 0.01, 0.02, 0.04, 0.1, 0.2, 0.3, 0., 0., 0., 0./
DATA RVZ / 0.01, 0.02, 0.04, 0.1, 0.2, 0.3, 0., 0., 0., 0./
C
END
C ***********************************************************************
C
SUBROUTINE DINPT
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CST/ PAI
COMMON / NAMEM / NMT( 2)
CHARACTER*6 NMT
COMMON /MTV/ MT( 50), TV( 14)
COMMON /NUM/ MTS, NT, NRX, NRY, NRZ, NKX, NKY, NKZ
C
C ----------------
C
CALL TTLE
C
C ----------------
C
PAI = 3.14159 26535 89793 D0
C
MT( 1) = 1
MT( 2) = 2
C
C ----------------
C
CALL PCON
C
C ----------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SMTHD ( P1, P2, P3, P4, RX, RY, RZ, VAX, VAY, VAZ,
1 I, T )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /MTV/ MT( 50), TV( 14)
C
GO TO ( 1, 2) , MT( I)
C
C ***** FEMs ***
C
C ----- Implicit MEthod ---
C
1 CALL FEM_IMP ( P1, P2, P3, P4, RX, RY, RZ, VAX, VAY, VAZ, T )
C
GO TO 99
C
C ----- Explicit MEthd ---
C
2 CALL FEM_EXP ( P1, P2, P3, P4, RX, RY, RZ, VAX, VAY, VAZ)
C
99 RETURN
END
C **********************************************************************
C
SUBROUTINE OUTP1 ( RX, RY, RZ, RMX, SGA, SGC, TVR, DAB2 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /DEB/ IDEB( 10)
C
WRITE(6,2000) RX, RMX, TVR
WRITE(7,2000) RX, RMX, TVR
2000 FORMAT(' R =',F5.2,3X,'Gzai =',F8.5,3X,'Et.Err =',F8.5)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE OUTP2 ( KR, OPET, OPR, OPT, OPGZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION OPET( 50), OPT( 50), OPR( 50), OPGZ( 50)
C
WRITE(7,2000)
2000 FORMAT(' *** Optimum Time Scheme Par.*** ')
C
DO 100 I = 1, KR
C
WRITE(6,2100) OPR( I), OPT( I), OPGZ( I), OPET( I)
WRITE(7,2100) OPR( I), OPT( I), OPGZ( I), OPET( I)
2100 FORMAT(3X,'R =',F5.2,3X,'OPT. Time Sc. Par.=',F5.2,
1 3X,'Gzai Max.=',F8.5,3X,'Et =',F8.5)
C
100 CONTINUE
C
RETURN
END
C **********************************************************************