# u/erica/ExactRiemannSolver: solver.f90

File solver.f90, 11.8 KB (added by , 12 years ago) |
---|

Line | |
---|---|

1 | PROGRAM ExactRiemannSolver |

2 | |

3 | !------------------------------------------------------------------------------------------------- |

4 | ! This program finds the solution of p*, u*, and rho*_L and rho*_R given initial data states wL |

5 | ! and wR of the Riemann problem |

6 | ! Author: Erica Kaminski, 12.17.12 |

7 | !------------------------------------------------------------------------------------------------- |

8 | |

9 | IMPLICIT NONE |

10 | |

11 | REAL :: DomLen, TimeOut, DisPos |

12 | INTEGER :: cells, i |

13 | REAL :: rho_L, p_L, u_L !Left data vars |

14 | REAL :: rho_R, p_R, u_R !Right data vars |

15 | REAL, dimension(3) :: wL !Left data state |

16 | REAL, dimension(3) :: wR !Right data state |

17 | REAL :: Cs_L, Cs_R !Sound speed Left and Right states |

18 | REAL, PARAMETER :: gamma = 1.4 ! Ratio specific heats |

19 | REAL :: p_Star, u_Star ! Output of program |

20 | REAL :: mpa ! Normalization constant for P |

21 | REAL :: u_diff, dx, xpos, s |

22 | REAL :: rhos, us, ps |

23 | |

24 | NAMELIST /InitialData/ DomLen, TimeOut, DisPos, cells, rho_L, p_L, u_L, rho_R, p_R, u_R, mpa |

25 | |

26 | !-------------------------------------------------------------------------------------------------- |

27 | ! Main body of program here |

28 | !-------------------------------------------------------------------------------------------------- |

29 | |

30 | |

31 | OPEN(UNIT = 3, File = "initial.data", status = 'old') !Initialize data |

32 | READ(3, NML = InitialData) |

33 | CLOSE(3) |

34 | |

35 | wL = (/rho_L, p_L, u_L/) !Populate w arrays |

36 | wR = (/rho_R, p_R, u_R/) |

37 | |

38 | Cs_L = Sqrt(Gamma*p_L/rho_L) |

39 | Cs_R = Sqrt(Gamma*p_R/rho_R) |

40 | |

41 | U_DIFF = u_R - u_L |

42 | |

43 | !PRESSURE POSITIVITY CONDITION --- |

44 | |

45 | IF (2.0*Cs_L/(gamma-1.0) + 2.0*Cs_R/(gamma-1.0).LE.U_DIFF) THEN |

46 | |

47 | PRINT *, 'PRESSURE POSITIVITY CONDITION VIOLATED!' |

48 | PRINT *, 'VACUUM STATE IS CREATED BY INITIAL CONDITIONS' |

49 | PRINT *, 'THIS RIEMANN SOLVER IS NOT CORRECT METHOD FOR SOLUTION' |

50 | PRINT *, 'PROGRAM STOPPED' |

51 | STOP |

52 | |

53 | END IF |

54 | |

55 | CALL Find_Pstar(p_Star, u_Star) |

56 | !check that p_Star = p from subprpgram with a print statement - I did, it does |

57 | |

58 | OPEN(UNIT=2, FILE='exact.out', STATUS='OLD', POSITION='APPEND') |

59 | |

60 | dx = DomLen/REAL(cells) |

61 | WRITE(2,*)'xpos rho u P iE' |

62 | DO 10 i = 1, cells |

63 | xpos = dx*(REAL(i)-0.5) !cell-centered |

64 | s = (xpos - DisPos)/TimeOut !characteristic speed |

65 | CALL SampleExact(p_Star, u_Star, s, rhos, us, ps) |

66 | WRITE(2, 20)xpos, rhos, us, ps/mpa, ps/rhos/(gamma-1.0)/mpa !the last may be internal energy...? |

67 | 10 Continue |

68 | |

69 | CLOSE(2) |

70 | |

71 | 20 FORMAT(5(F14.6,2X)) |

72 | |

73 | |

74 | |

75 | CONTAINS |

76 | |

77 | SUBROUTINE SampleExact(p_Star, u_Star, s, rhos, us, ps) |

78 | REAL,INTENT(IN) :: p_Star, u_Star, s |

79 | REAL,INTENT(OUT) ::rhos, us, ps !sampled density, velocity, pressure at point (x,t) |

