戻る

C**********************************************************************C
C C
C CVDF1D - FEM C
C C
C ( 1-Dimensional Transient Convection - Diffusion Equation Solver ) C
C C
C ( V.3.8 ) C
C C
C MK = 1 ( Modified Galerkin Method ) C
C = 0 ( Galerkin Method ) 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, KBC
C
DIMENSION DMT( 101,2), PMT( 101,2), DPM( 100,3),
1 BT ( 100,3), CMT( 101,2), DEL( 100,3),
2 PEL( 100,3), TMR( 100,3), NDE( 100,3),
3 EVC( 101), RVC( 101), RL ( 100),
4 BE ( 100,2), Q ( 100), C ( 101),
5 CO ( 101), KBC( 101)
C
COMMON /CST2/ MK, TSP
COMMON /CRS / RMS, AR( 101)
COMMON /BDCS/ CB( 2), NCTB( 2), NBC
COMMON /CNST/ CFD( 5), DT, LMT, U, RK, RRL
C
C ----- Attention : NNP < = 101 ---
C
OPEN ( UNIT= 7, FILE= 'CVDF1D_FEM_OUT.DATA' )
OPEN ( UNIT= 8, FILE= 'CVDF1D_ANAL.DATA' )
C
CALL TTLE ! Ex. NNP = 11, U = 0.1 K = 0.01, DT = 0.4 & 7
C
9999 CONTINUE
C
CALL SBIP ( NNP, NEL )
C
TIM = 0.D0
LOP = 0
C
CALL INIT ( RL, NDE, C, CO, Q, KBC, MB, MW, NNP, NEL )
C
C ----- Matrix ( D),( P),( DP) ---
C
CALL MTCN ( DMT, PMT, CMT, DPM, TMR, DEL, PEL, RL, BE, NDE,
1 NNP, NEL, MB, MW )
C
999 TIM = TIM + DT
LOP = LOP + 1
C
C ----- Right-Hand side Vector ---
C
IF ( MK.EQ.0) CALL VECG ( TMR, CO, RVC, NNP, MW, MB )
C
IF ( MK.EQ.1)
C
1CALL VCON ( CMT, CO, EVC, RVC, Q, RL, NDE, NNP, NEL, MB )
C
CALL SOLV ( DPM, RVC, C, BT, KBC, NNP, MW, MB, LOP )
C
CALL ANCM ( LOP, RL, NNP, NEL )
C
CALL ERCK ( C, NNP, LOP )
C
CALL POUT ( LOP, TIM, C, CO, KBC, NNP, *999 )
C
CALL FNAL ( C, NNP, LOP, *9999 )
C
STOP
END
C **********************************************************************
C
SUBROUTINE ANCM ( LOP, RL, NNP, NEL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CRS / RMS, AR( 101)
COMMON /CMPR/ RCF, X ( 101)
COMMON /BDCS/ CB ( 2), NCTB( 2), NBC
COMMON /CNST/ CFD( 5), DT, LMT, U, RK, RRL
C
DIMENSION RL( NEL)
C
PAI = 3.14159 2653 589793 D0
C
EL = 20.D0
TMX = DT* LOP
C
C --------------------------------------
DO 100 I = 1, NNP
C
AR( I) = 0.D0
100 CONTINUE
C --------------------------------------
C
AR( 1) = 1.D0
C
C ------------------------------------------------
DO 200 I = 2, NNP
C
IF ( X( I).LT.0.D0 ) AR( I) = 1.D0
IF ( ABS(X( I)).LT.1.D-4 ) AR( I) = 0.5D0
C -----------------------------------------------
C
DEM = 0.D0
C
C --------------------------------------
DO 250 K = 1, 50000

AK = 2* K - 1
DUM = AK* PAI / EL
SNE = DSIN( DUM* ( X( I) - U* TMX ))
DIM = - RCF* DUM* DUM* TMX
C
IF ( DIM .LT.-20.D0) GO TO 200
DEM = DEM + ( SNE/ AK )* DEXP( DIM)

250 CONTINUE
C --------------------------------------
C
200 AR( I) = 0.5D0 - 2.D0* DEM/ PAI
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ERCK ( C, NNP, LOP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION C( NNP)
C
COMMON /CRS / RMS, AR( 101)
COMMON /CNST/ CFD( 5), DT, LMT, U, RK, RRL

C ----- RMS : Root Mean Square ---
C
SUM = 0.D0
C
DO 100 I = 2, NNP-1
C
SUM = SUM + ( AR( I) - C( I))** 2
100 CONTINUE
C
RMS = DSQRT( SUM/ DFLOAT( NNP-2))
C
** IF ( LOP.EQ.LMT) WRITE(6,2300) ( I, C ( I),I = 1, NNP)
** IF ( LOP.EQ.LMT) WRITE(6,*) ' '
** IF ( LOP.EQ.LMT) WRITE(6,*) ' *** Analytical Solution ***'
** IF ( LOP.EQ.LMT) WRITE(6,2300) ( I, AR( I),I = 1, NNP)
C
** 2300 FORMAT(' ',10( I2,F7.3,1X))
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))
C
ERFC = 1.D0 - ERFC
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FNAL ( C, NNP, LOP, * )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CST2/ MK, TSP
COMMON /CRS / RMS, AR( 101)
COMMON /CNST/ CFD( 5), DT, LMT, U, RK, RRL
C
DIMENSION C( NNP)
C
RK = CFD( 1)
DX = RRL
Bx = U* DT/ DX
Rx = RK* DT/( DX* DX)
C
WRITE(7,2000) NNP, MK, TSP, U, DT, RK, DX, Bx, Rx, LMT
2000 FORMAT(I5,I2,7F10.3,I3)
WRITE(7,2100) LOP, RMS
2100 FORMAT(I3,D12.5)
WRITE(7,2200) ( C( I), I = 1, NNP)
2200 FORMAT(10F7.3)
C
WRITE(8,2200) ( AR( I), I = 1, NNP)
WRITE(6,2300)
2300 FORMAT(/,' Simulation to be continued ? ( 1: Yes 2: No)')
C
*** READ(5,*) KCT
KCT = 2
C
IF ( KCT.EQ.1) RETURN 1
C
RETURN
END
C **********************************************************************
C
INTEGER FUNCTION IJ1 ( I, J )
C
IMPLICIT REAL*8(A-H,O-Z)
C
IJ1 = I + J - 1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GDPMT ( DPM, DPMT, NNP, MB, MW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ------ 'DPM' from 'DPMT' ---
C
DIMENSION DPM( NNP, MW), DPMT( NNP, MB)
C
C ----------------------------------------------
DO 100 I = 1, MB
C
C --------------------------------------
DO 120 J = 1, I
C
DPM( I,J) = DPMT( J,I-J+1)
120 CONTINUE
C --------------------------------------
JS = I + 1
JE = MB + I - 1
C --------------------------------------
DO 140 J = JS, JE
C
DPM( I,J) = DPMT( I,J-I+1)
140 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
MBP1 = MB + 1
MBM1 = MB - 1
C ------------------------------------------------
DO 200 I = MBP1, NNP
C
C --------------------------------------
DO 220 J = 1, MBM1
C
DPM( I,J) = DPMT( I+J-MB, MB-J+1)
C
220 CONTINUE
C --------------------------------------
DO 240 J = MB, MW
C
DPM( I,J) = DPMT( I, J-MB+1)
C
240 CONTINUE
C --------------------------------------
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GLMT ( 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, 3), NDE( NEL, 3)
C
C -----------------------------------
C
CALL ZARY ( NNP, MB, A )
C
C -----------------------------------
DO 100 NE = 1, NEL
DO 100 I = 1, 2
C -----------------------------------
C
IG = NDE( NE,I)
C
C -----------------------------------
DO 100 J = I, 2
C
JG = NDE( NE,J)
IF( IG - JG) 120, 140, 140
120 A( IG,JG-IG+1) = A( IG,JG-IG+1) + B( NE,IJ1( I,J))
GO TO 100
C
140 A( JG,IG-JG+1) = A( JG,IG-JG+1) + B( NE,IJ1( 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 a banded unsymmetric matrix A( N,N) stored
C
C in A( N,MW) space using Gaussian eliminationmethod ---
C
C | * * * *
C | * * * * *
C | * * * * * *
C | * * * * * * *
C | * * * * * * *
C | * * * * * * * Unsymmetric banded matrix A( N,MW)
C N * ---MW---- *
C | * * * * * * *
C | * * * * * * *
C | * * * * * * *
C | * * * * * *
C | * * * * *
C | * * * *
C --MB---
C
************************************************************************
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
KMAX = MB - 1
ELSE
KMAX = 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
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
DO 120 J = K, KMAX + K
C
A( I,J) = A( I,J) - AIK* A( K,J)
120 CONTINUE
END IF
C
IF( I.GT.MB) THEN
C
DO 130 J = KK, KMAX + 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
130 CONTINUE
C
END IF
C
110 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
C ----- Back Substitution ---
C
K = N
B( K) = B( K)/ A( K, MB)
C
900 K = K - 1
C
IF ( K.LE.0) RETURN
K1 = K + 1
SUM = 0.D0
IF ( K.GE.MB ) THEN
MM = MIN0( ( N-K+MB), MW )
C
C --------------------------------------
DO 200 J = MB + 1, MM
C
K2 = K1 + J - MB - 1
SUM = SUM + A( K,J)* B( K2)
C
200 CONTINUE
C --------------------------------------
C
B( K) = ( B( K) - SUM )/ A( K,MB )
C
ELSE
C --------------------------------------
DO 300 J = K1, MB+K-1
C
SUM = SUM + A( K,J)* B( J)
300 CONTINUE
C --------------------------------------
B( K) = ( B( K) - SUM )/ A( K,K)
C
END IF
C
GO TO 900
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INIT ( RL, NDE, C, CO, Q, KBC, MB, MW, NNP, NEL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE,KBC
C
DIMENSION RL( NEL), NDE( NEL, 3), C( NNP), CO( NNP), KBC( NNP),
1 Q ( NEL)
C
COMMON /CST2/ MK, TSP
COMMON /CMPR/ RCF, X( 101)
COMMON /BDCS/ CB( 2),NCTB( 2),NBC
COMMON /CNST/ CFD( 5),DT,LMT,U,RK,RRL
C
WRITE(6,2000)
2000 FORMAT(' CFD ( Diffusion coefficient ) = ?')
C
*** READ(5,*) CFD( 1)
CFD( 1) = 0.01D0
C ======================
C
RCF = CFD( 1)
C
WRITE(6,2100)
2100 FORMAT(' U = ?')
C
*** READ(5,*) U
C
U = 0.1D0
C ==============
C
NBC = 1
NCTB( 1) = 1
CB( 1) = 1.D0
C
C ----- Computational Length ---
C
RRL = 2.D0/ DFLOAT( NNP-1)
C
WRITE(6,2200) NNP, MK, TSP, DT, LMT
2200 FORMAT(/,' *** NNP=',I3,' ** MK =',I2,' *** TSP =',F10.3,',
1 *** DT=',F7.3,' *** LOP =',I3)
C
CALL ZRVC ( NNP, Q )
C
CALL ZRVC ( NNP, CO )
C
CALL ZRVC ( NNP, C )
C
CALL ZVI2 ( NNP, KBC )
C
C ----- To set initial values ---
C
MNP = ( NNP + 1)/2
MMNP = MNP - 1
C
DO 50 I = 1, MMNP
C
CO( I) = 1.D0
50 CONTINUE
C
CO( MNP) = 0.5D0
C ( MNP) = 0.5D0
C
WRITE(6,2300) CFD( 1), U, RRL
2300 FORMAT(' * CFD( 1)=',F10.3,' * U =',F10.3,' * DX =',F10.3)
WRITE(6,2400) ( I, NCTB( I), CB( I), I = 1, NBC )
2400 FORMAT(' '/' NO. BD-NDE ''C'' '/(2I7,F8.3))
C
C --------------------------------------
DO 100 I = 1, NBC
C
C ( NCTB( I)) = CB( I)
CO ( NCTB( I)) = CB( I)
KBC( NCTB( I)) = 9999
C
100 CONTINUE
C --------------------------------------
C
RK = CFD( 1)
IF ( MK.EQ.1) CFD(1) = CFD(1) + 0.5D0* U* U* DT
IF ( MK.EQ.1) WRITE(6,2500) CFD(1)
2500 FORMAT(/,' * New Diffusion Coef.=',F9.5,/)
C
C ----------------------------
DO 200 I = 1, NEL
C
RL( I) = RRL
C
200 CONTINUE
C ----------------------------
DO 300 NE = 1, NEL
C
NDE( NE,1) = NE
NDE( NE,2) = NE + 1
NDE( NE,3) = 1
C
300 CONTINUE
C ----------------------------
C
RL( 1) = 2.D0/ DFLOAT( NNP-1)
C
DO 400 I = 1, NNP
C
X( I)= -1.D0 + RL( 1)* DFLOAT( I-1)
C
400 CONTINUE
C --------------------------------------
C
MB = 2
MW = 2* MB - 1
C ---------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MDTD ( DPM, BE, NDE, NNP, NEL, MB, MW, PPP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- To modify 'DPM' considering 'U* DC/DX' ---
C
INTEGER*2 NDE
C
DIMENSION DPM( NNP, MW), BE( NEL, 2), NDE( NEL, 3)
C
COMMON /CST2/ MK, TSP
COMMON /CNST/ CFD( 5), DT, LMT, U, RK, RRL
C
C --------------------------------------
DO 100 NE = 1, NEL
C
C ----------------------------
DO 150 I = 1, 2
C
IG = NDE( NE,I)
C
C ----------------------------
DO 150 J = 1, 2
C
JG = NDE( NE,J)
C
IF( IG.LE.MB)
1 DPM( IG,JG) = DPM( IG,JG) + DT* PPP* U* BE( NE,J)/ 2.D0
C
IF( IG.LE.MB) GO TO 150
JT = JG - IG + MB
C
DPM( IG,JT) = DPM( IG,JT) + DT* PPP* U* BE( NE,J)/ 2.D0
C
150 CONTINUE
C ----------------------------
C
100 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MTCN ( DMT, PMT, CMT, DPM, TMR, DEL, PEL,
1 RL, BE, NDE, NNP, NEL, MB, MW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- Matrix 'DMT', 'PMT' & 'DPM' ---
C
INTEGER*2 NDE
C
DIMENSION DMT( NNP, MB), PMT( NNP, MB), DPM( NNP, MW),
1 DEL ( NEL, 3), PEL ( NEL, 3), RL ( NEL),
2 BE ( NEL, 2), CMT( NNP, MB), TMR( NNP, MW),
3 NDE ( NEL, 3)
C
COMMON /CST2/ MK, TSP
COMMON /CNST/ CFD( 5), DT, LMT, U, RK, RRL
C
DO 100 NE = 1, NEL
C
BE( NE,1) = -1.D0
BE( NE,2) = 1.D0
C
CFDT = CFD( NDE( NE,3))/ RL( NE)
C
DEL( NE,1) = CFDT
DEL( NE,2) = - CFDT
DEL( NE,3) = CFDT
C
PEL( NE,1) = RL( NE)/ 3.D0
PEL( NE,2) = RL( NE)/ 6.D0
PEL( NE,3) = RL( NE)/ 3.D0
C
100 CONTINUE
C
C ---------------------------------------------------
C
CALL GLMT ( DMT, DEL, NDE, NNP, NEL, MB )
C
C ---------------------------------------------------
C
CALL GLMT ( PMT, PEL, NDE, NNP, NEL, MB )
C
C ---------------------------------------------------
DO 200 I = 1, NNP
DO 200 J = 1, MB
C
CMT( I,J) = PMT( I,J) + DT* ( TSP - 1.D0 )* DMT( I,J)
C
200 CONTINUE
C ---------------------------------------------------
DO 300 I = 1, NNP
DO 300 J = 1, MB
C
DMT( I,J) = PMT( I,J) + DT* TSP* DMT( I,J)
C
300 CONTINUE
C ---------------------------------------------------
C
CALL GDPMT ( DPM, DMT, NNP, MB, MW )
C
C ---------------------------------------------------
C
CALL GDPMT ( TMR, CMT, NNP, MB, MW )
C
C ---------------------------------------------------
IF ( MK.EQ.0) THEN
C
CALL MDTD ( DPM, BE, NDE, NNP, NEL, MB, MW, TSP )
C
TS2 = TSP - 1.D0
C
CALL MDTD ( TMR, BE, NDE, NNP, NEL, MB, MW, TS2 )
C
END IF
C ---------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE POUT ( LOP, TIM, C, CO, NBC, NNP, * )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- Stationarity Check ---
C
INTEGER*2 NBC
C
DIMENSION C( NNP), CO( NNP), NBC( NNP)
C
COMMON /CRS / RMS, AR( 101)
COMMON /CNST/ CFD( 5), DT, LMT, U, RK, RRL
C
RMXD = 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
DFRC = ( C( I) - CO( I))/ C( I)
IF ( ABS( DFRC).LT.RMXD) GO TO 100
RMXD = DABS( DFRC)
N = I
C
100 CONTINUE
C ------------------------------------------------
C
RMXD = RMXD* 100.D0
C
WRITE(6,2000) LOP, TIM, RMXD, N, RMS* 100.D0
2000 FORMAT(' LOP=',I4, 2X,'( T=',F7.3,') Vr_C =',F7.2,'(%) at Node ='
1,I3,' RMS* 100 =',F9.4)
C
** WRITE(6,2300) ( I,C( I),I = 1, NNP)
** 2300 FORMAT(' ',10( I2,F7.3,1X))
C
C --------------------------------------
DO 200 I = 1, NNP
C
CO( I) = C( I)
200 CONTINUE
C --------------------------------------
C
IF ( LOP.LT.LMT) RETURN 1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRUM ( DPM, RVC, C, NBC, M, MW, MB, LOP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NBC
C
DIMENSION DPM( M,MW), 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 900
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
DPM( I,I) = DPM( I,I)* 1.D+11
C ==================================
C
C( I) = C( I)* DPM( I,I)
GO TO 200
C
240 CONTINUE
DPM( I,MB) = DPM( I,MB)* 1.D+11
C ===================================
C
C( I) = C( I)* DPM( I,MB)
C
200 CONTINUE
C ------------------------------------------------
RETURN
C
900 CONTINUE
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)* DPM( I,I)
GO TO 300
C
340 C( I) = C( I)* DPM( I,MB)
C
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE QCON ( Q, CO, RL, U, NDE, NNP, NEL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION Q( NEL), CO( NNP), RL( NEL), NDE( NEL, 3)
C
CALL ZRVC ( NEL, Q )
C
DO 100 NE = 1, NEL
C
Q( NE) = Q( NE) + U*( CO( NDE( NE,1)) - CO( NDE( NE, 2)))/RL( NE)
C
100 CONTINUE
C
RETURN
END
C ************************************************************************
C
SUBROUTINE SBIP ( NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CST2/ MK, TSP
COMMON /CNST/ CFD( 5), DT, LMT, U, RK, RRL
C
WRITE(6,2000)
2000 FORMAT(/,' ** NNP = ? ')
*** READ(5,*) NNP
NNP = 11
C
WRITE(6,2100)
2100 FORMAT(/,' *** MK ? ( 1: Modified Gakerkin Method, 0: Galerkin Me
1thod ) ' )
*** READ(5,*) MK
MK = 0 ! GM
MK = 1 ! MGM
WRITE(6,2200) MK
2200 FORMAT(/,' ** MK =',I2 )
C
WRITE(6,2300)
2300 FORMAT(/,' Time Scheme Parameter = ? ' )
*** READ(5,*) TSP
C
TSP = 1.D0
C ---------------
C
WRITE(6,2400)
2400 FORMAT(' Time Step = ? & Loops = ?')
*** READ(5,*) DT, LMT ! DT= 0.4 LMT= 7
DT = 0.4D0
LMT = 7
C
NEL = NNP - 1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SOLV ( DPM, RVC, C, BT, NBC, M, MW, MB, LOP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- Solution of simultaneous linear equations ---
C
INTEGER*2 NBC
C
DIMENSION DPM( M, MW), RVC( M), C( M), NBC( M), BT( M, MW)
C
C -------------------------------------------------------
C
CALL PRUM ( DPM, RVC, C, NBC, M, MW, MB, LOP )
C
C -------------------------------------------------------
DO 100 I = 1, M
DO 100 J = 1, MW
C
BT( I,J) = DPM( I,J)
C
100 CONTINUE
C -------------------------------------
C
CALL GSOL ( BT, C, M, MB, MW )
C
C -------------------------------------
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 CVDF1D - FEM C'
WRITE(6,101)'C C'
WRITE(6,101)'C ( V.3.8 ) C'
WRITE(6,101)'C C'
WRITE(6,101)'C 1-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 VCON ( A, CO, EVC, RVC, Q, RL, NDE, NNP, NEL, MB )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- 'RVC' for Modified Galerkin Method ---
C
INTEGER*2 NDE
C
DIMENSION A( NNP, MB), CO( NNP), EVC( NNP), RVC( NNP),
1 Q( NEL), RL( NEL), NDE( NEL,3)
C
COMMON /CNST/ CFD( 5), DT, LMT, U, RK, RRL
C
C -----------------------------------------------------
C
CALL QCON ( Q, CO, RL, U, NDE, NNP, NEL )
C
C -----------------------------------------------------
C
CALL ZRVC ( NNP, EVC )
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
IF ( I.EQ.1) GO TO 105
JM = I - 1
IF( I.LE.MB) THEN
JS = 1
GO TO 110
END IF
JS = I - MB + 1
C
110 CONTINUE
C --------------------------------------
DO 120 J = JS, JM
C
EVC( I) = EVC( I) + A(J,I-J+1)* CO( J)
120 CONTINUE
C --------------------------------------
C
105 CONTINUE
C --------------------------------------
DO 140 J = 1, MB
C
EVC( I) = EVC( I) + A( I,J)* CO(J+I-1)
C
140 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
CALL ZRVC ( NNP, RVC )
C
C ------------------------------------------------
DO 200 NE = 1, NEL
DO 200 K = 1, 2
C
I = NDE( NE,K)
RVC( I) = RVC( I) - Q( NE)* RL( NE)/2.D0
C
200 CONTINUE
C ------------------------------------------------
DO 300 I = 1, NNP
C
RVC( I) = EVC( I) - DT* RVC( I)
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VECG ( TA, CO, RVC, NNP, MW, MB )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- 'RVC' for Galerkin Method ---
C
DIMENSION TA( NNP, MW), CO( NNP), RVC( NNP)
C
C ------------------------------------------------
DO 100 I = 1, MB
C
SUM = 0.D0
JE = MB + I - 1
C --------------------------------------
DO 150 J=1, JE
C
SUM = SUM + TA( I,J)* CO( J)
C
150 CONTINUE
C --------------------------------------
RVC( I) = SUM
C
100 CONTINUE
C ------------------------------------------------
IS = MB + 1
IE = NNP - MB
C ------------------------------------------------
DO 200 I = IS, IE
C
SUM = 0.D0
JS = I - MB + 1
JE = I + MB - 1
C --------------------------------------
DO 250 J = JS, JE
C
SUM = SUM + TA( I,J-I+MB)* CO( J)
C
250 CONTINUE
C --------------------------------------
RVC( I) = SUM
C
200 CONTINUE
C ------------------------------------------------
C
IS = NNP - MB + 1
C ------------------------------------------------
DO 300 I = IS, NNP
C
SUM = 0.D0
JS = I - MB + 1
C --------------------------------------
DO 350 J = JS, NNP
C
SUM = SUM + TA( I,J-I+MB)* CO( J)
C
350 CONTINUE
C --------------------------------------
RVC( I) = SUM
C
300 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
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ZRVC ( 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
C
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
C
100 CONTINUE
C
RETURN
END
C **********************************************************************