plotSolution.py 8.92 KB
Newer Older
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
###############################################################################
 # This file is part of SWIFT.
 # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 # 
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published
 # by the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 # 
 # This program is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU General Public License for more details.
 # 
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 # 
 ##############################################################################

# Computes the analytical solution of the Sod shock and plots the SPH answer
 

# Generates the analytical  solution for the Sod shock test case
# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
# This follows the solution given in (Toro, 2009)


# Parameters
gas_gamma = 5./3.      # Polytropic index
rho_L = 1.             # Density left state
rho_R = 0.125          # Density right state
v_L = 0.               # Velocity left state
v_R = 0.               # Velocity right state
P_L = 1.               # Pressure left state
P_R = 0.1              # Pressure right state


import matplotlib
matplotlib.use("Agg")
from pylab import *
import h5py

# Plot parameters
params = {'axes.labelsize': 10,
'axes.titlesize': 10,
'font.size': 12,
'legend.fontsize': 12,
'xtick.labelsize': 10,
'ytick.labelsize': 10,
'text.usetex': True,
 'figure.figsize' : (9.90,6.45),
'figure.subplot.left'    : 0.045,
'figure.subplot.right'   : 0.99,
'figure.subplot.bottom'  : 0.05,
'figure.subplot.top'     : 0.99,
'figure.subplot.wspace'  : 0.15,
'figure.subplot.hspace'  : 0.12,
'lines.markersize' : 6,
'lines.linewidth' : 3.,
'text.latex.unicode': True
}
rcParams.update(params)
rc('font',**{'family':'sans-serif','sans-serif':['Times']})


snap = int(sys.argv[1])


# Read the simulation data
70
sim = h5py.File("sodShock_%04d.hdf5"%snap, "r")
71
72
73
74
75
76
boxSize = sim["/Header"].attrs["BoxSize"][0]
time = sim["/Header"].attrs["Time"][0]
scheme = sim["/HydroScheme"].attrs["Scheme"]
kernel = sim["/HydroScheme"].attrs["Kernel function"]
neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
eta = sim["/HydroScheme"].attrs["Kernel eta"]
77
git = sim["Code"].attrs["Git Revision"]
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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227

x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0]
u = sim["/PartType0/InternalEnergy"][:]
S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:]
rho = sim["/PartType0/Density"][:]

N = 1000  # Number of points
x_min = -1.
x_max = 1.

x += x_min

# ---------------------------------------------------------------
# Don't touch anything after this.
# ---------------------------------------------------------------

c_L = sqrt(gas_gamma * P_L / rho_L)   # Speed of the rarefaction wave
c_R = sqrt(gas_gamma * P_R / rho_R)   # Speed of the shock front

# Helpful variable
Gama = (gas_gamma - 1.) / (gas_gamma + 1.)
beta = (gas_gamma - 1.) / (2. * gas_gamma)

# Characteristic function and its derivative, following Toro (2009)
def compute_f(P_3, P, c):
    u = P_3 / P
    if u > 1:
        term1 = gas_gamma*((gas_gamma+1.)*u + gas_gamma-1.)
        term2 = sqrt(2./term1)
        fp = (u - 1.)*c*term2
        dfdp = c*term2/P + (u - 1.)*c/term2*(-1./term1**2)*gas_gamma*(gas_gamma+1.)/P
    else:
        fp = (u**beta - 1.)*(2.*c/(gas_gamma-1.))
        dfdp = 2.*c/(gas_gamma-1.)*beta*u**(beta-1.)/P
    return (fp, dfdp)

# Solution of the Riemann problem following Toro (2009) 
def RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R):
    P_new = ((c_L + c_R + (v_L - v_R)*0.5*(gas_gamma-1.))/(c_L / P_L**beta + c_R / P_R**beta))**(1./beta)
    P_3 = 0.5*(P_R + P_L)
    f_L = 1.
    while fabs(P_3 - P_new) > 1e-6:
        P_3 = P_new
        (f_L, dfdp_L) = compute_f(P_3, P_L, c_L)
        (f_R, dfdp_R) = compute_f(P_3, P_R, c_R)
        f = f_L + f_R + (v_R - v_L)
        df = dfdp_L + dfdp_R
        dp =  -f/df
        prnew = P_3 + dp
    v_3 = v_L - f_L
    return (P_new, v_3)


# Solve Riemann problem for post-shock region
(P_3, v_3) = RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R)

# Check direction of shocks and wave
shock_R = (P_3 > P_R)
shock_L = (P_3 > P_L)

# Velocity of shock front and and rarefaction wave
if shock_R:
    v_right = v_R + c_R**2*(P_3/P_R - 1.)/(gas_gamma*(v_3-v_R))
else:
    v_right = c_R + 0.5*(gas_gamma+1.)*v_3 - 0.5*(gas_gamma-1.)*v_R

if shock_L:
    v_left = v_L + c_L**2*(P_3/p_L - 1.)/(gas_gamma*(v_3-v_L))
else:
    v_left = c_L - 0.5*(gas_gamma+1.)*v_3 + 0.5*(gas_gamma-1.)*v_L