80 | REAL :: SHL, STL !Speed of left rarefaction head and tail |

81 | REAL :: SHR, STR !Speed of right rarefaction head and tail |

82 | REAL :: CsL_Star, CsR_Star !Sound speed behind left and right rarefactions |

83 | REAL :: SL, SR !Left and right shock speeds |

84 | REAL :: G1, G2, G3, G4, G5, G6, ppl, ppr |

85 | |

86 | CsL_Star = Cs_L*(p_Star/p_L)**(G3/G2) |

87 | CsR_Star = Cs_R*(p_Star/p_R)**(G3/G2) |

88 | |

89 | G1 = gamma + 1.0 |

90 | G2 = 2.0*gamma |

91 | G3 = gamma - 1.0 |

92 | ppl = p_Star/p_L |

93 | ppr = p_Star/p_R |

94 | G4 = G3/G1 |

95 | G5 = G1/G2 |

96 | G6 = G3/G2 |

97 | |

98 | IF(s.LE.u_Star)THEN |

99 | ! (x,t) is to the left of contact |

100 | IF(p_L.GE.p_Star)THEN |

101 | !Left wave is a rarefaction, based on initial data |

102 | !3 possible left data states then - WL, WL_Fan, W_Star - next we determine which |

103 | ! PRINT*,'IT"S A RAREFACTION' |

104 | |

105 | !Are these negative quantities:? |

106 | SHL = u_L - Cs_L |

107 | STL = u_Star - CsL_Star |

108 | |

109 | IF(S.LE.SHL)THEN |

110 | ! PRINT *, 'LEFT DATA STATE' |

111 | !Left data state |

112 | rhos = rho_L |

113 | us = u_L |

114 | ps = p_L |

115 | ELSE |

116 | |

117 | ! IF((SHL.GE.S).AND.(S.LE.STL))THEN |

118 | IF((SHL.LE.S).AND.(S.LE.STL)) THEN |

119 | ! PRINT*, 'iNSIDE RARE WAVE' |

120 | !Inside rarefaction fan |

121 | rhos = rho_L*(2/G1 + G3/(G1*Cs_L)*(u_L - s))**(2/G3) |

122 | us = 2/G1*(Cs_L + (G3/2)*u_L + s) |

123 | ps = p_L*( 2/G1 + G3/(G1*Cs_L)*(u_L - s) )**(G2/G3) |

124 | ELSE |

125 | !In left star region |

126 | ! PRINT *, 'IN STAR REGION' |

127 | rhos = rho_L*ppl**(1.0/gamma) !follows from isentropic law |

128 | us = u_Star |

129 | ps = p_Star |

130 | END IF |

131 | |

132 | END IF |

133 | |

134 | ELSE |

135 | !Left wave is a shock |

136 | !Two states possible - pre-shock, post-shock |

137 | ! PRINT*,'IT"S A SHOCK ! ! !' |

138 | SL = (u_L - Cs_L)*SQRT(G1/G2*ppl + G3/G2) |

139 | IF(s.LE.SL)THEN |

140 | !Data state is left state/pre-shock |

141 | rhos = rho_L |

142 | us = u_L |

143 | ps = p_L |

144 | ELSE |

145 | !Data state is post-shock |

146 | rhos = rho_L*((ppl + G3/G1)/((G3/G1)*ppl + 1.0)) |

147 | us = u_Star |

148 | ps = p_Star |

149 | END IF |

150 | END IF |

151 | |

152 | ELSE |

153 | |

154 | !(x,t) is on right of contact |

155 | IF(p_Star.LE.p_R)THEN |

156 | !Rarefaction, and 3 wave states possible: wR, WRFan, W*R |

157 | SHR = Cs_R + u_R |

158 | STR = CsR_Star + u_Star |

159 | IF(S.LE.STR)THEN |

160 | !In r star region |

161 | rhos = rho_R*ppr**(1.0/gamma) !follows from isentropic law -- CHECK? |

162 | us = u_Star |

163 | ps = p_Star |

164 | ELSE |

165 | IF((STR.LE.S).AND.(S.LE.SHR))THEN |

166 | !Inside right rarefaction fan |

167 | rhos = rho_R*(2/G1 - G3/(G1*Cs_R)*(u_R - s))**(2/G3) |

168 | us = 2/G1*(-Cs_R + (G3/2)*u_R + s) |

