1 : C RACSTI SOURCE BECC 11/05/26 21:16:16 6981 2 : SUBROUTINE RACSTI( 3 : & PC_L, GAM_L, 4 : & PC_R, GAM_R, 5 : & PMIN, 6 : & RHO_L, P_L, U_L, 7 : & RHO_R, P_R, U_R, 8 : & D_L, D_B, D_A, D_R, 9 : & RHO_B, P_B, U_B, 10 : & RHO_A, P_A, U_A, 11 : & LOGDEB, LOGAN, LOGNC, LOGVAC) 12 : * 13 : ************************************************************************* 14 : * 15 : * project : 16 : * 17 : * name : racsti 18 : * 19 : * description : euler equations for a mixture of stiffened gases 20 : * flux in the non-reactive case. 21 : * 22 : * field-by-field decomposition (entropy-respecting). 23 : * 24 : * language : fortran 77 25 : * 26 : * author : a. beccantini den/dm2s/sfme/ltmf 27 : * 28 : ************************************************************************* 29 : * 29 : * called by : 30 : * 31 : * 32 : ************************************************************************* 33 : * 34 : ***** input 35 : * 36 : * pc_l, gam_l = properties of the gas in the left 37 : * 38 : * pc_r, gam_r = properties of the gas in the right 39 : * 40 : * rho_lr, p_lr, w_lr 41 : * = density, pressure, velocity 42 : * 43 : * pmin = if p < pmin => vacuum 44 : * 45 : * logdeb = debugging ? 46 : * 47 : ***** input / output 48 : * 49 : * logvac = vacuum ? 50 : * 51 : ***** output 52 : * 53 : * rho_ba, p_ba, w_ba 54 : * = density, pressure, velocity 55 : * 56 : * d_l, d_b, d_a, d_r 57 : * = wave speeds d_l <= d_b <= d_a <= d_r 58 : * 59 : * logan = anomaly 60 : * 61 : * lognc = if true, the Newton for the Riemann problem 62 : * solution did not converge 63 : * 64 : ************************************************************************* 65 : * 66 : * 26/11/2009 created 67 : * 25/05/2011 evolution in CAST3M 68 : * 69 : ************************************************************************* 70 : * 71 : * n.b.: all variables are declared 72 : * 73 : C IMPLICIT NONE 74 : IMPLICIT INTEGER(I-N) 75 : INTEGER ITER, NITER 76 : PARAMETER (NITER=50) 77 : REAL*8 78 : & PMIN, 79 : & PC_L, GAM_L, 80 : & PC_R, GAM_R, 81 : & D_L, D_B, D_A, D_R, 82 : & RHO_L, P_L, U_L, 83 : & RHO_B, P_B, U_B, 84 : & RHO_A, P_A, U_A, 85 : & RHO_R, P_R, U_R 86 : & , EPSERR 87 : PARAMETER (EPSERR = 1.0D-6) 88 : * 89 : REAL*8 P_VAC, A_L, A_R 90 : & , P_ES, ROES_B, UES_B, DES_B, DER_B 91 : & , ROES_A, UES_A, DES_A, DER_A 92 : & , P_INT 93 : * 94 : * debugging ? 95 : * 96 : LOGICAL LOGDEB, LOGAN, LOGNC, LOGVAC 97 : * 98 : LOGNC = .TRUE. 99 : * 100 : ***** initialisation of the speed of the wave on the left 101 : * and on the right 102 : * 102 : CALL ASTIFF(RHO_L, P_L, GAM_L, PC_L, A_L, LOGDEB, LOGAN) 102 : CALL ASTIFF(RHO_R, P_R, GAM_R, PC_R, A_R, LOGDEB, LOGAN) 103 : IF ( LOGAN ) THEN 104 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 105 : WRITE(*,*) 'INITIALISATION OF THE WAVE SPEEDS' 106 : WRITE(*,*) 'ANOMALY DETECTED' 107 : GOTO 9999 108 : ENDIF 109 : D_L = U_L - A_L 110 : D_R = U_R + A_R 111 : * 112 : ******** states 113 : * 114 : * l 115 : * *** 116 : * |* r 117 : * | * b ********* 118 : * | ******** *| 119 : * | | | * a * | 120 : * | | | ************* | 121 : * | | | | | | 122 : * | | | | | | 123 : * d_l d_b un_b un_a d_a d_r 124 : * 125 : * 126 : * if no vacuum, un_b = un_a 127 : * 128 : P_VAC = MIN( P_L, P_R ) 129 : P_VAC = MAX( P_VAC, PMIN) 130 : P_VAC = 1.0D-8 * P_VAC 130 : CALL FSTRAL( PC_L, GAM_L, P_VAC, RHO_L, P_L, U_L, 131 : & U_B, RHO_B, D_B, LOGDEB, LOGAN) 131 : CALL FSTRAR(PC_R, GAM_R, P_VAC, RHO_R, P_R, U_R, 132 : & U_A, RHO_A, D_A, LOGDEB, LOGAN) 133 : IF ( LOGAN ) THEN 134 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 135 : WRITE(*,*) 'COMPUTATION OF THE VACUUM STATE' 136 : WRITE(*,*) 'ANOMALY DETECTED' 137 : GOTO 9999 138 : ENDIF 139 : * 140 : IF (LOGVAC) THEN 141 : * 142 : * From elsewhere, we know that we are in the vacuum... 143 : * 144 : P_B = P_VAC 145 : P_A = P_VAC 146 : LOGNC = .FALSE. 147 : 148 : C 149 : ELSEIF ( U_B .LT. U_A) THEN 150 : * 151 : * vacuum 152 : * 153 : P_B = P_VAC 154 : P_A = P_VAC 155 : LOGVAC = .TRUE. 156 : LOGNC = .FALSE. 157 : ELSE 158 : LOGVAC = .FALSE. 159 : * 160 : ************************************************************************* 161 : ***** intersection (states l, b, a, r) ********************************** 162 : ************************************************************************* 163 : * 164 : * 165 : ***** initialization of p_int and of the states b, a 166 : * 167 : * 168 : P_INT = 0.5D0 * (P_R + P_L) 169 : * 169 : CALL FSTERL( PC_L, GAM_L, P_INT, RHO_L, P_L, U_L, 170 : & U_B, RHO_B, D_B, LOGDEB, LOGAN) 170 : CALL FSTERR( PC_R, GAM_R, P_INT, RHO_R, P_R, U_R, 171 : & U_A, RHO_A, D_A, LOGDEB, LOGAN) 172 : * 173 : IF ( LOGAN ) THEN 174 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 175 : WRITE(*,*) 'INITIALIZATION OF THE STATES' 176 : WRITE(*,*) 'ANOMALY DETECTED' 177 : GOTO 9999 178 : ENDIF 179 : * 180 : DO ITER = 1, NITER, 1 181 : * 182 : * pmin 183 : * 184 : P_VAC = 1D-8 * MIN (P_L, P_R, P_INT) 185 : P_VAC = MAX( P_VAC, PMIN) 186 : * 187 : IF (.TRUE.) THEN 188 : * 189 : * exact evaluation of the derivative 190 : * 190 : CALL DERL(PC_L, GAM_L, P_INT, RHO_L, P_L, U_L, DER_B, 191 : $ LOGDEB, LOGAN) 192 : IF ( LOGAN ) THEN 193 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 194 : WRITE(*,*) 'EVALUATION OF THE DERIVATES' 195 : WRITE(*,*) 'ANOMALY DETECTED' 196 : GOTO 9999 197 : ENDIF 197 : CALL DERR(PC_R, GAM_R, P_INT, RHO_R, P_R, U_R, DER_A, 198 : $ LOGDEB, LOGAN) 199 : * 200 : C WRITE(*,*) 'DER_B, DER_A ex = ', DER_B, DER_A 201 : C 202 : IF ( LOGAN ) THEN 203 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 204 : WRITE(*,*) 'EVALUATION OF THE DERIVATES' 205 : WRITE(*,*) 'ANOMALY DETECTED' 206 : GOTO 9999 207 : ENDIF 208 : ELSE 209 : * 210 : * numerical evaluation of the derivatives 211 : * 212 : P_ES = 1.001D0 * P_INT 212 : CALL FSTERL( PC_L, GAM_L, P_ES, RHO_L, P_L, U_L, 213 : & UES_B, ROES_B, DES_B, LOGDEB, LOGAN) 214 : IF ( LOGAN ) THEN 215 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 216 : WRITE(*,*) 'EVALUATION OF THE DERIVATES' 217 : WRITE(*,*) 'ANOMALY DETECTED' 218 : GOTO 9999 219 : ENDIF 220 : DER_B = (UES_B - U_B) / (P_ES - P_INT) 220 : CALL FSTERR( PC_R, GAM_R, P_ES, RHO_R, P_R, U_R, 221 : & UES_A, ROES_A, DES_A, LOGDEB, LOGAN) 222 : IF ( LOGAN ) THEN 223 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 224 : WRITE(*,*) 'EVALUATION OF THE DERIVATES' 225 : WRITE(*,*) 'ANOMALY DETECTED' 226 : GOTO 9999 227 : ENDIF 228 : DER_A = (UES_A - U_A) / (P_ES - P_INT) 229 : C WRITE(*,*) 'DER_B, DER_A num = ', DER_B, DER_A 230 : C STOP 231 : ENDIF 232 : IF ((DER_B .GT. 0.0D0)) THEN 233 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 234 : WRITE (*,*) 'DER_B =', DER_B 235 : LOGAN = .TRUE. 236 : GOTO 9999 237 : ENDIF 238 : IF ((DER_A .LT. 0.0D0)) THEN 239 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 240 : WRITE (*,*) 'DER_A =', DER_A 241 : LOGAN = .TRUE. 242 : GOTO 9999 243 : ENDIF 244 : IF ((DER_B - DER_A) .GE. 0.0D0) THEN 245 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 246 : WRITE (*,*) 'DER_B - DER_A =', (DER_B - DER_A) 247 : LOGAN = .TRUE. 248 : GOTO 9999 249 : ENDIF 250 : C 251 : P_ES = P_INT - (U_B - U_A) / (DER_B - DER_A) 252 : P_ES = MAX (P_ES, P_VAC) 252 : CALL FSTERL( PC_L, GAM_L, P_ES, RHO_L, P_L, U_L, 253 : & U_B, RHO_B, D_B, LOGDEB, LOGAN) 253 : CALL FSTERR( PC_R, GAM_R, P_ES, RHO_R, P_R, U_R, 254 : & U_A, RHO_A, D_A, LOGDEB, LOGAN) 255 : IF ( LOGAN ) THEN 256 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 257 : WRITE(*,*) 'COMPUTATION OF THE INTERSECTION' 258 : WRITE(*,*) 'ANOMALY DETECTED' 259 : GOTO 9999 260 : ENDIF 261 : IF ((ABS (P_ES - P_INT) / P_INT) .LT. EPSERR) THEN 262 : LOGNC = .FALSE. 263 : P_A = P_ES 264 : P_B = P_ES 265 : IF (P_B .GT. P_L) THEN 266 : D_L = D_B 267 : ENDIF 268 : IF (P_A .GT. P_R) THEN 269 : D_R = D_A 270 : ENDIF 271 : IF (LOGDEB) THEN 272 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 273 : WRITE(*,*) 'CONVERGENCE ACHIEVED' 274 : WRITE(*,*) 'ITER =', ITER 275 : ENDIF 276 : GOTO 9999 277 : ENDIF 278 : P_INT = P_ES 279 : IF (LOGDEB) THEN 280 : WRITE(*,*) 'SUBROUTINE RACSTI' 281 : WRITE(*,*) 'P_INT = ', P_INT 282 : WRITE(*,*) 'RHO_B, U_B', RHO_B, U_B 283 : WRITE(*,*) 'RHO_A, U_A', RHO_A, U_A 284 : ENDIF 285 : ENDDO 286 : P_A = P_ES 287 : P_B = P_ES 288 : IF (LOGNC) THEN 289 : WRITE(*,*) 'SUBROUTINE RACSTI.F' 290 : WRITE(*,*) 'CONVERGENCE NOT ACHIEVED' 291 : WRITE(*,*) 'ITER =', ITER 292 : WRITE(*,*) 'ERROR =', (ABS (P_ES - P_INT) / P_INT) 293 : ENDIF 294 : ENDIF 295 : * 296 : 9999 RETURN 297 : END 298 :
© Cast3M 2003 - All rights reserved.
Disclaimer