戻る


C ROGRAM MAIN
C ******************************************************************** C
C C
C BURG - FDM ( Shock Wave Simulation Solver ) C
C C
C Burgers' Equation Solver C
C C
C ( V.2.0 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION R ( 101), U( 101), UD( 101), UE( 101), AM( 5,101),
1 EF( 4), X( 101), AA( 2), BB( 2), CC( 2)
C
DIMENSION UO ( 101), UA1( 101), UA2( 101), UN1( 101),
1 UN2( 101)
C
CHARACTER*18 MTH
C
COMMON /CS2/ DFP, Q, MTH, DT, DX, LPE
COMMON /CST/ Cr, Fn, NTM, ME, JMP, JF, NNP, QC, SA, MQ, EM
C
C ----- Max. NNP = < 101 ---
C
OPEN ( 7, FILE='T_BURG_OUT.DATA' ) ! NNP= 21 Cr= 1 ( 41/ 0.3)
OPEN ( 8, FILE='T_BURG_ANAL.DATA' ) ! ( 81/ 0.1)
C
C --------------------------
C
CALL TTLE
C
C --------------------------
C
CALL DINP ( KYY )
C
C --------------------------------------------------------------------
DO 999 ME = 1, LPE
C
NK = 1
IF ( ME.EQ.3) NK = 3
C --------------------------------------------------------------------
DO 999 IJ = 1, NK
C
C ------------------------------------------------------------------
C
CALL SBRA ( X, U, UD, UE, UO, UA1, AM, AA, BB, CC, IJ,
1 KYY )
C
C ------------------------------------------------------------------
C
CALL CPLP ( U, UD, UE, UN1, UA2, UN2, EF, AA, BB, CC,
1 AM, R )
C
C ------------------------------------------------------------------
C
CALL CHCK ( U, UE )
C
999 CONTINUE
C --------------------------------------------------------------------
C
CLOSE ( 7)
CLOSE ( 8)
C
STOP
END
C **********************************************************************
C
SUBROUTINE DINP ( KYY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
CHARACTER*18 MTH
COMMON /CS2/ DFP, Q, MTH, DT, DX, LPE
COMMON /CST/ Cr, Fn, NTM, ME, JMP, JF, NNP, QC, SA, MQ, EM
C
KYY = 1
C
LPE = 4
C
WRITE(6,*) ' '
WRITE(6,*) '*** No of Nodes ? '
READ(6,*) NNP
C -------------------
C
WRITE(6,*) '**: Courant number( Cr) & Fourie number( Fn) = ? '
READ(6,*) Cr, Fn
C -----------------------
C
NTM = DFLOAT( NNP-1)/ 4.D0/ Cr
C
WRITE(6,*) ' '
WRITE(6,*) '** Input Data **'
WRITE(6,2000) NNP, NTM, Cr, Fn
2000 FORMAT(/,1X,' * No.of Nodes =',I5,' * No of Loops =',I4,
1 1X,' * Cr =',F8.5, ' * Fn =',F8.5 )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CHCK ( U, UE )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U( 101), UE( 101)
C
CHARACTER*18 MTH
COMMON /CS2/ DFP, Q, MTH, DT, DX, LPE
COMMON /CST/ Cr, Fn, NTM, ME, JMP, JF, NNP, QC, SA, MQ, EM
C
C -------------------------------------------------
SUM = 0.D0
DO 100 J = 2, JMP
C
SUM = SUM + ( U( J) - UE( J))** 2
100 CONTINUE
RMS = DSQRT( SUM/ ( DFLOAT( JMP) - 1.D0))
C -------------------------------------------------
C
IF ( ME.EQ.1) WRITE(6,*) ' '
IF ( ME.EQ.3) WRITE(6,2000) MTH, Q, RMS
2000 FORMAT(' *',A18,' ( Q =',F4.1,')',' * RMS =',D12.5)
IF ( ME.EQ.4) WRITE(6,*) ' '
IF ( ME.NE.3) WRITE(6,2100) MTH, RMS
2100 FORMAT(' *',A18,' * RMS =',D12.5,/)
C
C ------------------------------------------------
IF ( ME.EQ.LPE ) THEN
WRITE(6,*) ' *** Simulation End ***'
WRITE(8,2200) ( UE( J), J = 1, NNP)
END IF
C ------------------------------------------------
C
WRITE(7,2300) NNP, DT, DFP, DX, Cr, Fn, NTM
2300 FORMAT(I5,5F9.5,I3)
WRITE(7,2400) ME, RMS
2400 FORMAT(I3, E12.5)
WRITE(7,2200) ( U( J), J = 1, NNP)
2200 FORMAT(10F7.3)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SBRA ( X, U, UD, UE, UO, UA1, AM, AA, BB, CC,
1 IJ, KYY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION X ( 101), U ( 101), UD( 101), UE( 101),
1 AM( 5,101), AA( 2), BB( 2), CC( 2)
C
DIMENSION UO( 101), UA1( 101)
C
CHARACTER*18 MTH
COMMON /CS2/ DFP, Q, MTH, DT, DX, LPE
COMMON /CST/ Cr, Fn, NTM, ME, JMP, JF, NNP, QC, SA, MQ, EM
C
Q = 0.5D0* DFLOAT( IJ-1)
C
XMX = 1.D0
AD = 0.D0
EM = 0.D0
C
IF ( ME.EQ.4) Q = 0.D0
JF = NNP - 2
JMP = NNP - 1
C
DX = 2.D0* XMX/ DFLOAT( JMP)
DT = Cr* DX ! Courant Number :Cr = DT / DX
SA = AD* Cr* Cr
DFP = Fn* DX* DX/ DT ! Fourie Number :Fn = DFP* DT /( DX* DX)
C
IF ( DFP.LT.1.D-10) DFP = 1.D-10
C --------------------------------------
C
ATM = NTM
TIM = 0.D0
T1 = DT* ( NTM/ 2)
TMX = DT* ATM
C
IF ( ME.EQ.1) MTH = ' FTCS'
IF ( ME.EQ.2) MTH = ' Lax-Wendroff'
IF ( ME.EQ.3) MTH = ' Exp. 4pts Upwind'
IF ( ME.EQ.4) MTH = ' Crank-Nicolson'
C
IF ( KYY.NE.1) GO TO 330
C
WRITE(6,2000) DX, DT, DFP, TMX
2000 FORMAT(/,' * DX =',F5.2,' DT =',F5.2,' DFP =',D11.3,
1 ' F_Time=',F6.3)
WRITE(6,*) ' '
WRITE(6,*) ' *******************************************'
C
KYY = 0
330 CONTINUE
QC = Q* Cr/ 3.D0
IF ( ME.LT.3) QC = 0.D0
MQ = 0
IF ( DABS( QC).GT.0.0001D0) MQ = 1
IF ( ME.GT.3) GO TO 10
C
AA( 1) = 0.5D0* Cr + 3.D0* QC
BB( 1) = - 3.D0* QC
CC( 1) = - 0.5D0* Cr + QC
C
IF( ME.NE.2) GO TO 50
AA( 1) = 0.D0
BB( 1) = 0.5D0* Cr
CC( 1) = - 0.5D0* Cr
C
AA( 2) = Cr
BB( 2) = - Cr
CC( 2) = 0.D0
GO TO 50
C
10 CONTINUE
AA( 1) = EM + 0.5D0* Fn
BB( 1) = 1.D0 - 2.D0* EM -Fn
CC( 1) = EM + 0.5D0* Fn
C
50 CONTINUE
C
C ----- To initialise U & to evaluate UEX --------------------
C
CALL EXSH ( NNP, X, U, UA1, DX, XMX, DFP, T1 )
C
C ------------------------------------------------------------
C
CALL EXSH ( NNP, X, U, UE, DX, XMX, DFP, TMX )
C
C ------------------------------------------------------------
UD( 1) = U( 1)
UD( NNP) = U( NNP)
C
C --------------------------
DO 100 J = 1, NNP
C
UO( J) = U( J)
C --------------------------
DO 100 K = 1, 5
C
AM( K, J) = 0.D0
100 CONTINUE
C --------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CPLP ( U, UD, UE, UN1, UA2, UN2, EF, AA, BB,
1 CC, AM, R )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U ( 101), UD ( 101), UE( 101), UN1( 101),
1 UA2( 101), UN2( 101), EF( 4), AA ( 2),
2 BB ( 2), CC ( 2), AM( 5,101), R ( 101)
C
COMMON /CST/ Cr, Fn, NTM, ME, JMP, JF, NNP, QC, SA, MQ, EM
C
C -------------------------------------------------------------------
DO 100 N = 1, NTM
C
IF ( ME.GT.3) GO TO 50
C
C ----- Explicit Schems ---
C
IP = 1
60 CONTINUE
IF ( IP.EQ.1) EF( 1) = 0.5D0* U ( 1)* U ( 1)
IF ( IP.EQ.2) EF( 1) = 0.5D0* UD( 1)* UD( 1)
EF( 2) = EF( 1)
EF( 3) = EF( 2)
C
IF ( IP.EQ.1) EF( 4) = 0.5D0* U ( 2)* U ( 2)
IF ( IP.EQ.2) EF( 4) = 0.5D0* UD( 2)* UD( 2)
C
C ---------------------------------------------------
DO 200 J = 2, JMP
C
IF ( ME.EQ.3) EF( 1) = EF( 2)
EF( 2) = EF( 3)
EF( 3) = EF( 4)
C
IF ( IP.EQ.1) EF( 4) = 0.5D0* U ( J+1)* U ( J+1)
IF ( IP.EQ.2) EF( 4) = 0.5D0* UD( J+1)* UD( J+1)
DUM = Fn* ( U( J-1) - 2.D0* U( J) + U( J+1))
C
IF ( ME.NE.2 .OR. IP.EQ.2) GO TO 80
JP = J + 1
IF ( J.EQ.JMP) JP = J
DUM = 0.5D0* DUM + 0.5D0* Fn* ( U( J) - 2.D0* U( J+1) + U( JP+1))
1 + 0.5D0* ( U( JP) - U( J))
C
80 CONTINUE
UD( J) = U( J) + AA( IP)* EF( 2) + BB( IP)* EF( 3)
1 + CC( IP)* EF( 4) + DUM
C
IF ( ME.EQ.3) UD( J) = UD( J) - QC* EF( 1)
C
200 CONTINUE
C ---------------------------------------------------
C
IF ( ME.NE.2.OR.IP.EQ.2) GO TO 25
IP = IP + 1
GO TO 60
C
C ----- Tridiagonal system for Implicit Schemes ---
C
50 CONTINUE
IF ( MQ.EQ.1) DIM = U( 1)
C
C ------------------------------------------------
DO 300 J = 2, JMP
C
JM = J - 1
IF ( JM.NE.1) DIM = U( JM-1)
AM( 1,JM) = 0.5D0* QC* DIM
AM( 2,JM) = EM - 0.5D0* Fn - ( 0.25D0* Cr + 1.5D0* QC)* U( JM)
C
AM( 3,JM) = 1.D0 - 2.D0* EM + Fn + 1.5D0* QC* U( J)
AM( 4,JM) = EM - 0.5D0* Fn + ( 0.25D0* Cr - 0.5D0* QC)* U( J+1)
C
C ------------------
IF ( ME.NE.5 ) GO TO 350
AM( 2,JM) = AM( 2,JM) - 0.5D0* SA* U( JM)
AM( 3,JM) = AM( 3,JM) + SA* U( J)
AM( 4,JM) = AM( 4,JM) - 0.5D0* SA* U( J+1)
350 CONTINUE
C ------------------
C
R( JM) = AA( 1)* U( JM) + BB( 1)* U( J) + CC( 1)* U( J+1)
C
300 CONTINUE
C --------------------------------------------
C
R( 1) = R( 1) - AM( 2,1)* U( 1)
C
IF ( MQ.EQ.0 ) GO TO 70
R( 1) = R( 1) - AM(1,1)* U( 1)
R( 2) = R( 2) - AM(1,2)* U( 1)
C
C ----- To reduce A to Tridiagonal Form ---
C
C ---------------------------------------------
DO 400 JM = 3, JF
C
JMM = JM - 1
DUM = AM( 1,JM)/ AM( 2,JMM)
C
AM( 2,JM) =AM( 2,JM) -AM( 3,JMM)* DUM
AM( 3,JM) =AM( 3,JM) -AM( 4,JMM)* DUM
AM( 1,JM) = 0.D0
R( JM) = R( JM)- R( JMM)* DUM
C
400 CONTINUE
C ---------------------------------------------
C
AM( 1,1) = 0.D0
AM( 1,2) = 0.D0
C
70 CONTINUE
AM( 2,1) = 0.D0
AM( 4,JF) = 0.D0
C
C ---------------------------------------
C
CALL BNFC ( AM, JF, 1)
C
C ---------------------------------------
C
CALL BNSL ( R, UD, AM, JF, 1)
C
C ---------------------------------------
25 CONTINUE
C
C ---------------------------------------------
DO 500 J = 2, JMP
C
IF ( ME.LE.3) U( J) = UD( J)
IF ( ME.GT.3) U( J) = UD( J-1)
C
500 CONTINUE
C ---------------------------------------------
DO 600 J = 2, JMP
C
C -----------------
IF ( U( J).GE.2.0 .OR. U( J).LE.-1.D00) THEN
WRITE(6,*) ' '
IF ( ME.EQ.1) WRITE(6,*) ' FTCS'
IF ( ME.EQ.2) WRITE(6,*) ' Lax-Wendroff'
IF ( ME.EQ.3) WRITE(6,*) ' Exp. 4pts Upwind'
IF ( ME.EQ.4 ) WRITE(6,*) ' Crank-Nicolson'
WRITE(6,*) ' *** Divergence ! ***'
END IF
C ------------------
C
600 CONTINUE
C ---------------------------------------------
C
IF ( N.NE.NTM/ 2) GO TO 100
C ---------------------------------------------
DO 700 J = 1, NNP
C
UN1( J) = U ( J)
700 CONTINUE
C ----------------------------------------------
C
100 CONTINUE
C --------------------------------------------------------------------
C
C --------------------------------------
DO 410 J = 1, NNP
C
UA2( J) = UE( J)
UN2( J) = U ( J)
410 CONTINUE
C --------------------------------------
C
* ( U( J), J=1,NNP ) Numerical : ( UE( J), J=1,NNP ) Analytical
C
RETURN
END
C **********************************************************************
C
SUBROUTINE BNFC ( B, N, INT )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- To factorise Band Matrix arising Linear or Quadratic ---
C
DIMENSION B( 5,101)
C
IF ( INT.EQ.2) GO TO 50
C
C --- INT = 1, Linear elements = Tridaiagonal System ---
C
NP = N - 1
C ------------------------------
DO 100 J = 1, NP
C
JP = J + 1
B( 2,JP) = B( 2,JP)/ B( 3, J)
B( 3,JP) = B( 3,JP) - B( 2,JP)* B( 4,J)
C
100 CONTINUE
C -----------------------------
RETURN
C
C ---- INT = 2 : Quadratic elements ---
C
C --- To assume first equation formed at midside node ---
C
50 CONTINUE
NH = N/ 2
C
C -----------------------------------------
DO 200 I = 1, 2
C
JS = 3 - I
C
C ----------------------------------
DO 250 J = JS, NH
C
JA = 2* ( J-1)
IF ( I.EQ.2) GO TO 70
JB = JA + 2
C
B( 1,JB) = B( 1,JB)/ B( 2,JB-1)
B( 2,JB) = B( 2,JB) - B( 1,JB)* B( 3,JB-1)
B( 3,JB) = B( 3,JB) - B( 1,JB)* B( 4,JB-1)
C
GO TO 250
C
C ----- I = 2 ---
C
70 CONTINUE
C
JB = JA + 3
C
B( 2,JB-1) = B( 2,JB-1)/ B( 3,JB-2)
B( 3,JB-1) = B( 3,JB-1) - B( 2,JB-1)* B( 4,JB-2)
C
IF ( JB.GT.N) GO TO 250
C
B( 2,JB) = B( 2,JB)/ B( 3,JB-1)
B( 3,JB) = B( 3,JB) - B( 2,JB)* B( 4,JB-1)
B( 4,JB) = B( 4,JB) - B( 2,JB)* B( 5,JB-1)
C
250 CONTINUE
C ---------------------------------
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE BNSL ( R, X, B, N, INT )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- To use L, U factorisation to solve X, given R ---
C
DIMENSION R( 101), X( 101), B( 5, 101)
C
IF ( INT.EQ.2) GO TO 50
C
C ----- INT=1, R Tridiagonal system ---
C
NP = N - 1
C ------------------------------
DO 100 J = 1, NP
C
JP = J + 1
R( JP) = R( JP) - B( 2,JP)* R( J)
C
100 CONTINUE
C ------------------------------
C
X( N) = R( N)/ B( 3,N)
C ------------------------------
DO 200 J = 1, NP
C
JA = N - J
X( JA) = ( R( JA) - B( 4,JA)* X( JA+1))/ B( 3,JA)
C
200 CONTINUE
C ------------------------------
RETURN
C
50 CONTINUE
IBC = 1
NH = N/ 2
IF ( 2* NH.EQ.N) R( N+1) = 0.D0
IF ( 2* NH.EQ.N) B( 2,N+1) = 0.D0
C
C --------------------------------------------------
DO 300 I = 1, 2
C
C ---------------------------------------
DO 320 J = 1, NH
C
JA = 2*J
C
C ------------------------------
DO 340 K = 1, I
C
JB = JA - 1 + K
R( JB) = R( JB) - B( I,JB)* R( JB-1)
340 CONTINUE
C ------------------------------
C
320 CONTINUE
C ---------------------------------------
C
300 CONTINUE
C -------------------------------------------------
C
NEN = NH - IBC
X( N) = R( N)/ B( 3,N)
C
IF ( IBC.EQ.1)
1 X(N-1) = ( R( N-1) - B( 4,N-1)* X( N))/ B( 3,N-1)
C
C ----------------------------------
DO 400 J = 1, NEN
C
JA = N - 2*J + 1 - IBC
C
X( JA) = ( R( JA) - B( 4,JA)* X( JA+1) - B( 5,JA)* X( JA+2))
1 / B( 3,JA)
C
X( JA-1) = ( R( JA-1) - B( 4,JA-1)* X( JA))/ B( 3,JA-1)
C
400 CONTINUE
C -----------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE ERFC ( X, ERC )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- Complimentary Error Function ---
C
B = DABS( X)
C
IF ( B.LT.4) GO TO 50
DV = 1.D00
GO TO 100
C
50 CONTINUE
CC = DEXP( - B* B )
C
TV = 1.D0/( 1.D0 + 0.3275 911D0* B)
C
DV = 0.25483D0* TV - 0.284497D0* TV** 2 + 1.421414D0* TV** 3
1 -1.453152D0* TV** 4 + 1.061405D0* TV** 5
C
DV = 1.D0 - DV* CC
C
100 CONTINUE
C
IF ( X.LT.0.D0) DV = - DV
C
ERC = 1.D0 - DV
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXSH ( NNP, X, U, UE, DX, XMX, DFP, TMX )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ---- To set the initial U solution & Finalexact( UEX) solution ---
C
DIMENSION U( 101), UE( 101), X( 101)
C
JMP = NNP - 1
C
PI = 3.1415 921013 58979 3 D0
XST = - XMX
C
C --------------------------------
DO 100 J = 1, NNP
C
AJ = J - 1
X ( J) = XST + AJ* DX
U ( J) = 0.D0
UE( J) = 0.D0
100 CONTINUE
C --------------------------------
C
U ( 1) = 1.D0
UE( 1) = 1.D0
C
C --------------------------------------------
DO 200 J = 2, JMP
C
IF ( X( J).LT.0.D0) U( J) = 1.D00
IF ( DABS( X( J)).LT.0.0001D0) U( J) = 0.5D0
C
AJ = X( J)
XB = - AJ
XA = AJ - TMX
PQ = DSQRT( 0.25D0* PI)
C
PP = DSQRT( PI* TMX* DFP)
XB = XB* PQ/ PP
XA = XA* PQ/ PP
C
C ----- Complimentary Error Function ---
C
CALL ERFC ( XA, EXA )
C
C --------------------------------------
C
CALL ERFC ( XB, EXB )
C
C --------------------------------------
SUMA = PP* EXA
SUMB = PP* EXB
C
DUM = 0.5D0* ( AJ - 0.5D0* TMX)/ DFP
IF ( DUM.GT. 20.D0) DUM = 20.D0
IF ( DUM.LT.-20.D0) DUM = - 20.D0
SUMC = DEXP( DUM)
UE( J) = SUMA/( SUMA + SUMC* SUMB )
C
200 CONTINUE
C -----------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8(A-H,O-Z)
C
WRITE(6,*) ' '
WRITE(6,*) '*****************************************************'
WRITE(6,*) '* *'
WRITE(6,*) '* B U R G - F D M *'
WRITE(6,*) '* *'
WRITE(6,*) '* Burgers Equation Solver *'
WRITE(6,*) '* *'
WRITE(6,*) '* ( Shock Wave Analysis ) *'
WRITE(6,*) '* *'
WRITE(6,*) '* ( V.2.0 ) *'
WRITE(6,*) '* *'
WRITE(6,*) '* Copyright 2011 : Yasuhiro MATSUDA *'
WRITE(6,*) '* *'
WRITE(6,*) '*****************************************************'
C
RETURN
END
C **********************************************************************