1 : C FSTIFF SOURCE BECC 11/05/26 21:15:28 6981
2 : SUBROUTINE FSTIFF(NDIM, NPAR,
3 : & PC_L, GAM_L,
4 : & PC_R, GAM_R,
5 : & PMIN,
6 : & RHO_L, P_L, U_L, UT1_L, UT2_L, W_L,
7 : & RHO_R, P_R, U_R, UT1_R, UT2_R, W_R,
8 : & RHO_B, P_B, U_B, D_B,
9 : & RHO_A, P_A, U_A, D_A,
10 : & Z,
11 : & UINT,
12 : & F_LAG, CELLT,
13 : & LOGAN, LOGNC, LOGVAC)
14 : *
15 : *************************************************************************
16 : *
17 : * project : CAST3M, Europlexus...
18 : *
19 : * name : fstiff
20 : *
21 : * description : stiffened gas mixture.
22 : *
23 : * language : fortran 77
24 : *
25 : * author : a. beccantini den/dm2s/sfme/ltmf
26 : *
27 : *************************************************************************
28 : *
28 : * called sub :
29 : *
29 : * called by :
30 : *
31 : *************************************************************************
32 : *
33 : ***** input
34 : *
35 : * ndim = dimension (1, 2 or 3)
36 : *
37 : * npar = number of parameters involved in the
38 : * equations
39 : *
40 : * pc_, gam_ = parameters involved in the eos
41 : *
42 : * rho_, p_, u_, ut1_, ut2_, w_,
43 : * = density, pressure, normal and tangential
44 : * velocities, parameter vector
45 : * nb dim(w_) = npar
46 : *
47 : * z = value of the characteristic variable x/t in
48 : * which we want to evaluate the flux f_lag
49 : *
50 : ***** output
51 : *
52 : * uint = speed of the contact discontinuity
53 : *
54 : * f_lag = ale interfacial flux in (n,t1,t2), i.e.
55 : * rho*(un - x/t) mass
56 : * rho*(un - x/t)*un + p momentum along n
57 : * rho*(un - x/t)*ut1 momentum along t1
58 : * rho*(un - x/t)*ut2 momentum along t2
59 : * rho*(un - x/t)*ht energy
60 : * rho*(un - x/t)*w(1) parameter 1
61 : * ...
62 : *
63 : * nb
64 : * according to nkonga, comput methods appl. mech engnr 190, 2000
65 : * \dep{u}{x} + \dep{f(u)}{x} = 0
66 : * z = speed of the surface
67 : * f_z = f - z u
68 : *
69 : *
70 : * cellt = stability condition, i.e.
71 : *
72 : * dt/dx < cellt (dimension 1/velocity)
73 : *
74 : * logan = if true, anomaly detected
75 : *
76 : * lognc = if true, the Newton for the Riemann problem
77 : * solution did not converge
78 : *
79 : *************************************************************************
80 : *
81 : * 02/12/2009 created
82 : * 25/05/2011 evolution in CAST3M
83 : *
84 : *************************************************************************
85 : *
86 : * n.b.: all variables are declared
87 : *
88 : C IMPLICIT NONE
89 : IMPLICIT INTEGER(I-N)
90 : INTEGER NDIM, NPAR, IPAR
91 : REAL*8 PMIN
92 : REAL*8 PC_L, GAM_L, PC_R, GAM_R, Z
93 : & , RHO_L, P_L, U_L, UT1_L, UT2_L, W_L(NPAR)
94 : & , RHO_R, P_R, U_R, UT1_R, UT2_R, W_R(NPAR)
95 : & , F_LAG(NDIM + 3 + NPAR)
96 : & , CELLT
97 : & , D_L, D_B, D_A, D_R, RHO_B, P_B, U_B
98 : & , RHO_A, P_A, U_A
99 : & , FLURHO, FLURU, FLURT1, FLURT2, FLURET, COEF
100 : & , UINT
101 : LOGICAL LOGAN, LOGDEB, LOGNC, LOGVAC
102 : PARAMETER (LOGDEB = .FALSE.)
103 : *
104 : *************************************************************************
105 : ******** field-by-field decomposition ***********************************
106 : *************************************************************************
107 : *
107 : CALL RACSTI(
108 : & PC_L, GAM_L,
109 : & PC_R, GAM_R,
110 : & PMIN,
111 : & RHO_L, P_L, U_L,
112 : & RHO_R, P_R, U_R,
113 : & D_L, D_B, D_A, D_R,
114 : & RHO_B, P_B, U_B,
115 : & RHO_A, P_A, U_A,
116 : & LOGDEB, LOGAN, LOGNC, LOGVAC)
117 : *
118 : IF (LOGDEB .OR. LOGNC) THEN
119 : * IF (LOGDEB .OR. LOGNC .OR. LOGVAC) THEN
120 : WRITE(*,*)
121 : WRITE(*,*) 'RHO_L, P_L, U_L, D_L'
122 : WRITE(*,*) RHO_L, P_L, U_L, D_L
123 : WRITE(*,*) 'RHO_B, P_B, U_B, D_B'
124 : WRITE(*,*) RHO_B, P_B, U_B, D_B
125 : WRITE(*,*) 'RHO_A, P_A, U_A, D_A'
126 : WRITE(*,*) RHO_A, P_A, U_A, D_A
127 : WRITE(*,*) 'RHO_R, P_R, U_R, D_R'
128 : WRITE(*,*) RHO_R, P_R, U_R, D_R
129 : WRITE(*,*)
130 : ENDIF
131 : IF (LOGAN) THEN
132 : WRITE(*,*) 'SUBROUTINE FSTIFF.F'
133 : WRITE(*,*) 'ANOMALY DETECTED'
134 : GOTO 9999
135 : ENDIF
136 : *
137 : *************************************************************************
138 : ******** flux ***********************************************************
139 : *************************************************************************
140 : *
140 : CALL FLUSTI(
141 : & PC_L, GAM_L,
142 : & PC_R, GAM_R,
143 : & Z,
144 : & D_L, D_B, D_A, D_R,
145 : & RHO_L, P_L, U_L, UT1_L, UT2_L,
146 : & RHO_B, P_B, U_B,
147 : & RHO_A, P_A, U_A,
148 : & RHO_R, P_R, U_R, UT1_R, UT2_R,
149 : & UINT,
150 : & FLURHO, FLURU, FLURT1, FLURT2, FLURET,
151 : & LOGDEB, LOGAN)
152 : IF (LOGDEB) THEN
153 : C IF (LOGVAC) THEN
154 : WRITE(*,*)
155 : WRITE(*,*) 'FLURHO, FLURU, FLURT1, FLURT2, FLURET'
156 : WRITE(*,*) FLURHO, FLURU, FLURT1, FLURT2, FLURET
157 : WRITE(*,*)
158 : stop
159 : ENDIF
160 : IF (LOGAN) THEN
161 : WRITE(*,*) 'SUBROUTINE FSTIFF.F'
162 : WRITE(*,*) 'ANOMALY DETECTED'
163 : GOTO 9999
164 : ENDIF
165 : *
166 : ***** condition for the cfl
167 : *
168 : CELLT = MAX (ABS(D_L), ABS(D_R))
169 : *
170 : F_LAG(1) = FLURHO
171 : F_LAG(2) = FLURU
172 : F_LAG(3) = FLURT1
173 : IF (NDIM .EQ. 3) THEN
174 : * nb if ndim = 1 and npar = 0, index 4 does not
175 : * exist
176 : F_LAG(4) = FLURT2
177 : ENDIF
178 : F_LAG(2 + NDIM) = FLURET
179 : COEF = 1.0D0
180 : IF (FLURHO .LT. 0) THEN
181 : COEF = -1.0D0
182 : ENDIF
183 : *
184 : DO IPAR = 1, NPAR, 1
185 : F_LAG(2 + NDIM + IPAR) = 0.5D0 * FLURHO * (
186 : & ((COEF + 1.0D0) * W_L(IPAR)) +
187 : & ((1.0D0 - COEF) * W_R(IPAR))
188 : & )
189 : ENDDO
190 : 9999 CONTINUE
191 : RETURN
192 : END
193 :
194 :
© Cast3M 2003 - All rights reserved.
Disclaimer