-
Notifications
You must be signed in to change notification settings - Fork 11
/
riemann_exact.py
177 lines (131 loc) · 5.66 KB
/
riemann_exact.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
"""An exact Riemann solver for the Euler equations with a gamma-law
gas. The left and right states are stored as State objects. We then
create a RiemannProblem object with the left and right state:
> rp = RiemannProblem(left_state, right_state)
Next we solve for the star state:
> rp.find_star_state()
Finally, we sample the solution to find the interface state, which
is returned as a State object:
> q_int = rp.sample_solution()
"""
import numpy as np
import scipy.optimize as optimize
class State(object):
""" a simple object to hold a primitive variable state """
def __init__(self, p=1.0, u=0.0, rho=1.0):
self.p = p
self.u = u
self.rho = rho
def __str__(self):
return "rho: {}; u: {}; p: {}".format(self.rho, self.u, self.p)
class RiemannProblem(object):
""" a class to define a Riemann problem. It takes a left
and right state. Note: we assume a constant gamma """
def __init__(self, left_state, right_state, gamma=1.4):
self.left = left_state
self.right = right_state
self.gamma = gamma
self.ustar = None
self.pstar = None
def u_hugoniot(self, p, side):
"""define the Hugoniot curve, u(p)."""
if side == "left":
state = self.left
s = 1.0
elif side == "right":
state = self.right
s = -1.0
c = np.sqrt(self.gamma*state.p/state.rho)
if p < state.p:
# rarefaction
u = state.u + s*(2.0*c/(self.gamma-1.0))* \
(1.0 - (p/state.p)**((self.gamma-1.0)/(2.0*self.gamma)))
else:
# shock
beta = (self.gamma+1.0)/(self.gamma-1.0)
u = state.u + s*(2.0*c/np.sqrt(2.0*self.gamma*(self.gamma-1.0)))* \
(1.0 - p/state.p)/np.sqrt(1.0 + beta*p/state.p)
return u
def find_star_state(self, p_min=0.001, p_max=1000.0):
""" root find the Hugoniot curve to find ustar, pstar """
# we need to root-find on
self.pstar = optimize.brentq(
lambda p: self.u_hugoniot(p, "left") - self.u_hugoniot(p, "right"),
p_min, p_max)
self.ustar = self.u_hugoniot(self.pstar, "left")
def shock_solution(self, sgn, state):
"""return the interface solution considering a shock"""
p_ratio = self.pstar/state.p
c = np.sqrt(self.gamma*state.p/state.rho)
# Toro, eq. 4.52 / 4.59
S = state.u + sgn*c*np.sqrt(0.5*(self.gamma + 1.0)/self.gamma*p_ratio +
0.5*(self.gamma - 1.0)/self.gamma)
# are we to the left or right of the shock?
if (self.ustar < 0 and S < 0) or (self.ustar > 0 and S > 0):
# R/L region
solution = state
else:
# * region -- get rhostar from Toro, eq. 4.50 / 4.57
gam_fac = (self.gamma - 1.0)/(self.gamma + 1.0)
rhostar = state.rho * (p_ratio + gam_fac)/(gam_fac * p_ratio + 1.0)
solution = State(rho=rhostar, u=self.ustar, p=self.pstar)
return solution
def rarefaction_solution(self, sgn, state):
"""return the interface solution considering a rarefaction wave"""
# find the speed of the head and tail of the rarefaction fan
# isentropic (Toro eq. 4.54 / 4.61)
p_ratio = self.pstar/state.p
c = np.sqrt(self.gamma*state.p/state.rho)
cstar = c*p_ratio**((self.gamma-1.0)/(2*self.gamma))
lambda_head = state.u + sgn*c
lambda_tail = self.ustar + sgn*cstar
gam_fac = (self.gamma - 1.0)/(self.gamma + 1.0)
if (sgn > 0 and lambda_head < 0) or (sgn < 0 and lambda_head > 0):
# R/L region
solution = state
elif (sgn > 0 and lambda_tail > 0) or (sgn < 0 and lambda_tail < 0):
# * region, we use the isentropic density (Toro 4.53 / 4.60)
solution = State(rho = state.rho*p_ratio**(1.0/self.gamma),
u = self.ustar, p = self.pstar)
else:
# we are in the fan -- Toro 4.56 / 4.63
rho = state.rho * (2/(self.gamma + 1.0) -
sgn*gam_fac*state.u/c)**(2.0/(self.gamma-1.0))
u = 2.0/(self.gamma + 1.0) * ( -sgn*c + 0.5*(self.gamma - 1.0)*state.u)
p = state.p * (2/(self.gamma + 1.0) -
sgn*gam_fac*state.u/c)**(2.0*self.gamma/(self.gamma-1.0))
solution = State(rho=rho, u=u, p=p)
return solution
def sample_solution(self):
"""given the star state (ustar, pstar), find the state on the interface"""
if self.ustar < 0:
# we are in the R* or R region
state = self.right
sgn = 1.0
else:
# we are in the L* or L region
state = self.left
sgn = -1.0
# is the non-contact wave a shock or rarefaction?
if self.pstar > state.p:
# compression! we are a shock
solution = self.shock_solution(sgn, state)
else:
# rarefaction
solution = self.rarefaction_solution(sgn, state)
return solution
def cons_flux(state, v):
""" given an interface state, return the conservative flux"""
flux = np.zeros((v.nvar), dtype=np.float64)
flux[v.urho] = state.rho * state.u
flux[v.umx] = flux[v.urho] * state.u + state.p
flux[v.uener] = (0.5 * state.rho * state.u**2 +
state.p/(v.gamma - 1.0) + state.p) * state.u
return flux
if __name__ == "__main__":
q_l = State(rho=1.0, u=0.0, p=1.0)
q_r = State(rho=0.125, u=0.0, p=0.1)
rp = RiemannProblem(q_l, q_r, gamma=1.4)
rp.find_star_state()
q_int = rp.sample_solution()
print(q_int)