1 | import math
|
---|
2 | import numpy as np
|
---|
3 | import matplotlib.pyplot as plt
|
---|
4 |
|
---|
5 | gamma = 1.4
|
---|
6 | xmin = -.5
|
---|
7 | xmax = .5
|
---|
8 | step = .001
|
---|
9 |
|
---|
10 | def func(p, U):
|
---|
11 | A = 2/((gamma+1)*U[0])
|
---|
12 | B = (gamma - 1)/(gamma + 1) *U[2]
|
---|
13 | if p > U[2]: #shock
|
---|
14 | return ((p -U[2])*(A/(p+B))**(.5))
|
---|
15 | else: #rarefaction
|
---|
16 | return ((2*(gamma*U[2]/U[0])**(.5)/(gamma-1))*((p/U[2])**((gamma-1)/(2*gamma)) -1))
|
---|
17 |
|
---|
18 | def dfunc(p,U):
|
---|
19 | A = 2/((gamma+1)*U[0])
|
---|
20 | B = ((gamma - 1)/(gamma + 1))*U[2]
|
---|
21 | if p > U[2]: #shock
|
---|
22 | return ( 1 - (p-U[2])/(2*(B +p)))* (A/(B +p))**(.5)
|
---|
23 | else: #rarefaction
|
---|
24 | return ((gamma*U[0]*U[2])**(-.5))*((p/U[2])**(-(gamma+1)/(2*gamma)))
|
---|
25 |
|
---|
26 | def init_p(UL, UR):
|
---|
27 | num = ((gamma*UL[2]/UL[0])**(.5) + (gamma*UR[2]/UR[0])**(.5) - .5*(gamma-1)*(UR[1] - UL[1]))
|
---|
28 | den = ((gamma*UL[2]/UL[0])**(.5))*(UL[2])**((1-gamma)/(2*gamma)) + ((gamma*UR[2]/UR[0])**(.5))*(UR[2])**((1-gamma)/(2*gamma))
|
---|
29 | return (num/den)**(2*gamma/(gamma-1))
|
---|
30 |
|
---|
31 | def findp(p, UL, UR):
|
---|
32 | diff = .2
|
---|
33 | TOL = 10**(-4)
|
---|
34 | while diff > TOL:
|
---|
35 | pold = p
|
---|
36 | p += -(func(p,UR) + func(p,UL) + UR[1] - UL[1])/(dfunc(p,UR) + dfunc(p, UL))
|
---|
37 | diff = math.fabs(p - pold) / (.5 * (p + pold))
|
---|
38 | return p
|
---|
39 |
|
---|
40 | def findu(p, UL, UR):
|
---|
41 | return .5*(UL[1] + UR[1]) + .5*(func(p, UR) - func(p, UL))
|
---|
42 |
|
---|
43 | def findrho(p, U):
|
---|
44 | if p > U[2]:
|
---|
45 | return U[0]* ( (p/U[2]) + (gamma -1)/(gamma +1))/(((gamma -1)/(gamma +1))*p/U[2] +1)
|
---|
46 | else:
|
---|
47 | return U[0]*(p/U[2])**(1/gamma)
|
---|
48 |
|
---|
49 | def getFlux(rho, u, p):
|
---|
50 | return np.array([ rho*u, rho*u**2 + p, u*(.5*rho*u**2 + gamma*p/(gamma -1))])
|
---|
51 |
|
---|
52 | def ExactSolve(U, i):
|
---|
53 | if U[i,0] != U[i+1, 0]:
|
---|
54 | UL = np.array([U[i,0], U[i,1], U[i,2]])
|
---|
55 | UR = np.array([U[i+1, 0], U[i+1,1], U[i+1, 2]])
|
---|
56 | p0 = .5*(UL[0] + UR[0])
|
---|
57 | pstar = findp(p0, UL, UR)
|
---|
58 | ustar = findu(pstar, UL, UR)
|
---|
59 | if ustar > 0:
|
---|
60 | rhostar = findrho(pstar, UL)
|
---|
61 | elif ustar < 0:
|
---|
62 | rhostar = findrho(pstar, UR)
|
---|
63 | else:
|
---|
64 | rhostar = UL[0] # equivalent to UR[0]
|
---|
65 | return pstar, rhostar, ustar
|
---|
66 | else:
|
---|
67 | return 0, 0, 0
|
---|
68 |
|
---|
69 |
|
---|
70 | if __name__ == "__main__":
|
---|
71 | main()
|
---|