169 | ps = p_R*( 2/G1 - G3/(G1*Cs_R)*(u_R - s) )**(G2/G3) |

170 | !Check these expressions |

171 | ELSE |

172 | !In right data state |

173 | rhos = rho_R |

174 | us = u_R |

175 | ps = p_R |

176 | END IF |

177 | END IF |

178 | ELSE |

179 | !Right wave is a shock |

180 | SR = (u_R + Cs_R)*SQRT(G1/G2*ppr + G3/G2) |

181 | IF(s.LE.SR)THEN |

182 | !In post-shock region |

183 | rhos = rho_R*((ppr + G3/G1)/((G3/G1)*ppr + 1.0)) |

184 | us = u_Star |

185 | ps = p_Star |

186 | ELSE |

187 | !In pre-shock region |

188 | rhos = rho_R |

189 | us = u_R |

190 | ps = p_R |

191 | END IF |

192 | END IF |

193 | END IF |

194 | |

195 | |

196 | END SUBROUTINE SampleExact |

197 | |

198 | |

199 | SUBROUTINE Find_Pstar(p, u) |

200 | !Does mpa NEED to be read into this subroutine, if it is already initialized in main body of program? |

201 | !----------------------------------------------------------------------------------------------- |

202 | ! Implements Newton's iterative method to solve for p_star |

203 | !----------------------------------------------------------------------------------------------- |

204 | REAL, INTENT(OUT) :: p, u |

205 | REAL :: p0 ! Initial estimate of P_star |

206 | REAL, PARAMETER :: tol = 0.0000001 |

207 | REAL :: change ! Relative change in pressure due to iteration |

208 | REAL :: pOld, pNew! Cs_L, Cs_R, rho_L, rho_R, p_L, p_R, u_L, u_R -- these should be all accessible as they are global vars |

209 | INTEGER :: i, NmrItr |

210 | REAL :: f_p, f_pD, FL, FLD, FR, FRD ! pressure function in Toro |

211 | REAL :: Delta |

212 | NAMELIST /Itr/ NmrItr |

213 | |

214 | !statement labels are not global: |

215 | OPEN(UNIT = 3, FILE = 'initial.data') |

216 | READ(3, NML = Itr) |

217 | CLOSE(3) |

218 | |

219 | |

220 | !Get estimate for P_star |

221 | CALL P_estimate(p0) |

222 | |

223 | |

224 | pOld = p0 ! Initialize pOld |

225 | |

226 | !Check that the vars are global: |

227 | ! PRINT *, 'u_L, p_L, rho_L, cs_L =', u_L, p_L, rho_L, cs_L |

228 | |

229 | OPEN(UNIT = 2, FILE = 'exact.out', STATUS = 'UNKNOWN') |

230 | WRITE(2,*) '--------------------------------------------------------------------------' |

231 | WRITE(2,*) ' Iteration number Change ' |

232 | WRITE(2,*) '--------------------------------------------------------------------------' |

233 | |

234 | !Need define Cs_L here too or main program enough? |

235 | DO 10 i = 1,NmrItr |

236 | CALL PressFunct(FL, FLD, pOld, rho_L, p_L, Cs_L) |

237 | CALL PressFunct(FR, FRD, pOld, rho_R, p_R, Cs_R) |

238 | f_p = FL + FR + U_diff |

239 | f_pD = FLD + FRD |

240 | Delta = f_p/f_pD |

241 | pNew = pOld - Delta ! This is corrected p. |

242 | change = 2.0*ABS((pNew - pOld)/(pNew + pOld)) |

243 | WRITE(2,*) i, change |

244 | IF (change.LE.tol) GOTO 20 |

245 | IF (pNew.LT.0.0) pNew = tol |

246 | pOld = pNew |

247 | 10 Continue |

248 | |

249 | WRITE(2,*) 'Newton Rhapson Method Diverged!' |

250 | |

251 | 20 WRITE(2,*) 'Pressure in Star Region and U in star region' |

252 | p = pNew/MPA |

253 | u = (1.0/2.0)*(u_L + u_R) + (1.0/2.0)*(fR - fL) |

254 | WRITE(2,*) 'p_Star = ', p |

255 | WRITE(2,*) 'u_Star = ', u |

256 | |

257 | CLOSE(2) |

258 | |

259 | END SUBROUTINE Find_Pstar |

260 | |

