# Exact Riemann solver, based on Toro, Suresh # Philip Mocz, 2016 import matplotlib.pyplot as plt import numpy as np from matplotlib.widgets import Slider, Button """ Create Your Own Riemann Solver (With Python) Philip Mocz (2023), @PMocz Solve the Riemann problem for the Euler equations. Interactive! Slightly modified by M.Sato. added plot for internal energy, simplified pressure guess, moved all slider bar at the right side of canvas, etc Reference: [1]: E.F. Toro. Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer, Berlin, New York, 2009. """ class RiemannSolver: def __init__(self, rho_L, vx_L, P_L, rho_R, vx_R, P_R, gamma, t): """ Initialize Riemann Solver """ self.gamma = gamma # Gamma-related constants self.g1 = (gamma - 1.0) / (2.0 * gamma) self.g2 = (gamma + 1.0) / (2.0 * gamma) self.g3 = 2.0 * gamma / (gamma - 1.0) self.g4 = 2.0 / (gamma - 1.0) self.g5 = 2.0 / (gamma + 1.0) self.g6 = (gamma - 1.0) / (gamma + 1.0) self.g7 = (gamma - 1.0) / 2.0 self.g8 = gamma - 1.0 self.g9 = 1.0 / gamma self.set_state(rho_L, vx_L, P_L, rho_R, vx_R, P_R, t) Lbox = 1.0 # domain [-Lbox/2,Lbox/2], shock initially at x=0 N = 1000 # resolution self.x = np.linspace(-Lbox / 2, Lbox / 2, N) self.rho = np.zeros(N) self.vx = np.zeros(N) self.P = np.zeros(N) self.ene = np.zeros(N) def set_state(self, rho_L, vx_L, P_L, rho_R, vx_R, P_R, t): """ Set Left and Right States """ self.rho_L = rho_L self.vx_L = vx_L self.P_L = P_L self.rho_R = rho_R self.vx_R = vx_R self.P_R = P_R self.t = t self.P_star = 0 # pressure solution in star region self.vx_star = 0 # velocity solution in star region self.success = False # solve successful? # sound speeds self.c_L = np.sqrt(self.gamma * P_L / rho_L) self.c_R = np.sqrt(self.gamma * P_R / rho_R) def calc_star_P_and_vx(self): """ Compute solution for pressure & velocity in the star region """ tolerance = 1.0e-8 max_iter = 100 P_old = self.guess_P() vx_diff = self.vx_R - self.vx_L # compute pressure in star region via Newton-Raphson iteration for i in np.arange(max_iter): (f_L, df_L) = self.pressure_function_and_derivative( P_old, self.rho_L, self.P_L, self.c_L ) (f_R, df_R) = self.pressure_function_and_derivative( P_old, self.rho_R, self.P_R, self.c_R ) P = P_old - (f_L + f_R + vx_diff) / (df_L + df_R) change = 2.0 * abs((P - P_old) / (P + P_old)) if change < tolerance: break if P < 0.0: P = tolerance P_old = P # compute velocity in Star Region vx = 0.5 * (self.vx_L + self.vx_R + f_R - f_L) return (P, vx) def guess_P(self): # primitive variable riemann solver cup = 0.25 * (self.rho_L + self.rho_R) * (self.c_L + self.c_R) P_pv = 0.5 * (self.P_L + self.P_R) + 0.5 * (self.vx_L - self.vx_R) * cup P_guess = max(1e-6, P_pv) return P_guess def pressure_function_and_derivative(self, P, rho_k, P_k, c_k): """ Evaluate functions to solve for pressure in Newton-Raphson iterator """ if P <= P_k: # Rarefaction Wave q = P / P_k f = self.g4 * c_k * (q**self.g1 - 1.0) df = (1.0 / (rho_k * c_k)) * q ** (-self.g2) else: # Shock Wave ak = self.g5 / rho_k bk = self.g6 * P_k qrt = np.sqrt(ak / (bk + P)) f = (P - P_k) * qrt df = (1.0 - 0.5 * (P - P_k) / (bk + P)) * qrt return (f, df) def sample(self, P_star, vx_star, s): """ Sample the solution, given the Star region pressure and velocity, in terms of s = x/t """ if s <= vx_star: # Sampling point lies to the left of the contact discontinuity if P_star <= self.P_L: # Left Rarefaction sh_L = self.vx_L - self.c_L if s <= sh_L: # Sampled point is left data state rho = self.rho_L vx = self.vx_L P = self.P_L else: cmL = self.c_L * (P_star / self.P_L) ** self.g1 st_L = vx_star - cmL if s > st_L: # Sampled point is Star Left state rho = self.rho_L * (P_star / self.P_L) ** self.g9 vx = vx_star P = P_star else: # Sampled point is inside left fan vx = self.g5 * (self.c_L + self.g7 * self.vx_L + s) c = self.g5 * (self.c_L + self.g7 * (self.vx_L - s)) rho = self.rho_L * (c / self.c_L) ** self.g4 P = self.P_L * (c / self.c_L) ** self.g3 else: # Left shock P_starL = P_star / self.P_L s_L = self.vx_L - self.c_L * np.sqrt(self.g2 * P_starL + self.g1) if s <= s_L: # Sampled point is left data state rho = self.rho_L vx = self.vx_L P = self.P_L else: # Sampled point is Star Left state rho = self.rho_L * (P_starL + self.g6) / (P_starL * self.g6 + 1.0) vx = vx_star P = P_star else: # Sampling point lies to the right of the contact discontinuity if P_star > self.P_R: # Right Shock P_starR = P_star / self.P_R s_R = self.vx_R + self.c_R * np.sqrt(self.g2 * P_starR + self.g1) if s >= s_R: # Sampled point is right data state rho = self.rho_R vx = self.vx_R P = self.P_R else: # Sampled point is Star Right state rho = self.rho_R * (P_starR + self.g6) / (P_starR * self.g6 + 1.0) vx = vx_star P = P_star else: # Right Rarefaction sh_R = self.vx_R + self.c_R if s >= sh_R: # Sampled point is right data state rho = self.rho_R vx = self.vx_R P = self.P_R else: cmR = self.c_R * (P_star / self.P_R) ** self.g1 st_R = vx_star + cmR if s <= st_R: # Sampled point is Star Right state rho = self.rho_R * (P_star / self.P_R) ** self.g9 vx = vx_star P = P_star else: # Sampled point is inside left fan vx = self.g5 * (-self.c_R + self.g7 * self.vx_R + s) c = self.g5 * (self.c_R - self.g7 * (self.vx_R - s)) rho = self.rho_R * (c / self.c_R) ** self.g4 P = self.P_R * (c / self.c_R) ** self.g3 return (rho, vx, P) def solve(self): # Check pressure positivity condition if self.g4 * (self.c_L + self.c_R) < (self.vx_R - self.vx_L): print("Error: initial data is such that the vacuum is generated!") self.success = False # Find exact solution for pressure & velocity in star region (P_star, vx_star) = self.calc_star_P_and_vx() # print(f"P_star={P_star} vx_star={vx_star}") for i in np.arange(len(self.x)): s = self.x[i] / self.t (rho, vx, P) = self.sample(P_star, vx_star, s) self.rho[i] = rho self.vx[i] = vx self.P[i] = P self.ene[i] = P/(rho*self.g8) self.success = True return (self.x, self.rho, self.vx, self.P, self.ene) def main(): """Interactive Riemann Solver""" #----------------------------- test-1 # Left State rho_L = 1.0 vx_L = 0.0 P_L = 1.0 # Right State rho_R = 0.125 vx_R = 0.0 P_R = 0.1 # time t = 0.25 title = "test1" #----------------------------- test-2 # # Left State # rho_L = 1.0 # vx_L = -2.0 # P_L = 0.4 # # # Right State # rho_R = 1.0 # vx_R = 2.0 # P_R = 0.4 # # # time # t = 0.15 # title = "test2" #----------------------------- test-3 # # Left State # rho_L = 1.0 # vx_L = 0.0 # P_L = 1000.0 # # # Right State # rho_R = 1.0 # vx_R = 0.0 # P_R = 0.01 # # # time # t = 0.012 # title = "test3" #----------------------------- test-4 # # Left State # rho_L = 1.0 # vx_L = 0.0 # P_L = 0.01 # # # Right State # rho_R = 1.0 # vx_R = 0.0 # P_R = 100.0 # # # time # t = 0.035 # title = "test4" #----------------------------- test-5 # # Left State # rho_L = 5.99924 # vx_L = 19.5975 # P_L = 460.894 # # # Right State # rho_R = 5.99242 # vx_R = -6.19633 # P_R = 46.0950 # # # time # t = 0.035 # title = "test5" # ideal gas gamma gamma = 1.4 # Riemann Solver rs = RiemannSolver(rho_L, vx_L, P_L, rho_R, vx_R, P_R, gamma, t) x, rho, vx, P, ene = rs.solve() # --- 1. 画面サイズとレイアウトの調整 --- fig = plt.figure(figsize=(8, 6), dpi=150) plt.subplots_adjust(left=0.1, right=0.50, top=0.9, bottom=0.1, hspace=0.3) ax1 = fig.add_subplot(411) (l1,) = plt.plot(x, rho_L * (x <= 0) + rho_R * (x > 0), "--") (l1b,) = plt.plot(x, rho, linewidth=2) ax1.set_xticks([]) ax1.set_xlim(-0.5, 0.5) ax1.set_ylabel("rho") ax1.set_ylim(0, 1.1 * np.max(rho)) ax2 = fig.add_subplot(412) (l2,) = plt.plot(x, vx_L * (x <= 0) + vx_R * (x > 0), "--") (l2b,) = plt.plot(x, vx, linewidth=2) ax2.set_xticks([]) ax2.set_xlim(-0.5, 0.5) ax2.set_ylabel("vx") ax2.set_ylim(min(0, 1.1 * np.min(vx)), max(0, 1.1 * np.max(vx))) ax3 = fig.add_subplot(413) (l3,) = plt.plot(x, P_L * (x <= 0) + P_R * (x > 0), "--") (l3b,) = plt.plot(x, P, linewidth=2) ax3.set_xticks([]) ax3.set_xlim(-0.5, 0.5) ax3.set_ylabel("P") ax3.set_ylim(0, 1.1 * np.max(P)) ax4 = fig.add_subplot(414) (l4,) = plt.plot(x, (P_L/(gamma-1)/rho_L) * (x <= 0) + (P_R/(gamma-1)/rho_R) * (x > 0), "--") (l4b,) = plt.plot(x, ene, linewidth=2) ax4.set_xlim(-0.5, 0.5) ax4.set_ylabel("e") ax4.set_ylim(0, 1.1 * np.max(ene)) # --- 2. スライダーの座標設定 --- # [左端のX座標, 下端のY座標, 横幅, 高さ] s_x = 0.6 s_w = 0.25 s_h = 0.025 ax_t = plt.axes([s_x, 0.80, s_w, s_h]) ax_rho_L = plt.axes([s_x, 0.70, s_w, s_h]) ax_vx_L = plt.axes([s_x, 0.65, s_w, s_h]) ax_P_L = plt.axes([s_x, 0.60, s_w, s_h]) ax_rho_R = plt.axes([s_x, 0.50, s_w, s_h]) ax_vx_R = plt.axes([s_x, 0.45, s_w, s_h]) ax_P_R = plt.axes([s_x, 0.40, s_w, s_h]) # --- 3. スライダーの作成（可動範囲を初期値に合わせて調整） --- sl_t = Slider(ax_t, "time", 0.0, 0.5, valinit=t) sl_rho_L = Slider(ax_rho_L, "rho_L", 0.001, rho_L*2.0, valinit=rho_L) sl_vx_L = Slider(ax_vx_L, "vx_L", vx_L-10.0, vx_L+10.0, valinit=vx_L) sl_P_L = Slider(ax_P_L, "P_L", 0.001, P_L*2.0, valinit=P_L) sl_rho_R = Slider(ax_rho_R, "rho_R", 0.001, rho_R*2.0, valinit=rho_R) sl_vx_R = Slider(ax_vx_R, "vx_R", vx_R-10.0, vx_R+10.0, valinit=vx_R) sl_P_R = Slider(ax_P_R, "P_R", 0.001, P_R*2.0, valinit=P_R) # update plot function def update(val): t = sl_t.val rho_L = sl_rho_L.val vx_L = sl_vx_L.val P_L = sl_P_L.val rho_R = sl_rho_R.val vx_R = sl_vx_R.val P_R = sl_P_R.val rs.set_state(rho_L, vx_L, P_L, rho_R, vx_R, P_R, t) x, rho, vx, P, ene = rs.solve() l1.set_ydata(rho_L * (x <= 0) + rho_R * (x > 0)) l2.set_ydata(vx_L * (x <= 0) + vx_R * (x > 0)) l3.set_ydata(P_L * (x <= 0) + P_R * (x > 0)) l4.set_ydata(P_L/(rs.g8*rho_L) * (x <= 0) + P_R/(rs.g8*rho_R) * (x > 0)) l1b.set_ydata(rho) l2b.set_ydata(vx) l3b.set_ydata(P) l4b.set_ydata(ene) ax1.set_ylim(0, 1.1 * np.max(rho)) ax2.set_ylim(min(0, 1.1 * np.min(vx)), max(0, 1.1 * np.max(vx))) ax3.set_ylim(0, 1.1 * np.max(P)) ax4.set_ylim(0, 1.1 * np.max(ene)) fig.canvas.draw_idle() # Call update function when a slider is changed sl_t.on_changed(update) sl_rho_L.on_changed(update) sl_vx_L.on_changed(update) sl_P_L.on_changed(update) sl_rho_R.on_changed(update) sl_vx_R.on_changed(update) sl_P_R.on_changed(update) # Save figure plt.savefig(f"{title}.png", dpi=240) plt.show() return 0 if __name__ == "__main__": main()