-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathinvestigate_wave_pattern.py
151 lines (128 loc) · 5.37 KB
/
investigate_wave_pattern.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
"""
Set up range of initial data for the reactive Riemann problem and see if the
wave pattern changes with tangential velocity.
"""
from r3d2 import eos_defns, State, RiemannProblem
from itertools import combinations
import numpy as np
def check_wave_pattern(U_l, U_r, vt_side, vts=[-0.9,0.5,0.0,-0.5,0.9]):
"""
Given an initial reactive state and burnt state, will run the
reactive Riemann problem with the reactive state having a (given) range
of tangential velocities. Shall print output to screen where the
resulting wave patterns are different.
"""
flat_patterns = make_flat_patterns(U_l, U_r, vts, vt_side)
# generate list of pairs
pairs = list(combinations(range(len(vts)), 2))
# check to see if patterns match
for i, j in pairs:
if not flat_patterns[i] == flat_patterns[j]:
print('vt = {}, {} are different'.format(vts[i], vts[j]))
print(', '.join(flat_patterns[i]))
print(', '.join(flat_patterns[j]))
def make_flat_patterns(U_l, U_r, vts, vt_side):
"""
Save some code repetition. Given reactive and burnt states, produces a list of lists of wave patterns for a given list of tangential velocities.
"""
eos_l = eos_defns.eos_gamma_law(5.0/3.0)
eos_r = eos_defns.eos_gamma_law_react(5.0/3.0, 0.1, 1.0, 1.0, eos_l)
wave_patterns = []
if vt_side == 'l':
rho_l = U_l.rho
v_l = U_l.v
eps_l = U_l.eps
eos_l = U_l.eos
else:
rho_r = U_r.rho
v_r = U_r.v
eps_r = U_r.eps
eos_r = U_r.eos
for vt in vts:
# first change the vt
if vt_side == 'l':
U_l = State(rho_l, v_l, vt, eps_l, eos_l)
#U_r = State(rho_r, v_r, 0.0, eps_r, eos_r)
else:
U_r = State(rho_r, v_r, vt, eps_r, eos_r)
rp = RiemannProblem(U_l, U_r)
wave_patterns.append(rp.waves)
# now check if all the wave patterns are the same
# flatten patterns
flat_patterns = []
for i, p in enumerate(wave_patterns):
flat_patterns.append([])
for w in p:
for s in w.wave_sections:
if not s.trivial:
flat_patterns[i].append(s.type)
return flat_patterns
def find_critical_vt(U_l, U_r, vt_side):
"""
Pretty sure that only magnitude of vt matters. As the function that
would be used for root finding is discontinuous, shall just use a very
basic bisection method. If too coarse a grid of vt values is used in
the initial pass, then a root may be missed if the wave pattern changes
to a different pattern then back again.
"""
vts = np.linspace(0., 0.9, num=100)
tolerance = 1.e-4
def bisect(vt0, vtend, tol=tolerance):
"""
Can be super lazy and only evaluate RiemannProblem for vt0 and vthalf
as if pattern(vt0) == pattern(vthalf), must be that
pattern(vthalf) != pattern(vtend) unless something has gone really
wrong.
"""
maxIts = 100
nIts = 0
while (vtend-vt0) > tol and nIts < maxIts:
vthalf = 0.5 * (vt0 + vtend)
flat_patterns = make_flat_patterns(U_l, U_r, [vt0, vthalf], vt_side)
if not flat_patterns[0] == flat_patterns[1]:
vtend = vthalf
else:
vt0 = vthalf
nIts += 1
flat_patterns = make_flat_patterns(U_l, U_r, [vt0, vtend], vt_side)
return 0.5 * (vt0 + vtend), flat_patterns
# do a first pass to find where pattern changes
flat_patterns = make_flat_patterns(U_l, U_r, vts, vt_side)
#print(flat_patterns)
critical_vts = []
critical_patterns = []
for i in range(len(vts) - 1):
if not flat_patterns[i] == flat_patterns[i+1]:
vt, pattern = bisect(vts[i], vts[i+1])
critical_vts.append(vt)
critical_patterns.append(pattern)
if len(critical_vts) == 1:
print('There is one critical tangential velocity ')
else:
print('There are {} critical tangential velocities '.format(len(critical_vts)))
for i, v in enumerate(critical_vts):
print('vt: {}, patterns: {} -> {}'.format(v, ', '.join(critical_patterns[i][0]), ', '.join(critical_patterns[i][1])))
if __name__ == "__main__":
eos = eos_defns.eos_gamma_law(5.0/3.0)
eos_reactive = eos_defns.eos_gamma_law_react(5.0/3.0, 0.1, 1.0, 1.0, eos)
#U_reactive = State(5.0, 0.0, 0.0, 2.0, eos_reactive)
# detonation wave
#U_burnt = State(8.113665227084942, -0.34940431910454606, 0.0,
# 2.7730993786742353, eos)
# cj detonation wave
#U_burnt = State(5.1558523350586452, -0.031145176327346744, 0.0,
# 2.0365206985013153, eos)
# deflagration
#U_burnt = State(0.10089486779791534, 0.97346270073482888, 0.0,
# 0.14866950243842186, eos)
# deflagration with precursor shock
#U_burnt = State(0.24316548798524526, 0.39922932397353039, 0.0,
# 0.61686385086179807, eos)
# FIXME: there is a really weird bug where this breaks if U_r has a
# non-zero normal velocity and try to give U_l tangential velocity.
U_l = State(1.0, 0.0, 0.0, 1.6, eos)
U_r = State(0.125, 0.0, 0.0, 1.2, eos_reactive)
#check_wave_pattern(U_l, U_r, 'l', vts=[-0.5,-0.1, 0.0, 0.1, 0.5, 0.87])
#check_wave_pattern(U_l, U_r, 'l', vts=[0.5, 0.55])
#check_wave_pattern(U_r, U_l, 'l')
find_critical_vt(U_l, U_r, 'l')