261 | SUBROUTINE P_estimate(pEst) |

262 | !Employs the subroutine from Toro for estimating initial P_star |

263 | |

264 | REAL,INTENT(OUT) :: pEst |

265 | REAL :: G1, G2, G3, G4, G5, G6, G7, G8 |

266 | REAL :: CUP, GEL, GER, pMax, pMin, PPV, PQ, PTL, PTR, QMAX, QUSER, UM |

267 | |

268 | G1 = (gamma - 1.0)/(2.0*gamma) |

269 | G2 = (gamma + 1.0)/(2.0*gamma) |

270 | G3 = 2.0*gamma/(gamma - 1.0) |

271 | G4 = 2.0/(gamma - 1.0) |

272 | G5 = 2.0/(gamma + 1.0) |

273 | G6 = (gamma - 1.0)/(gamma + 1.0) |

274 | G7 = (gamma - 1.0)/2.0 |

275 | G8 = (gamma - 1.0) |

276 | |

277 | QUSER = 2.0 |

278 | |

279 | ! Guess P using the PVRS solver |

280 | |

281 | CUP = 0.25*(rho_L + rho_R)*(Cs_L + Cs_R) |

282 | PPV = 0.5*(p_L + p_R) + 0.5*(u_L - u_R)*CUP |

283 | PPV = MAX(0.0, PPV) |

284 | PMIN = MIN(p_L, p_R) |

285 | PMAX = MAX(p_L, p_R) |

286 | QMAX = PMAX/PMIN |

287 | |

288 | IF (QMAX.LE.QUSER.AND.(PMIN.LE.PPV.AND.PPV.LE.PMAX)) THEN |

289 | !SELECT PVRS SOLVER |

290 | pEst = PPV |

291 | |

292 | ELSE |

293 | |

294 | IF(PPV.LT.PMIN) THEN |

295 | !SELECT 2-RAREFACTION RS |

296 | PQ = (p_L/p_R)**G1 |

297 | UM = (PQ*u_L/Cs_L + u_R/Cs_R + G4*(PQ - 1.0))/(PQ/Cs_L + 1.0/Cs_R) |

298 | PTL = 1.0 + G7*(u_L - UM)/Cs_L |

299 | PTR = 1.0 + G7*(UM - u_R)/Cs_R |

300 | pEst = 0.5*(p_L*PTL**G3 + p_R*PTR**G3) |

301 | |

302 | ELSE |

303 | !Select the 2-shock RS with PVRS as estimate |

304 | |

305 | GEL = SQRT((G5/rho_L)/(G6*p_L + PPV)) |

306 | GER = SQRT((G5/rho_R)/(G6*p_R + PPV)) |

307 | pEst = (GEL*p_L + GER*p_R - (u_R - u_L))/(GEL + GER) |

308 | |

309 | END IF |

310 | |

311 | END IF |

312 | |

313 | |

314 | END SUBROUTINE P_estimate |

315 | |

316 | SUBROUTINE PressFunct(f, fD, pItr, rho, p, Cs) |

317 | |

318 | |

319 | REAL, INTENT(OUT) :: f, fd ! Pressure function and its derivative in Toro |

320 | REAL, INTENT(IN) :: pItr ! Pressure correction as computed by Find_Pstar |

321 | REAL, INTENT(IN) :: rho, p, Cs !Initial data |

322 | REAL :: A, B ! Constants in pressure function |

323 | |

324 | |

325 | A = 2.0/((gamma+1.0)*rho) |

326 | B = ((gamma-1.0)/(gamma+1.0))*p |

327 | |

328 | |

329 | IF(pItr.LT.p)THEN !Rarefaction wave |

330 | |

331 | f = (2.0*Cs/(gamma - 1.0))*((pItr/p)**((gamma-1.0)/(2.0*gamma)) - 1.0) |

332 | fd = (1.0/(rho*Cs))*(pItr/p)**(-(gamma+1)/(2*gamma)) |

333 | |

334 | ELSE |

335 | |

336 | !Shock |

337 | f = (pItr - p)*SQRT(A/(pItr+B)) |

338 | fd = Sqrt(A /(B + pItr))*(1 - (pItr - p)/(2*(B+pItr))) |

339 | |

340 | END IF |

341 | |

342 | END SUBROUTINE PressFunct |

343 | |

344 | END PROGRAM ExactRiemannSolver |