u/erica/ExactRiemannSolver: solver.f90

File solver.f90, 11.8 KB (added by Erica Kaminski, 12 years ago)
Line 
1PROGRAM 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...?
6710 Continue
68
69 CLOSE(2)
70
7120 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
24710 Continue
248
249 WRITE(2,*) 'Newton Rhapson Method Diverged!'
250
25120 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