Blog: Meeting 3/24: ERS.py

File ERS.py, 1.9 KB (added by smurugan, 11 years ago)
Line 
1import math
2import numpy as np
3import matplotlib.pyplot as plt
4
5gamma = 1.4
6xmin = -.5
7xmax = .5
8step = .001
9
10def 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
18def 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
26def 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
31def 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
40def findu(p, UL, UR):
41 return .5*(UL[1] + UR[1]) + .5*(func(p, UR) - func(p, UL))
42
43def 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
49def getFlux(rho, u, p):
50 return np.array([ rho*u, rho*u**2 + p, u*(.5*rho*u**2 + gamma*p/(gamma -1))])
51
52def 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
70if __name__ == "__main__":
71 main()