戻る


PROGRAM MAIN
C ******************************************************************** C
C C
C DIF3D - FDM C
C C
C ( Explicit Method : FTCS ) 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
DIMENSION X ( 9), Y( 4), Z ( 4),
1 CO( 9,4,4), C( 9,4,4), CAA( 9,4,4)
C
COMMON /CST3/ NPR,IP
COMMON /RBVL/ RX, RY, RZ
COMMON /CST1/ RKX, RKY, RKZ
COMMON /CST2/ DT, LMT, LOP, ER
C
OPEN ( 5, FILE='FTCS.DAT' )
OPEN ( 7, FILE='OUTFTC.DAT')
C
C ----- Prob. Example ---
C
NX = 9
NY = 4
NZ = 4
C
C -----------------------
C
CALL TTLE
C
C ----------------------------------------------------
C
CALL INIT ( X, Y, Z, C, CO, TIM, NX, NY, NZ )
C
C ----------------------------------------------------
900 LOP = LOP + 1
TIM = DT* DFLOAT( LOP)
C
C -------------------------------------------------
C
CALL FTCS ( CO, C, NX, NY, NZ )
C
C -------------------------------------------------
C
CALL ANSL ( X, Y, Z, CAA, TIM, NX, NY, NZ )
C
C ----------------------------------------------------------
C
CALL ROUT ( C, CAA, CO, X, TIM, NX, NY, NZ, *900 )
C
C ----------------------------------------------------------
C
CLOSE ( 5)
CLOSE ( 7)
C
STOP
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8 (A-H,O-Z)
C
WRITE(6,100)
WRITE(6,101) 'C *********************************************** C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C DIF3D - FDM C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C ( FTCS ) C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C Three-Dimensional Diffusion Analysis by FDM C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C ( V.3.0 ) 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 ANSL ( X, Y, Z, CAA, TIM, NX, NY, NZ )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X( NX), Y( NY), Z( NZ), CAA( NX, NY, NZ)
C
COMMON /CST1/ RKX, RKY, RKZ
C
T = TIM
DKT = 0.0192D0
C
DDX = 0.05D0
DDY = 0.05D0
DDZ = 0.05D0
C
C ----------------------------------------------------------
DO 100 I = 1, NX
DO 100 J = 1, NY
DO 100 K = 1, NZ
C ----------------------------------------------------------
C
CF1 = ( DDX - X( I))/( 2.D0* DSQRT( RKX* T + DKT))
DF1 = ERF( CF1)
C
CF2 = ( DDX + X( I))/( 2.D0* DSQRT( RKX* T + DKT))
DF2 = ERF( CF2)
C
CF3 = ( DDY - Y( J))/( 2.D0* DSQRT( RKY* T + DKT))
DF3 = ERF( CF3)
C
CF4 = ( DDY + Y( J))/( 2.D0* DSQRT( RKY* T + DKT))
DF4 = ERF( CF4)
C
CF5 = ( DDZ - Z( K))/( 2.D0* DSQRT( RKZ* T + DKT))
DF5 = ERF( CF5)
C
CF6 = ( DDZ + Z( K))/( 2.D0* DSQRT( RKZ* T + DKT))
DF6 = ERF( CF6)
C
CAA( I,J,K) = 12.5D0* ( DF1 + DF2)* ( DF3 + DF4)* ( DF5 + DF6)
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
FUNCTION ERF( A)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
L = 100
DERF = 0.D0
DH = A / DFLOAT( L)
C
C ----------------------------------------------------------
DO 100 I = 1, L
C
DU = DH* DFLOAT( I)
DV = DH* DFLOAT( I-1)
C
DERF = DERF + ( DEXP( - DU* DU) + DEXP( - DV* DV))/ 2.D0* DH
C
100 CONTINUE
C ----------------------------------------------------------
C
ERF = DERF* 2.D0/ DSQRT( 3.1415 92653 58979 323 D0)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ROUT ( C, CAA, CO, X, TIM, NX, NY, NZ, * )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION C( NX,NY,NZ), CAA( NX,NY,NZ), CO( NX,NY,NZ), X( NX)
C
COMMON /CST3/ NPR, IP
COMMON /CST2/ DT, LMT, LOP, ER
C
C --------------------------------------
CAMX = CAA( 1,1,1)
C
DO 100 I = 1, NX
DO 100 J = 1, NY
DO 100 K = 1, NZ
C --------------------------------------
C
CAMX = DMAX1( CAA( I,J,K), CAMX )
C
100 CONTINUE
C --------------------------------------
C
C ----- Computational Error ---
C
ER2 = 0.D0
C ------------------------------------------
DO 200 I = 1, NX
DO 200 J = 1, NY
DO 200 K = 1, NZ
C
ER2 = ER2 + DABS( CAA( I,J,K) - C( I,J,K))/ CAMX* 100.D0
C
200 CONTINUE
C ------------------------------------------
ER = ER + ER2/ DFLOAT( NX* NY* NZ )
C
C ------------------------------------------
IF ( MOD( LOP,NPR).EQ.0) THEN
WRITE(7,2000) LOP, TIM, CAMX, ER/ DFLOAT( LOP)
WRITE(6,2000) LOP, TIM, CAMX, ER/ DFLOAT( LOP)
2000 FORMAT(/,' * LOP =',I4,' T =',F6.3,' Can_max.=',F8.3,
1 ' Error =',F8.3,' (%)')
C
C ---------------------------
DO 300 I = ( NX+1)/2, NX
C
WRITE(7,2300) X( I), CAA( I,1,1), C( I,1,1)
WRITE(6,2300) X( I), CAA( I,1,1), C( I,1,1)
2300 FORMAT(5X,' X =',F6.3,2(3X,F11.7))
C
300 CONTINUE
C ----------------------------
END IF
C ------------------------------------------
C
C --------------------------------------
DO 400 I = 1, NX
DO 400 J = 1, NY
DO 400 K = 1, NZ
C --------------------------------------
C
CO( I,J,K) = C( I,J,K)
C
400 CONTINUE
C --------------------------------------
IF ( LOP.LT.LMT) RETURN 1
C
WRITE(7,2100) ER/ DFLOAT( LMT)
WRITE(6,2100) ER/ DFLOAT( LMT)
2100 FORMAT(/,5X,'*** Computational Error =',F11.6,' (%)')
C
WRITE(7,2200)
WRITE(6,2200)
2200 FORMAT(/,' ** Normal End **')
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FTCS ( CO, C, NX, NY, NZ )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION CO( NX, NY, NZ), C( NX, NY, NZ), EC( 7)
C
COMMON /RBVL/ RX, RY, RZ
C
EC( 1) = RX
EC( 2) = RY
EC( 3) = RZ
C
EC( 4) = 1.D0 - 2.D0* ( RX + RY + RZ )
C
EC( 5) = RZ
EC( 6) = RY
EC( 7) = RX
C
C ----------------------------
DO 100 I = 1, NX
DO 100 J = 1, NY
DO 100 K = 1, NZ
C ----------------------------
C
C( I,J,K) = 0.D0
100 CONTINUE
C ----------------------------
DO 200 I = 1, NX
DO 200 J = 1, NY
DO 200 K = 1, NZ
C ----------------------------
C
C ----------------------------
IF ( I.EQ.1) THEN
CX = ( EC( 1) + EC( 7))* CO( I+1,J,K)
C
ELSE IF ( I.EQ.NX) THEN
CX = ( EC( 1) + EC( 7))* CO( I-1,J,K)
ELSE
CX = EC( 1)* CO( I-1,J,K) + EC( 7)* CO( I+1,J,K)
END IF
C ----------------------------
IF (J.EQ.1) THEN
CY = ( EC( 2) + EC( 6))* CO( I,J+1,K)
C
ELSE IF (J.EQ.NY) THEN
CY = ( EC( 2) + EC( 6))* CO( I,J-1,K)
ELSE
CY = EC( 2)* CO( I,J-1,K) + EC( 6)* CO( I,J+1,K)
END IF
C ----------------------------
IF ( K.EQ.1) THEN
CZ = ( EC( 3) + EC( 5))* CO( I,J,K+1)
C
ELSE IF (K.EQ.NZ) THEN
CZ = ( EC( 3) + EC( 5))* CO( I,J,K-1)
ELSE
CZ = EC( 3)* CO( I,J,K-1) + EC( 5)* CO( I,J,K+1)
END IF
C ----------------------------
C
C( I,J,K) = EC( 4)* CO( I,J,K) + CX + CY + CZ
C
200 CONTINUE
C ----------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INIT ( X, Y, Z, C, CO, TIM, NX, NY, NZ )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X( NX), Y( NY), Z( NZ), C( NX,NY,NZ), CO( NX,NY,NZ)
C
COMMON /CST3/ NPR, IP
COMMON /RBVL/ RX, RY, RZ
COMMON /CST1/ RKX, RKY, RKZ
COMMON /CST2/ DT, LMT, LOP, ER
C
LOP = 0
TIM = 0.D0
C
ER = 0.D0
C ----------------------------
DO 100 I = 1, NX
DO 100 J = 1, NY
DO 100 K = 1, NZ
C ----------------------------
C
CO( I,J,K) = 0.D0
C ( I,J,K) = 0.D0
C
100 CONTINUE
C ----------------------------
REWIND 5
C
READ(5,1000) DX, DY, DZ
READ(5,1000) RKX, RKY, RKZ
1000 FORMAT(3F10.5)
READ(5,1100) DT, LMT, NPR, IP
1100 FORMAT(F10.5,3I5)
C
FKX = RKX
FKY = RKY
FKZ = RKZ
C
RX = FKX* DT/ DX/ DX
RY = FKY* DT/ DY/ DY
RZ = FKZ* DT/ DZ/ DZ
C
C --------------------------------------
DO 200 I = 1, NX
C
X( I) = - 1.D0 + DFLOAT( I-1)* DX
200 CONTINUE
C --------------------------------------
DO 300 J = 1, NY
C
Y( J) = DFLOAT( J-1)* DY
300 CONTINUE
C --------------------------------------
DO 400 K = 1, NZ
C
Z( K) = DFLOAT( K-1)* DZ
400 CONTINUE
C --------------------------------------
C
C ----- Initial Condition -------------------------
C
CALL ANSL ( X, Y, Z, CO, TIM, NX, NY, NZ )
C
C -------------------------------------------------
C
WRITE(7,2200) 'FTCS',DX,DY,DZ,RKX,FKX,RKY,FKY,RKZ,FKZ
WRITE(6,2200) 'FTCS',DX,DY,DZ,RKX,FKX,RKY,FKY,RKZ,FKZ
2200 FORMAT(/,5X,' *** ',A4,' ***',//
1 5X,' * DX =',F10.5,5X,'* DY =',F10.5,5X,'* DZ =',F10.5/
2 5X,' * RKX =',F10.5,' FKX =',F10.5/5X,' * RKY =',F10.5,
3 ' FKY =',F10.5/5X,' * RKZ =',F10.5,' FKZ =',F10.5 )
C
WRITE(7,2250) DT,LMT,NPR,RX,RY,RZ
WRITE(6,2250) DT,LMT,NPR,RX,RY,RZ
2250 FORMAT(/,5X,' ** DT=',F10.5,/
1 5X,' ** LMT=',I4,10X,' ** NPR=',I4,/
2 5X,' ** RX =',F10.5,5X,' ** RY =',F10.5,5X,' ** RZ =',F10.5)
C
WRITE(7,2300) TIM
WRITE(6,2300) TIM
2300 FORMAT(/'*** Initial Condition ( Y = Z = 0 ) * T =',F6.2)
C
WRITE(7,2400) 'X =','C='
WRITE(6,2400) 'X =','C='
2400 FORMAT(6X,A2,7X,A2)
C
DO 500 I = 1, NX
C
IF ( X( I).GT.-1.D-9) THEN
WRITE(7,2500) X( I), CO( I,1,1)
WRITE(6,2500) X( I), CO( I,1,1)
2500 FORMAT(5X,F6.3,3X,F9.6)
END IF
C
500 CONTINUE
C
RETURN
END
C *******************************************************************