# Compute position of the transitions
x_23 = -fabs(v_left) * time
if shock_L :
    x_12 = -fabs(v_left) * time
else:
    x_12 = -(c_L - v_L) * time

x_34 = v_3 * time

x_45 = fabs(v_right) * time
if shock_R:
    x_56 = fabs(v_right) * time
else:
    x_56 = (c_R + v_R) * time


# Prepare arrays
delta_x = (x_max - x_min) / N
x_s = arange(x_min, x_max, delta_x)
rho_s = zeros(N)
P_s = zeros(N)
v_s = zeros(N)

# Compute solution in the different regions
for i in range(N):
    if x_s[i] <= x_12:
        rho_s[i] = rho_L
        P_s[i] = P_L
        v_s[i] = v_L
    if x_s[i] >= x_12 and x_s[i] < x_23:
        if shock_L:
            rho_s[i] = rho_L*(Gama + P_3/P_L)/(1. + Gama * P_3/P_L)
            P_s[i] = P_3
            v_s[i] = v_3
        else:
            rho_s[i] = rho_L*(Gama * (0. - x_s[i])/(c_L * time) + Gama * v_L/c_L + (1.-Gama))**(2./(gas_gamma-1.))
            P_s[i] = P_L*(rho_s[i] / rho_L)**gas_gamma
            v_s[i] = (1.-Gama)*(c_L -(0. - x_s[i]) / time) + Gama*v_L
    if x_s[i] >= x_23 and x_s[i] < x_34:
        if shock_L:
            rho_s[i] = rho_L*(Gama + P_3/P_L)/(1+Gama * P_3/p_L)
        else:
            rho_s[i] = rho_L*(P_3 / P_L)**(1./gas_gamma)
        P_s[i] = P_3
        v_s[i] = v_3
    if x_s[i] >= x_34 and x_s[i] < x_45:
        if shock_R:
            rho_s[i] = rho_R*(Gama + P_3/P_R)/(1. + Gama * P_3/P_R)
        else:
            rho_s[i] = rho_R*(P_3 / P_R)**(1./gas_gamma)
        P_s[i] = P_3
        v_s[i] = v_3
    if x_s[i] >= x_45 and x_s[i] < x_56:
        if shock_R:
            rho_s[i] = rho_R
            P_s[i] = P_R
            v_s[i] = v_R
        else:
            rho_s[i] = rho_R*(Gama*(x_s[i])/(c_R*time) - Gama*v_R/c_R + (1.-Gama))**(2./(gas_gamma-1.))
            P_s[i] = p_R*(rho_s[i]/rho_R)**gas_gamma
            v_s[i] = (1.-Gama)*(-c_R - (-x_s[i])/time) + Gama*v_R
    if x_s[i] >= x_56:
        rho_s[i] = rho_R
        P_s[i] = P_R
        v_s[i] = v_R


# Additional arrays
u_s = P_s / (rho_s * (gas_gamma - 1.))  #internal energy
s_s = P_s / rho_s**gas_gamma # entropic function
        

# Plot the interesting quantities
figure()

# Velocity profile --------------------------------
subplot(231)
228
plot(x, v, '.', color='r', ms=4.0)
229
plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
230
231
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
232
233
234
235
236
xlim(-0.5, 0.5)
ylim(-0.1, 0.95)

# Density profile --------------------------------
subplot(232)
237
plot(x, rho, '.', color='r', ms=4.0)
238
plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
239
240
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
241
242
243
244
245
xlim(-0.5, 0.5)
ylim(0.05, 1.1)

# Pressure profile --------------------------------
subplot(233)
246
plot(x, P, '.', color='r', ms=4.0)
247
plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
248
249
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Pressure}}~P$", labelpad=0)
250
251
252
253
254
xlim(-0.5, 0.5)
ylim(0.01, 1.1)

# Internal energy profile -------------------------
subplot(234)
255
plot(x, u, '.', color='r', ms=4.0)
256
plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
257
258
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
259
260
261
262
263
xlim(-0.5, 0.5)
ylim(0.8, 2.2)

# Entropy profile ---------------------------------
subplot(235)
264
plot(x, S, '.', color='r', ms=4.0)
265
plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
266
267
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Entropy}}~S$", labelpad=0)
268
269
270
271
272
273
274
275
276
277
xlim(-0.5, 0.5)
ylim(0.8, 3.8)

# Information -------------------------------------
subplot(236, frameon=False)

text(-0.49, 0.9, "Sod shock with  $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
text(-0.49, 0.8, "Left:~~ $(P_L, \\rho_L, v_L) = (%.3f, %.3f, %.3f)$"%(P_L, rho_L, v_L), fontsize=10)
text(-0.49, 0.7, "Right: $(P_R, \\rho_R, v_R) = (%.3f, %.3f, %.3f)$"%(P_R, rho_R, v_R), fontsize=10)
plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
278
279
280
281
text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
text(-0.49, 0.4, scheme, fontsize=10)
text(-0.49, 0.3, kernel, fontsize=10)
text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
282
283
284
285
286
287
288
xlim(-0.5, 0.5)
ylim(0, 1)
xticks([])
yticks([])


savefig("SodShock.png", dpi=200)