mulimoen
/
SummationByParts
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

eulerplot slider

This commit is contained in:
Magnus Ulimoen 2020-04-08 19:25:44 +02:00
parent b2eb79e71d
commit e395c9c740
1 changed files with 124 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

143
sbp/eulerplot
View File

@ -25,38 +25,42 @@ def plot_all(grids, save: bool, filename="figure.png"):
f, axarr = plt.subplots(2, 2) f, axarr = plt.subplots(2, 2)
min_rho = min(np.min(g["rho"]) for g in grids) min_rho = min(np.min(g["rho"][-1, :, :]) for g in grids)
max_rho = max(np.max(g["rho"]) for g in grids) max_rho = max(np.max(g["rho"][-1, :, :]) for g in grids)
r = 1.2 * max(abs(min_rho - 1), abs(max_rho - 1)) r = 1.2 * max(abs(min_rho - 1), abs(max_rho - 1))
rho_levels = np.linspace(1 - r, 1 + r, 34) rho_levels = np.linspace(1 - r, 1 + r, 34)
min_rhou = min(np.min(g["rhou"]) for g in grids) min_rhou = min(np.min(g["rhou"][-1, :, :]) for g in grids)
max_rhou = max(np.max(g["rhov"]) for g in grids) max_rhou = max(np.max(g["rhov"][-1, :, :]) for g in grids)
r = 1.2 * max(abs(min_rhou - 1), abs(max_rhou - 1)) r = 1.2 * max(abs(min_rhou - 1), abs(max_rhou - 1))
rhou_levels = np.linspace(1 - r, 1 + r, 20) rhou_levels = np.linspace(1 - r, 1 + r, 20)
min_rhov = min(np.min(g["rhov"]) for g in grids) min_rhov = min(np.min(g["rhov"][-1, :, :]) for g in grids)
max_rhov = max(np.max(g["rhov"]) for g in grids) max_rhov = max(np.max(g["rhov"][-1, :, :]) for g in grids)
r = 1.2 * max(abs(min_rhov), abs(max_rhov)) r = 1.2 * max(abs(min_rhov), abs(max_rhov))
rhov_levels = np.linspace(-r, r, 20) rhov_levels = np.linspace(-r, r, 20)
min_e = min(np.min(g["e"]) for g in grids) min_e = min(np.min(g["e"][-1, :, :]) for g in grids)
max_e = max(np.max(g["e"]) for g in grids) max_e = max(np.max(g["e"][-1, :, :]) for g in grids)
e_levels = np.linspace(min_e, max_e) e_levels = np.linspace(min_e, max_e)
for g in grids: for g in grids:
x = g["x"] x = g["x"]
y = g["y"] y = g["y"]
axarr[0, 0].contourf(x, y, g["rho"], cmap=sym_cmap, levels=rho_levels) axarr[0, 0].contourf(x, y, g["rho"][-1, :, :], cmap=sym_cmap, levels=rho_levels)
gridlines(axarr[0, 0], x, y) gridlines(axarr[0, 0], x, y)
axarr[0, 1].contourf(x, y, g["rhou"], cmap=sym_cmap, levels=rhou_levels) axarr[0, 1].contourf(
x, y, g["rhou"][-1, :, :], cmap=sym_cmap, levels=rhou_levels
)
gridlines(axarr[0, 1], x, y) gridlines(axarr[0, 1], x, y)
axarr[1, 0].contourf(x, y, g["rhov"], cmap=sym_cmap, levels=rhov_levels) axarr[1, 0].contourf(
x, y, g["rhov"][-1, :, :], cmap=sym_cmap, levels=rhov_levels
)
gridlines(axarr[1, 0], x, y) gridlines(axarr[1, 0], x, y)
axarr[1, 1].contourf(x, y, g["e"], cmap=e_cmap, levels=e_levels) axarr[1, 1].contourf(x, y, g["e"][-1, :, :], cmap=e_cmap, levels=e_levels)
gridlines(axarr[1, 1], x, y) gridlines(axarr[1, 1], x, y)
axarr[0, 0].set_title(r"$\rho$") axarr[0, 0].set_title(r"$\rho$")
@ -97,13 +101,23 @@ def plot_all(grids, save: bool, filename="figure.png"):
plt.show() plt.show()
def pressure(rho, rhou, rhov, e):
gamma = 1.4
return (gamma - 1) * (e - (rhou ** 2 + rhov ** 2) / (2 * rho))
def plot_pressure(grids, save: bool, filename="figure.png"): def plot_pressure(grids, save: bool, filename="figure.png"):
cmap = plt.get_cmap("RdGy") cmap = plt.get_cmap("RdGy")
gamma = 1.4 # Assumption might be wrong
Mach = 0.5 Mach = 0.5
gamma = 1.4
p = [ p = [
(gamma - 1) * (g["e"] - (g["rhou"] ** 2 + g["rhov"] ** 2) / (2 * g["rho"])) pressure(
g["rho"][-1, :, :],
g["rhou"][-1, :, :],
g["rhov"][-1, :, :],
g["e"][-1, :, :],
)
for g in grids for g in grids
] ]
@ -142,6 +156,93 @@ def plot_pressure(grids, save: bool, filename="figure.png"):
plt.show() plt.show()
def plot_pressure_slider(grids, save: bool, filename="figure.png"):
cmap = plt.get_cmap("RdGy")
gamma = 1.4 # Assumption might be wrong
Mach = 0.5
def p(itime):
return [
pressure(
g["rho"][itime, :, :],
g["rhou"][itime, :, :],
g["rhov"][itime, :, :],
g["e"][itime, :, :],
)
for g in grids
]
max_p = 3.0
min_p = 1.75
p_inf = 1 / (gamma * Mach ** 2)
r = max(max_p - p_inf, p_inf - min_p)
levels = np.linspace(p_inf - r, p_inf + r, 30)
fig = plt.figure()
gs = mpl.gridspec.GridSpec(
2, 2, figure=fig, width_ratios=[1, 0.02], height_ratios=[1, 0.02]
)
ax = fig.add_subplot(gs[0, 0])
slider_ax = fig.add_subplot(gs[1, 0])
cbar_ax = fig.add_subplot(gs[0, 1])
xmin, xmax = np.inf, -np.inf
ymin, ymax = np.inf, -np.inf
for g in grids:
x = g["x"]
xmin = min(xmin, x.min())
xmax = max(xmax, x.max())
y = g["y"]
ymin = min(ymin, y.min())
ymax = max(ymax, y.max())
gridlines(ax, x, y)
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
plt.title("Pressure")
norm = mpl.colors.Normalize(vmin=levels[0], vmax=levels[-1])
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
plt.colorbar(sm, cax=cbar_ax)
plt.xlabel("x")
plt.ylabel("y")
itime = len(t) - 1
slider = mpl.widgets.Slider(
slider_ax, "itime", 0, itime, valinit=itime, valstep=1, valfmt="%0.0f"
)
contours = []
def update(itime):
global contours
contours = []
itime = int(itime)
for ct in contours:
for coll in ct.collections:
coll.remove()
for g in grids:
ct = ax.contourf(
g["x"],
g["y"],
pressure(
g["rho"][itime, :, :],
g["rhou"][itime, :, :],
g["rhov"][itime, :, :],
g["e"][itime, :, :],
),
cmap=cmap,
levels=levels,
)
contours.append(ct)
slider.valtext.set_text(t[itime])
update(itime)
slider.on_changed(update)
plt.show()
def read_from_file(filename): def read_from_file(filename):
grids = [] grids = []
@ -155,10 +256,10 @@ def read_from_file(filename):
{ {
"x": group["x"][:], "x": group["x"][:],
"y": group["y"][:], "y": group["y"][:],
"rho": group["rho"][-1, :, :], "rho": group["rho"][:],
"rhou": group["rhou"][-1, :, :], "rhou": group["rhou"][:],
"rhov": group["rhov"][-1, :, :], "rhov": group["rhov"][:],
"e": group["e"][-1, :, :], "e": group["e"][:],
} }
) )
@ -179,6 +280,7 @@ if __name__ == "__main__":
parser.add_argument( parser.add_argument(
"-a", help="Show all four variables", action="store_true", dest="all" "-a", help="Show all four variables", action="store_true", dest="all"
) )
parser.add_argument("--slider", help="Add slider", action="store_true")
args = parser.parse_args() args = parser.parse_args()
filename = args.filename filename = args.filename
@ -188,4 +290,7 @@ if __name__ == "__main__":
if args.all: if args.all:
plot_all(grids, args.save, args.output) plot_all(grids, args.save, args.output)
else: else:
plot_pressure(grids, args.save, args.output) if args.slider:
plot_pressure_slider(grids, t)
else:
plot_pressure(grids, args.save, args.output)