14 | | |
15 | | I. Declarations |
16 | | 1. Wl=() |
17 | | 2. Wr=() |
18 | | 3. tol = |
19 | | |
20 | | II. Main Program Body |
21 | | 1. Call P* routine - find and record p* |
22 | | 2. Call Identification routine - find types of waves present, record |
23 | | 3. Call U*() - reads in vars and returns U* |
24 | | 4. Call rho* - reads in vars and returns rho* |
25 | | 2. Print to File - u*, rho*, p* |
26 | | |
27 | | III. Functions |
28 | | 1. FL = () |
29 | | 2. FR = () |
30 | | 3. U* = () |
31 | | 4. rho*L = {} |
32 | | 5. rho*R = {} |
33 | | |
34 | | IV. Subroutines |
35 | | 1. Find p* using initial estimate |
36 | | a. type of approximation (1)-(5) for different options (on screen options) |
37 | | b. let p0=(approximation) |
38 | | c. Begin iterative process, that ends once | | < tol |
39 | | |
40 | | 2. Determine types of waves present |
41 | | a. comparing p* to pL and pR |
42 | | b. series of if statements corresponding to flow chart |
43 | | c. left_wave = 1 for shock, 2 for rarefaction |
| 12 | PROGRAM ExactRiemannSolver |
| 13 | |
| 14 | !------------------------------------------------------------------------------------------------- |
| 15 | ! This program finds the solution of p*, u*, and rho*_L and rho*_R given initial data states wL |
| 16 | ! and wR of the Riemann problem |
| 17 | ! Author: Erica Kaminski, 12.17.12 |
| 18 | !------------------------------------------------------------------------------------------------- |
| 19 | |
| 20 | IMPLICIT NONE |
| 21 | |
| 22 | REAL :: DomLen, TimeOut, DisPos |
| 23 | INTEGER :: cells, i |
| 24 | REAL :: rho_L, p_L, u_L !Left data vars |
| 25 | REAL :: rho_R, p_R, u_R !Right data vars |
| 26 | REAL, dimension(3) :: wL !Left data state |
| 27 | REAL, dimension(3) :: wR !Right data state |
| 28 | REAL :: Cs_L, Cs_R !Sound speed Left and Right states |
| 29 | REAL, PARAMETER :: gamma = 1.4 ! Ratio specific heats |
| 30 | REAL :: p_Star, u_Star ! Output of program |
| 31 | REAL :: mpa ! Normalization constant for P |
| 32 | REAL :: u_diff, dx, xpos, s |
| 33 | REAL :: rhos, us, ps |
| 34 | |
| 35 | NAMELIST /InitialData/ DomLen, TimeOut, DisPos, cells, rho_L, p_L, u_L, rho_R, p_R, u_R, mpa |
| 36 | |
| 37 | !-------------------------------------------------------------------------------------------------- |
| 38 | ! Main body of program here |
| 39 | !-------------------------------------------------------------------------------------------------- |
| 40 | |
| 41 | |
| 42 | OPEN(UNIT = 3, File = "initial.data", status = 'old') !Initialize data |
| 43 | READ(3, NML = InitialData) |
| 44 | CLOSE(3) |
| 45 | |
| 46 | wL = (/rho_L, p_L, u_L/) !Populate w arrays |
| 47 | wR = (/rho_R, p_R, u_R/) |
| 48 | |
| 49 | Cs_L = Sqrt(Gamma*p_L/rho_L) |
| 50 | Cs_R = Sqrt(Gamma*p_R/rho_R) |
| 51 | |
| 52 | U_DIFF = u_R - u_L |
| 53 | |
| 54 | !PRESSURE POSITIVITY CONDITION --- |
| 55 | |
| 56 | IF (2.0*Cs_L/(gamma-1.0) + 2.0*Cs_R/(gamma-1.0).LE.U_DIFF) THEN |
| 57 | |
| 58 | PRINT *, 'PRESSURE POSITIVITY CONDITION VIOLATED!' |
| 59 | PRINT *, 'VACUUM STATE IS CREATED BY INITIAL CONDITIONS' |
| 60 | PRINT *, 'THIS RIEMANN SOLVER IS NOT CORRECT METHOD FOR SOLUTION' |
| 61 | PRINT *, 'PROGRAM STOPPED' |
| 62 | STOP |
| 63 | |
| 64 | END IF |
| 65 | |
| 66 | CALL Find_Pstar(p_Star, u_Star) |
| 67 | !check that p_Star = p from subprpgram with a print statement - I did, it does |
| 68 | |
| 69 | OPEN(UNIT=2, FILE='exact.out', STATUS='OLD', POSITION='APPEND') |
| 70 | |
| 71 | dx = DomLen/REAL(cells) |
| 72 | WRITE(2,*)'xpos rho u P iE' |
| 73 | DO 10 i = 1, cells |
| 74 | xpos = dx*(REAL(i)-0.5) !cell-centered |
| 75 | s = (xpos - DisPos)/TimeOut !characteristic speed |
| 76 | CALL SampleExact(p_Star, u_Star, s, rhos, us, ps) |
| 77 | WRITE(2, 20)xpos, rhos, us, ps/mpa, ps/rhos/(gamma-1.0)/mpa !the last may be internal energy...? |
| 78 | 10 Continue |
| 79 | |
| 80 | CLOSE(2) |
| 81 | |
| 82 | 20 FORMAT(5(F14.6,2X)) |
| 83 | |
| 84 | |
| 85 | |
| 86 | CONTAINS |
| 87 | |
| 88 | SUBROUTINE SampleExact(p_Star, u_Star, s, rhos, us, ps) |
| 89 | REAL,INTENT(IN) :: p_Star, u_Star, s |
| 90 | REAL,INTENT(OUT) ::rhos, us, ps !sampled density, velocity, pressure at point (x,t) |
| 91 | REAL :: SHL, STL !Speed of left rarefaction head and tail |
| 92 | REAL :: SHR, STR !Speed of right rarefaction head and tail |
| 93 | REAL :: CsL_Star, CsR_Star !Sound speed behind left and right rarefactions |
| 94 | REAL :: SL, SR !Left and right shock speeds |
| 95 | REAL :: G1, G2, G3, G4, G5, G6, ppl, ppr |
| 96 | |
| 97 | CsL_Star = Cs_L*(p_Star/p_L)**(G3/G2) |
| 98 | CsR_Star = Cs_R*(p_Star/p_R)**(G3/G2) |
| 99 | |
| 100 | G1 = gamma + 1.0 |
| 101 | G2 = 2.0*gamma |
| 102 | G3 = gamma - 1.0 |
| 103 | ppl = p_Star/p_L |
| 104 | ppr = p_Star/p_R |
| 105 | G4 = G3/G1 |
| 106 | G5 = G1/G2 |
| 107 | G6 = G3/G2 |
| 108 | |
| 109 | IF(s.LE.u_Star)THEN |
| 110 | ! (x,t) is to the left of contact |
| 111 | IF(p_L.GE.p_Star)THEN |
| 112 | !Left wave is a rarefaction, based on initial data |
| 113 | !3 possible left data states then - WL, WL_Fan, W_Star - next we determine which |
| 114 | ! PRINT*,'IT"S A RAREFACTION' |
| 115 | |
| 116 | !Are these negative quantities:? |
| 117 | SHL = u_L - Cs_L |
| 118 | STL = u_Star - CsL_Star |
| 119 | |
| 120 | IF(S.LE.SHL)THEN |
| 121 | ! PRINT *, 'LEFT DATA STATE' |
| 122 | !Left data state |
| 123 | rhos = rho_L |
| 124 | us = u_L |
| 125 | ps = p_L |
| 126 | ELSE |
| 127 | |
| 128 | ! IF((SHL.GE.S).AND.(S.LE.STL))THEN |
| 129 | IF((SHL.LE.S).AND.(S.LE.STL)) THEN |
| 130 | ! PRINT*, 'iNSIDE RARE WAVE' |
| 131 | !Inside rarefaction fan |
| 132 | rhos = rho_L*(2/G1 + G3/(G1*Cs_L)*(u_L - s))**(2/G3) |
| 133 | us = 2/G1*(Cs_L + (G3/2)*u_L + s) |
| 134 | ps = p_L*( 2/G1 + G3/(G1*Cs_L)*(u_L - s) )**(G2/G3) |
| 135 | ELSE |
| 136 | !In left star region |
| 137 | ! PRINT *, 'IN STAR REGION' |
| 138 | rhos = rho_L*ppl**(1.0/gamma) !follows from isentropic law |
| 139 | us = u_Star |
| 140 | ps = p_Star |
| 141 | END IF |
| 142 | |
| 143 | END IF |
| 144 | |
| 145 | ELSE |
| 146 | !Left wave is a shock |
| 147 | !Two states possible - pre-shock, post-shock |
| 148 | ! PRINT*,'IT"S A SHOCK ! ! !' |
| 149 | SL = (u_L - Cs_L)*SQRT(G1/G2*ppl + G3/G2) |
| 150 | IF(s.LE.SL)THEN |
| 151 | !Data state is left state/pre-shock |
| 152 | rhos = rho_L |
| 153 | us = u_L |
| 154 | ps = p_L |
| 155 | ELSE |
| 156 | !Data state is post-shock |
| 157 | rhos = rho_L*((ppl + G3/G1)/((G3/G1)*ppl + 1.0)) |
| 158 | us = u_Star |
| 159 | ps = p_Star |
| 160 | END IF |
| 161 | END IF |
| 162 | |
| 163 | ELSE |
| 164 | |
| 165 | !(x,t) is on right of contact |
| 166 | IF(p_Star.LE.p_R)THEN |
| 167 | !Rarefaction, and 3 wave states possible: wR, WRFan, W*R |
| 168 | SHR = Cs_R + u_R |
| 169 | STR = CsR_Star + u_Star |
| 170 | IF(S.LE.STR)THEN |
| 171 | !In r star region |
| 172 | rhos = rho_R*ppr**(1.0/gamma) !follows from isentropic law -- CHECK? |
| 173 | us = u_Star |
| 174 | ps = p_Star |
| 175 | ELSE |
| 176 | IF((STR.LE.S).AND.(S.LE.SHR))THEN |
| 177 | !Inside right rarefaction fan |
| 178 | rhos = rho_R*(2/G1 - G3/(G1*Cs_R)*(u_R - s))**(2/G3) |
| 179 | us = 2/G1*(-Cs_R + (G3/2)*u_R + s) |
| 180 | ps = p_R*( 2/G1 - G3/(G1*Cs_R)*(u_R - s) )**(G2/G3) |
| 181 | !Check these expressions |
| 182 | ELSE |
| 183 | !In right data state |
| 184 | rhos = rho_R |
| 185 | us = u_R |
| 186 | ps = p_R |
| 187 | END IF |
| 188 | END IF |
| 189 | ELSE |
| 190 | !Right wave is a shock |
| 191 | SR = (u_R + Cs_R)*SQRT(G1/G2*ppr + G3/G2) |
| 192 | IF(s.LE.SR)THEN |
| 193 | !In post-shock region |
| 194 | rhos = rho_R*((ppr + G3/G1)/((G3/G1)*ppr + 1.0)) |
| 195 | us = u_Star |
| 196 | ps = p_Star |
| 197 | ELSE |
| 198 | !In pre-shock region |
| 199 | rhos = rho_R |
| 200 | us = u_R |
| 201 | ps = p_R |
| 202 | END IF |
| 203 | END IF |
| 204 | END IF |
| 205 | |
| 206 | |
| 207 | END SUBROUTINE SampleExact |
| 208 | |
| 209 | |
| 210 | SUBROUTINE Find_Pstar(p, u) |
| 211 | !Does mpa NEED to be read into this subroutine, if it is already initialized in main body of program? |
| 212 | !----------------------------------------------------------------------------------------------- |
| 213 | ! Implements Newton's iterative method to solve for p_star |
| 214 | !----------------------------------------------------------------------------------------------- |
| 215 | REAL, INTENT(OUT) :: p, u |
| 216 | REAL :: p0 ! Initial estimate of P_star |
| 217 | REAL, PARAMETER :: tol = 0.0000001 |
| 218 | REAL :: change ! Relative change in pressure due to iteration |
| 219 | 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 |
| 220 | INTEGER :: i, NmrItr |
| 221 | REAL :: f_p, f_pD, FL, FLD, FR, FRD ! pressure function in Toro |
| 222 | REAL :: Delta |
| 223 | NAMELIST /Itr/ NmrItr |
| 224 | |
| 225 | !statement labels are not global: |
| 226 | OPEN(UNIT = 3, FILE = 'initial.data') |
| 227 | READ(3, NML = Itr) |
| 228 | CLOSE(3) |
| 229 | |
| 230 | |
| 231 | !Get estimate for P_star |
| 232 | CALL P_estimate(p0) |
| 233 | |
| 234 | |
| 235 | pOld = p0 ! Initialize pOld |
| 236 | |
| 237 | !Check that the vars are global: |
| 238 | ! PRINT *, 'u_L, p_L, rho_L, cs_L =', u_L, p_L, rho_L, cs_L |
| 239 | |
| 240 | OPEN(UNIT = 2, FILE = 'exact.out', STATUS = 'UNKNOWN') |
| 241 | WRITE(2,*) '--------------------------------------------------------------------------' |
| 242 | WRITE(2,*) ' Iteration number Change ' |
| 243 | WRITE(2,*) '--------------------------------------------------------------------------' |
| 244 | |
| 245 | !Need define Cs_L here too or main program enough? |
| 246 | DO 10 i = 1,NmrItr |
| 247 | CALL PressFunct(FL, FLD, pOld, rho_L, p_L, Cs_L) |
| 248 | CALL PressFunct(FR, FRD, pOld, rho_R, p_R, Cs_R) |
| 249 | f_p = FL + FR + U_diff |
| 250 | f_pD = FLD + FRD |
| 251 | Delta = f_p/f_pD |
| 252 | pNew = pOld - Delta ! This is corrected p. |
| 253 | change = 2.0*ABS((pNew - pOld)/(pNew + pOld)) |
| 254 | WRITE(2,*) i, change |
| 255 | IF (change.LE.tol) GOTO 20 |
| 256 | IF (pNew.LT.0.0) pNew = tol |
| 257 | pOld = pNew |
| 258 | 10 Continue |
| 259 | |
| 260 | WRITE(2,*) 'Newton Rhapson Method Diverged!' |
| 261 | |
| 262 | 20 WRITE(2,*) 'Pressure in Star Region and U in star region' |
| 263 | p = pNew/MPA |
| 264 | u = (1.0/2.0)*(u_L + u_R) + (1.0/2.0)*(fR - fL) |
| 265 | WRITE(2,*) 'p_Star = ', p |
| 266 | WRITE(2,*) 'u_Star = ', u |
| 267 | |
| 268 | CLOSE(2) |
| 269 | |
| 270 | END SUBROUTINE Find_Pstar |
| 271 | |
| 272 | SUBROUTINE P_estimate(pEst) |
| 273 | !Employs the subroutine from Toro for estimating initial P_star |
| 274 | |
| 275 | REAL,INTENT(OUT) :: pEst |
| 276 | REAL :: G1, G2, G3, G4, G5, G6, G7, G8 |
| 277 | REAL :: CUP, GEL, GER, pMax, pMin, PPV, PQ, PTL, PTR, QMAX, QUSER, UM |
| 278 | |
| 279 | G1 = (gamma - 1.0)/(2.0*gamma) |
| 280 | G2 = (gamma + 1.0)/(2.0*gamma) |
| 281 | G3 = 2.0*gamma/(gamma - 1.0) |
| 282 | G4 = 2.0/(gamma - 1.0) |
| 283 | G5 = 2.0/(gamma + 1.0) |
| 284 | G6 = (gamma - 1.0)/(gamma + 1.0) |
| 285 | G7 = (gamma - 1.0)/2.0 |
| 286 | G8 = (gamma - 1.0) |
| 287 | |
| 288 | QUSER = 2.0 |
| 289 | |
| 290 | ! Guess P using the PVRS solver |
| 291 | |
| 292 | CUP = 0.25*(rho_L + rho_R)*(Cs_L + Cs_R) |
| 293 | PPV = 0.5*(p_L + p_R) + 0.5*(u_L - u_R)*CUP |
| 294 | PPV = MAX(0.0, PPV) |
| 295 | PMIN = MIN(p_L, p_R) |
| 296 | PMAX = MAX(p_L, p_R) |
| 297 | QMAX = PMAX/PMIN |
| 298 | |
| 299 | IF (QMAX.LE.QUSER.AND.(PMIN.LE.PPV.AND.PPV.LE.PMAX)) THEN |
| 300 | !SELECT PVRS SOLVER |
| 301 | pEst = PPV |
| 302 | |
| 303 | ELSE |
| 304 | |
| 305 | IF(PPV.LT.PMIN) THEN |
| 306 | !SELECT 2-RAREFACTION RS |
| 307 | PQ = (p_L/p_R)**G1 |
| 308 | UM = (PQ*u_L/Cs_L + u_R/Cs_R + G4*(PQ - 1.0))/(PQ/Cs_L + 1.0/Cs_R) |
| 309 | PTL = 1.0 + G7*(u_L - UM)/Cs_L |
| 310 | PTR = 1.0 + G7*(UM - u_R)/Cs_R |
| 311 | pEst = 0.5*(p_L*PTL**G3 + p_R*PTR**G3) |
| 312 | |
| 313 | ELSE |
| 314 | !Select the 2-shock RS with PVRS as estimate |
| 315 | |
| 316 | GEL = SQRT((G5/rho_L)/(G6*p_L + PPV)) |
| 317 | GER = SQRT((G5/rho_R)/(G6*p_R + PPV)) |
| 318 | pEst = (GEL*p_L + GER*p_R - (u_R - u_L))/(GEL + GER) |
| 319 | |
| 320 | END IF |
| 321 | |
| 322 | END IF |
| 323 | |
| 324 | |
| 325 | END SUBROUTINE P_estimate |
| 326 | |
| 327 | SUBROUTINE PressFunct(f, fD, pItr, rho, p, Cs) |
| 328 | |
| 329 | |
| 330 | REAL, INTENT(OUT) :: f, fd ! Pressure function and its derivative in Toro |
| 331 | REAL, INTENT(IN) :: pItr ! Pressure correction as computed by Find_Pstar |
| 332 | REAL, INTENT(IN) :: rho, p, Cs !Initial data |
| 333 | REAL :: A, B ! Constants in pressure function |
| 334 | |
| 335 | |
| 336 | A = 2.0/((gamma+1.0)*rho) |
| 337 | B = ((gamma-1.0)/(gamma+1.0))*p |
| 338 | |
| 339 | |
| 340 | IF(pItr.LT.p)THEN !Rarefaction wave |
| 341 | |
| 342 | f = (2.0*Cs/(gamma - 1.0))*((pItr/p)**((gamma-1.0)/(2.0*gamma)) - 1.0) |
| 343 | fd = (1.0/(rho*Cs))*(pItr/p)**(-(gamma+1)/(2*gamma)) |
| 344 | |
| 345 | ELSE |
| 346 | |
| 347 | !Shock |
| 348 | f = (pItr - p)*SQRT(A/(pItr+B)) |
| 349 | fd = Sqrt(A /(B + pItr))*(1 - (pItr - p)/(2*(B+pItr))) |
| 350 | |
| 351 | END IF |
| 352 | |
| 353 | END SUBROUTINE PressFunct |
| 354 | |
| 355 | END PROGRAM ExactRiemannSolver |
| 356 | |