From d1e355aeb6c4f5fd417ed88f148fe2741b0e9602 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sat, 14 Feb 2026 20:11:06 +0100 Subject: [PATCH 1/7] [Doc] add sod tube example to sphinx --- doc/sphinx/examples/ramses/run_sod.py | 246 ++++++++++++++++++++++++++ src/shamphys/src/SodTube.cpp | 3 +- 2 files changed, 248 insertions(+), 1 deletion(-) create mode 100644 doc/sphinx/examples/ramses/run_sod.py diff --git a/doc/sphinx/examples/ramses/run_sod.py b/doc/sphinx/examples/ramses/run_sod.py new file mode 100644 index 0000000000..7422bc8d07 --- /dev/null +++ b/doc/sphinx/examples/ramses/run_sod.py @@ -0,0 +1,246 @@ +""" +Advection test in RAMSES solver +============================================= + +Compare advection with all slope limiters & Riemann solvers +""" + +import os + +import matplotlib.pyplot as plt +import numpy as np + +import shamrock + +tmax = 0.245 +timestamps = 40 +gamma = 1.4 + +multx = 4 +multy = 1 +multz = 1 + +sz = 1 << 1 +base = 16 + +positions = [(x, 0, 0) for x in np.linspace(0.5, 1.5, 256).tolist()[:-1]] + + +def run_advect(slope_limiter: str, riemann_solver: str, only_last_step: bool = True): + ctx = shamrock.Context() + ctx.pdata_layout_new() + + model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3") + + cfg = model.gen_default_config() + scale_fact = 2 / (sz * base * multx) + cfg.set_scale_factor(scale_fact) + cfg.set_eos_gamma(gamma) + + if slope_limiter == "none": + cfg.set_slope_lim_none() + elif slope_limiter == "vanleer": + cfg.set_slope_lim_vanleer_f() + elif slope_limiter == "vanleer_std": + cfg.set_slope_lim_vanleer_std() + elif slope_limiter == "vanleer_sym": + cfg.set_slope_lim_vanleer_sym() + elif slope_limiter == "minmod": + cfg.set_slope_lim_minmod() + else: + raise ValueError(f"Invalid slope limiter: {slope_limiter}") + + if riemann_solver == "rusanov": + cfg.set_riemann_solver_rusanov() + elif riemann_solver == "hll": + cfg.set_riemann_solver_hll() + elif riemann_solver == "hllc": + cfg.set_riemann_solver_hllc() + else: + raise ValueError(f"Invalid Riemann solver: {riemann_solver}") + + cfg.set_face_time_interpolation(True) + model.set_solver_config(cfg) + + model.init_scheduler(int(1e7), 1) + model.make_base_grid((0, 0, 0), (sz, sz, sz), (base * multx, base * multy, base * multz)) + + + def rho_map(rmin, rmax): + x, y, z = rmin + if x < 1: + return 1 + else: + return 0.125 + + + etot_L = 1.0 / (gamma - 1) + etot_R = 0.1 / (gamma - 1) + + + def rhoetot_map(rmin, rmax): + rho = rho_map(rmin, rmax) + + x, y, z = rmin + if x < 1: + return etot_L + else: + return etot_R + + + def rhovel_map(rmin, rmax): + rho = rho_map(rmin, rmax) + + return (0, 0, 0) + + model.set_field_value_lambda_f64("rho", rho_map) + model.set_field_value_lambda_f64("rhoetot", rhoetot_map) + model.set_field_value_lambda_f64_3("rhovel", rhovel_map) + + results = [] + + def analysis(iplot: int): + rho_vals = model.render_slice("rho", "f64", positions) + rhov_vals = model.render_slice("rhovel", "f64_3", positions) + rhoetot_vals = model.render_slice("rhoetot", "f64", positions) + + + vx = np.array(rhov_vals)[:,0] / np.array(rho_vals) + P = (np.array(rhoetot_vals) - 0.5 * np.array(rho_vals) * vx**2) * (gamma - 1) + results_dic = { + "rho": np.array(rho_vals), + "vx": np.array(vx), + "P": np.array(P), + } + results.append(results_dic) + + if only_last_step: + model.evolve_until(tmax) + analysis(timestamps) + else: + dt_evolve = tmax / timestamps + + for i in range(timestamps + 1): + model.evolve_until(dt_evolve * i) + analysis(i) + + return results + + +# %% +data = {} +data["none_rusanov"] = run_advect("none", "rusanov") +data["none_hll"] = run_advect("none", "hll") +data["none_hllc"] = run_advect("none", "hllc") +data["minmod_rusanov"] = run_advect("minmod", "rusanov") +data["minmod_hll"] = run_advect("minmod", "hll") +data["minmod_hllc"] = run_advect("minmod", "hllc", only_last_step=False) + +# %% +# Plot 1: Comparison against analytical solution +riemann_solvers = ["rusanov", "hll", "hllc"] +slope_limiters = ["none", "minmod"] + +xref = 1.0 +xrange = 0.5 +sod = shamrock.phys.SodTube(gamma=gamma, rho_1=1, P_1=1, rho_5=0.125, P_5=0.1) + +#### add analytical soluce +arr_x = [x[0] for x in positions] + +arr_rho = [] +arr_P = [] +arr_vx = [] + +for i in range(len(arr_x)): + x_ = arr_x[i] - xref + + _rho, _vx, _P = sod.get_value(tmax, x_) + #print(x_,_rho, _vx, _P) + arr_rho.append(_rho) + arr_vx.append(_vx) + arr_P.append(_P) + +arr_rho = np.array(arr_rho) +arr_vx = np.array(arr_vx) +arr_P = np.array(arr_P) + +fig, axes = plt.subplots(3, 1, figsize=(6, 15)) +fig.suptitle(f"t={tmax} (Last Step)", fontsize=14) + +for i in range(3): + axes[i].set_xlabel("$x$") + axes[i].set_yscale("log") + axes[i].grid(True, alpha=0.3) + +ax1, ax2, ax3 = axes +ax1.set_xlabel("$x$") +ax1.set_ylabel("$\\rho$") + +ax2.set_xlabel("$x$") +ax2.set_ylabel("$v_x$") + +ax3.set_xlabel("$x$") +ax3.set_ylabel("$P$") + + +# ax1.plot(arr_x, arr_rho, color="black", label="analytic") +# ax2.plot(arr_x, arr_vx, color="black", label="analytic") +# ax3.plot(arr_x, arr_P, color="black", label="analytic") + + +for limiter in slope_limiters: + for solver in riemann_solvers: + key = f"{limiter}_{solver}" + if key in data: + # Get the last timestep + delta_rho = np.abs(data[key][-1]["rho"] - arr_rho) + delta_vx = np.abs(data[key][-1]["vx"] - arr_vx) + delta_P = np.abs(data[key][-1]["P"] - arr_P) + + ax1.plot(arr_x, delta_rho, label=f"{limiter} {solver} (rho)", linewidth=1) + ax2.plot(arr_x, delta_vx, label=f"{limiter} {solver} (vx)", linewidth=1) + ax3.plot(arr_x, delta_P, label=f"{limiter} {solver} (P)", linewidth=1) + + +ax1.legend() +ax2.legend() +ax3.legend() + +plt.tight_layout() +plt.show() + +# %% +# Plot 2: Animation of vanleer_sym_hllc configuration + +# sphinx_gallery_thumbnail_number = 2 + +from matplotlib.animation import FuncAnimation + +fig2, ax2 = plt.subplots() +ax2.set_xlabel("$x$") +ax2.set_ylabel("$\\rho$") +ax2.set_ylim(-0.1, 1.1) +ax2.set_xlim(0.5, 1.5) +ax2.grid(True, alpha=0.3) + +(line_rho,) = ax2.plot(arr_x, data["minmod_hllc"][0]["rho"], label="rho", linewidth=2) +(line_vx,) = ax2.plot(arr_x, data["minmod_hllc"][0]["vx"], label="vx", linewidth=2) +(line_P,) = ax2.plot(arr_x, data["minmod_hllc"][0]["P"], label="P", linewidth=2) + +ax2.legend() +ax2.set_title(f"minmod_hllc - t = {0.0:.3f} s") + + +def animate(frame): + t = tmax * frame / timestamps + line_rho.set_ydata(data["minmod_hllc"][frame]["rho"]) + line_vx.set_ydata(data["minmod_hllc"][frame]["vx"]) + line_P.set_ydata(data["minmod_hllc"][frame]["P"]) + ax2.set_title(f"minmod_hllc - t = {t:.3f} s") + return (line_rho, line_vx, line_P) + + +anim = FuncAnimation(fig2, animate, frames=timestamps + 1, interval=50, blit=False, repeat=True) +plt.tight_layout() +plt.show() diff --git a/src/shamphys/src/SodTube.cpp b/src/shamphys/src/SodTube.cpp index 8bd24758d2..9977ca2d9f 100644 --- a/src/shamphys/src/SodTube.cpp +++ b/src/shamphys/src/SodTube.cpp @@ -14,6 +14,7 @@ */ #include "shamphys/SodTube.hpp" +#include "shamcomm/logs.hpp" #include "shammath/derivatives.hpp" #include "shammath/solve.hpp" #include @@ -22,7 +23,7 @@ f64 shamphys::SodTube::soundspeed(f64 P, f64 rho) { return std::sqrt(gamma * P / rho); }; shamphys::SodTube::SodTube(f64 _gamma, f64 _rho_1, f64 _P_1, f64 _rho_5, f64 _P_5) - : gamma(_gamma), rho_1(_rho_1), P_1(_P_1), rho_5(_rho_5), P_5(_P_5) { + : gamma(_gamma), rho_1(_rho_1), P_1(_P_1), rho_5(_rho_5), P_5(_P_5), vx_1(0), vx_5(0) { c_1 = soundspeed(P_1, rho_1); c_5 = soundspeed(P_5, rho_5); From 4755e54fff50fd021cc660a7d19132face96ec08 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sat, 14 Feb 2026 20:11:29 +0100 Subject: [PATCH 2/7] [Doc] add sod tube example to sphinx --- doc/sphinx/examples/ramses/run_sod.py | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/doc/sphinx/examples/ramses/run_sod.py b/doc/sphinx/examples/ramses/run_sod.py index 7422bc8d07..b5299aef66 100644 --- a/doc/sphinx/examples/ramses/run_sod.py +++ b/doc/sphinx/examples/ramses/run_sod.py @@ -1,8 +1,8 @@ """ -Advection test in RAMSES solver +Sod tube test in RAMSES solver ============================================= -Compare advection with all slope limiters & Riemann solvers +Compare Sod tube with all slope limiters & Riemann solvers """ import os @@ -65,7 +65,6 @@ def run_advect(slope_limiter: str, riemann_solver: str, only_last_step: bool = T model.init_scheduler(int(1e7), 1) model.make_base_grid((0, 0, 0), (sz, sz, sz), (base * multx, base * multy, base * multz)) - def rho_map(rmin, rmax): x, y, z = rmin if x < 1: @@ -73,11 +72,9 @@ def rho_map(rmin, rmax): else: return 0.125 - etot_L = 1.0 / (gamma - 1) etot_R = 0.1 / (gamma - 1) - def rhoetot_map(rmin, rmax): rho = rho_map(rmin, rmax) @@ -87,7 +84,6 @@ def rhoetot_map(rmin, rmax): else: return etot_R - def rhovel_map(rmin, rmax): rho = rho_map(rmin, rmax) @@ -104,8 +100,7 @@ def analysis(iplot: int): rhov_vals = model.render_slice("rhovel", "f64_3", positions) rhoetot_vals = model.render_slice("rhoetot", "f64", positions) - - vx = np.array(rhov_vals)[:,0] / np.array(rho_vals) + vx = np.array(rhov_vals)[:, 0] / np.array(rho_vals) P = (np.array(rhoetot_vals) - 0.5 * np.array(rho_vals) * vx**2) * (gamma - 1) results_dic = { "rho": np.array(rho_vals), @@ -139,7 +134,7 @@ def analysis(iplot: int): # %% # Plot 1: Comparison against analytical solution riemann_solvers = ["rusanov", "hll", "hllc"] -slope_limiters = ["none", "minmod"] +slope_limiters = ["none", "minmod"] xref = 1.0 xrange = 0.5 @@ -156,7 +151,7 @@ def analysis(iplot: int): x_ = arr_x[i] - xref _rho, _vx, _P = sod.get_value(tmax, x_) - #print(x_,_rho, _vx, _P) + # print(x_,_rho, _vx, _P) arr_rho.append(_rho) arr_vx.append(_vx) arr_P.append(_P) @@ -201,7 +196,7 @@ def analysis(iplot: int): ax1.plot(arr_x, delta_rho, label=f"{limiter} {solver} (rho)", linewidth=1) ax2.plot(arr_x, delta_vx, label=f"{limiter} {solver} (vx)", linewidth=1) ax3.plot(arr_x, delta_P, label=f"{limiter} {solver} (P)", linewidth=1) - + ax1.legend() ax2.legend() From 1ff5c9e1680db04c2553e2dea59278629ddd5142 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sat, 14 Feb 2026 20:34:57 +0100 Subject: [PATCH 3/7] [Doc] add double expansion example to sphinx --- .../examples/ramses/run_double_expansion.py | 234 ++++++++++++++++++ 1 file changed, 234 insertions(+) create mode 100644 doc/sphinx/examples/ramses/run_double_expansion.py diff --git a/doc/sphinx/examples/ramses/run_double_expansion.py b/doc/sphinx/examples/ramses/run_double_expansion.py new file mode 100644 index 0000000000..f9869b6109 --- /dev/null +++ b/doc/sphinx/examples/ramses/run_double_expansion.py @@ -0,0 +1,234 @@ +""" +Sod tube test in RAMSES solver +============================================= + +Compare Sod tube with all slope limiters & Riemann solvers +""" + +import os + +import matplotlib.pyplot as plt +import numpy as np + +import shamrock + +tmax = 0.15 +timestamps = 40 +gamma = 1.4 + +multx = 4 +multy = 1 +multz = 1 + +sz = 1 << 1 +base = 16 + +positions = [(x, 0, 0) for x in np.linspace(0.5, 1.5, 256).tolist()[:-1]] + + +def run_advect(slope_limiter: str, riemann_solver: str, only_last_step: bool = True): + ctx = shamrock.Context() + ctx.pdata_layout_new() + + model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3") + + cfg = model.gen_default_config() + scale_fact = 2 / (sz * base * multx) + cfg.set_scale_factor(scale_fact) + cfg.set_eos_gamma(gamma) + + if slope_limiter == "none": + cfg.set_slope_lim_none() + elif slope_limiter == "vanleer": + cfg.set_slope_lim_vanleer_f() + elif slope_limiter == "vanleer_std": + cfg.set_slope_lim_vanleer_std() + elif slope_limiter == "vanleer_sym": + cfg.set_slope_lim_vanleer_sym() + elif slope_limiter == "minmod": + cfg.set_slope_lim_minmod() + else: + raise ValueError(f"Invalid slope limiter: {slope_limiter}") + + if riemann_solver == "rusanov": + cfg.set_riemann_solver_rusanov() + elif riemann_solver == "hll": + cfg.set_riemann_solver_hll() + elif riemann_solver == "hllc": + cfg.set_riemann_solver_hllc() + else: + raise ValueError(f"Invalid Riemann solver: {riemann_solver}") + + cfg.set_face_time_interpolation(True) + model.set_solver_config(cfg) + + model.init_scheduler(int(1e7), 1) + model.make_base_grid((0, 0, 0), (sz, sz, sz), (base * multx, base * multy, base * multz)) + + rho_l = 1 + rho_r = 1 + + v_l = -2 + v_r = 2 + + P_l = 0.4 + P_r = 0.4 + + def rho_map(rmin, rmax): + x, y, z = rmin + if x < 1: + return rho_l + else: + return rho_r + + etot_L = P_l / (gamma - 1) + 0.5 * rho_l * v_l**2 + etot_R = P_r / (gamma - 1) + 0.5 * rho_r * v_r**2 + + def rhoetot_map(rmin, rmax): + rho = rho_map(rmin, rmax) + + x, y, z = rmin + if x < 1: + return etot_L + else: + return etot_R + + def rhovel_map(rmin, rmax): + rho = rho_map(rmin, rmax) + + x, y, z = rmin + if x < 1: + return (rho_l * v_l, 0, 0) + else: + return (rho_r * v_r, 0, 0) + + model.set_field_value_lambda_f64("rho", rho_map) + model.set_field_value_lambda_f64("rhoetot", rhoetot_map) + model.set_field_value_lambda_f64_3("rhovel", rhovel_map) + + results = [] + + def analysis(iplot: int): + rho_vals = model.render_slice("rho", "f64", positions) + rhov_vals = model.render_slice("rhovel", "f64_3", positions) + rhoetot_vals = model.render_slice("rhoetot", "f64", positions) + + vx = np.array(rhov_vals)[:, 0] / np.array(rho_vals) + P = (np.array(rhoetot_vals) - 0.5 * np.array(rho_vals) * vx**2) * (gamma - 1) + results_dic = { + "rho": np.array(rho_vals), + "vx": np.array(vx), + "P": np.array(P), + } + results.append(results_dic) + + print(f"running {slope_limiter} {riemann_solver} with only_last_step={only_last_step}") + + if only_last_step: + model.evolve_until(tmax) + analysis(timestamps) + else: + dt_evolve = tmax / timestamps + + for i in range(timestamps + 1): + model.evolve_until(dt_evolve * i) + analysis(i) + + return results + + +# %% + +cases = [ + ("none", "rusanov"), + ("none", "hll"), + ("none", "hllc"), + ("minmod", "rusanov"), + ("minmod", "hll"), + # ("minmod", "hllc"), # Crash for now +] + +case_anim = "minmod_hll" + +data = {} +for slope_limiter, riemann_solver in cases: + print(f"running {slope_limiter} {riemann_solver}") + key = f"{slope_limiter}_{riemann_solver}" + only_last_step = not (case_anim == key) + data[key] = run_advect(slope_limiter, riemann_solver, only_last_step=only_last_step) + +# %% +# Plot 1: Comparison against analytical solution + +arr_x = [x[0] for x in positions] + +fig, axes = plt.subplots(3, 1, figsize=(6, 15)) +fig.suptitle(f"t={tmax} (Last Step)", fontsize=14) + +for i in range(3): + axes[i].set_xlabel("$x$") + # axes[i].set_yscale("log") + axes[i].grid(True, alpha=0.3) + +ax1, ax2, ax3 = axes +ax1.set_xlabel("$x$") +ax1.set_ylabel("$\\rho$") + +ax2.set_xlabel("$x$") +ax2.set_ylabel("$v_x$") + +ax3.set_xlabel("$x$") +ax3.set_ylabel("$P$") + +for limiter, solver in cases: + key = f"{limiter}_{solver}" + if key in data: + # Get the last timestep + + ax1.plot(arr_x, data[key][-1]["rho"], label=f"{limiter} {solver} (rho)", linewidth=1) + ax2.plot(arr_x, data[key][-1]["vx"], label=f"{limiter} {solver} (vx)", linewidth=1) + ax3.plot(arr_x, data[key][-1]["P"], label=f"{limiter} {solver} (P)", linewidth=1) + + +ax1.legend() +ax2.legend() +ax3.legend() + +plt.tight_layout() +plt.show() + +# %% +# Plot 2: Animation of vanleer_sym_hllc configuration + +# sphinx_gallery_thumbnail_number = 2 + +from matplotlib.animation import FuncAnimation + +fig2, ax2 = plt.subplots() +ax2.set_xlabel("$x$") +ax2.set_ylabel("$\\rho$") +ax2.set_ylim(-2, 2) +ax2.set_xlim(0.5, 1.5) +ax2.grid(True, alpha=0.3) + +(line_rho,) = ax2.plot(arr_x, data[case_anim][0]["rho"], label="rho", linewidth=2) +(line_vx,) = ax2.plot(arr_x, data[case_anim][0]["vx"], label="vx", linewidth=2) +(line_P,) = ax2.plot(arr_x, data[case_anim][0]["P"], label="P", linewidth=2) + + +ax2.legend() +ax2.set_title(f"{case_anim} - t = {0.0:.3f} s") + + +def animate(frame): + t = tmax * frame / timestamps + line_rho.set_ydata(data[case_anim][frame]["rho"]) + line_vx.set_ydata(data[case_anim][frame]["vx"]) + line_P.set_ydata(data[case_anim][frame]["P"]) + ax2.set_title(f"{case_anim} - t = {t:.3f} s") + return (line_rho, line_vx, line_P) + + +anim = FuncAnimation(fig2, animate, frames=timestamps + 1, interval=50, blit=False, repeat=True) +plt.tight_layout() +plt.show() From 1fc5f92a58524744ddb751474ca499c7b9b1be3d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sat, 14 Feb 2026 20:36:55 +0100 Subject: [PATCH 4/7] [Doc] add double expansion example to sphinx --- doc/sphinx/examples/ramses/run_double_expansion.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/doc/sphinx/examples/ramses/run_double_expansion.py b/doc/sphinx/examples/ramses/run_double_expansion.py index f9869b6109..80c657b728 100644 --- a/doc/sphinx/examples/ramses/run_double_expansion.py +++ b/doc/sphinx/examples/ramses/run_double_expansion.py @@ -3,6 +3,16 @@ ============================================= Compare Sod tube with all slope limiters & Riemann solvers +Toro 3rd initial condition + +rho_l = 1 +rho_r = 1 + +v_l = -2 +v_r = 2 + +P_l = 0.4 +P_r = 0.4 """ import os From 4acecb83f3cb51105a9b00bb66dc17f9025bd22a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David--Cl=C3=A9ris=20Timoth=C3=A9e?= Date: Sun, 22 Feb 2026 23:24:40 +0100 Subject: [PATCH 5/7] Delete doc/sphinx/examples/ramses/run_double_expansion.py --- .../examples/ramses/run_double_expansion.py | 244 ------------------ 1 file changed, 244 deletions(-) delete mode 100644 doc/sphinx/examples/ramses/run_double_expansion.py diff --git a/doc/sphinx/examples/ramses/run_double_expansion.py b/doc/sphinx/examples/ramses/run_double_expansion.py deleted file mode 100644 index 80c657b728..0000000000 --- a/doc/sphinx/examples/ramses/run_double_expansion.py +++ /dev/null @@ -1,244 +0,0 @@ -""" -Sod tube test in RAMSES solver -============================================= - -Compare Sod tube with all slope limiters & Riemann solvers -Toro 3rd initial condition - -rho_l = 1 -rho_r = 1 - -v_l = -2 -v_r = 2 - -P_l = 0.4 -P_r = 0.4 -""" - -import os - -import matplotlib.pyplot as plt -import numpy as np - -import shamrock - -tmax = 0.15 -timestamps = 40 -gamma = 1.4 - -multx = 4 -multy = 1 -multz = 1 - -sz = 1 << 1 -base = 16 - -positions = [(x, 0, 0) for x in np.linspace(0.5, 1.5, 256).tolist()[:-1]] - - -def run_advect(slope_limiter: str, riemann_solver: str, only_last_step: bool = True): - ctx = shamrock.Context() - ctx.pdata_layout_new() - - model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3") - - cfg = model.gen_default_config() - scale_fact = 2 / (sz * base * multx) - cfg.set_scale_factor(scale_fact) - cfg.set_eos_gamma(gamma) - - if slope_limiter == "none": - cfg.set_slope_lim_none() - elif slope_limiter == "vanleer": - cfg.set_slope_lim_vanleer_f() - elif slope_limiter == "vanleer_std": - cfg.set_slope_lim_vanleer_std() - elif slope_limiter == "vanleer_sym": - cfg.set_slope_lim_vanleer_sym() - elif slope_limiter == "minmod": - cfg.set_slope_lim_minmod() - else: - raise ValueError(f"Invalid slope limiter: {slope_limiter}") - - if riemann_solver == "rusanov": - cfg.set_riemann_solver_rusanov() - elif riemann_solver == "hll": - cfg.set_riemann_solver_hll() - elif riemann_solver == "hllc": - cfg.set_riemann_solver_hllc() - else: - raise ValueError(f"Invalid Riemann solver: {riemann_solver}") - - cfg.set_face_time_interpolation(True) - model.set_solver_config(cfg) - - model.init_scheduler(int(1e7), 1) - model.make_base_grid((0, 0, 0), (sz, sz, sz), (base * multx, base * multy, base * multz)) - - rho_l = 1 - rho_r = 1 - - v_l = -2 - v_r = 2 - - P_l = 0.4 - P_r = 0.4 - - def rho_map(rmin, rmax): - x, y, z = rmin - if x < 1: - return rho_l - else: - return rho_r - - etot_L = P_l / (gamma - 1) + 0.5 * rho_l * v_l**2 - etot_R = P_r / (gamma - 1) + 0.5 * rho_r * v_r**2 - - def rhoetot_map(rmin, rmax): - rho = rho_map(rmin, rmax) - - x, y, z = rmin - if x < 1: - return etot_L - else: - return etot_R - - def rhovel_map(rmin, rmax): - rho = rho_map(rmin, rmax) - - x, y, z = rmin - if x < 1: - return (rho_l * v_l, 0, 0) - else: - return (rho_r * v_r, 0, 0) - - model.set_field_value_lambda_f64("rho", rho_map) - model.set_field_value_lambda_f64("rhoetot", rhoetot_map) - model.set_field_value_lambda_f64_3("rhovel", rhovel_map) - - results = [] - - def analysis(iplot: int): - rho_vals = model.render_slice("rho", "f64", positions) - rhov_vals = model.render_slice("rhovel", "f64_3", positions) - rhoetot_vals = model.render_slice("rhoetot", "f64", positions) - - vx = np.array(rhov_vals)[:, 0] / np.array(rho_vals) - P = (np.array(rhoetot_vals) - 0.5 * np.array(rho_vals) * vx**2) * (gamma - 1) - results_dic = { - "rho": np.array(rho_vals), - "vx": np.array(vx), - "P": np.array(P), - } - results.append(results_dic) - - print(f"running {slope_limiter} {riemann_solver} with only_last_step={only_last_step}") - - if only_last_step: - model.evolve_until(tmax) - analysis(timestamps) - else: - dt_evolve = tmax / timestamps - - for i in range(timestamps + 1): - model.evolve_until(dt_evolve * i) - analysis(i) - - return results - - -# %% - -cases = [ - ("none", "rusanov"), - ("none", "hll"), - ("none", "hllc"), - ("minmod", "rusanov"), - ("minmod", "hll"), - # ("minmod", "hllc"), # Crash for now -] - -case_anim = "minmod_hll" - -data = {} -for slope_limiter, riemann_solver in cases: - print(f"running {slope_limiter} {riemann_solver}") - key = f"{slope_limiter}_{riemann_solver}" - only_last_step = not (case_anim == key) - data[key] = run_advect(slope_limiter, riemann_solver, only_last_step=only_last_step) - -# %% -# Plot 1: Comparison against analytical solution - -arr_x = [x[0] for x in positions] - -fig, axes = plt.subplots(3, 1, figsize=(6, 15)) -fig.suptitle(f"t={tmax} (Last Step)", fontsize=14) - -for i in range(3): - axes[i].set_xlabel("$x$") - # axes[i].set_yscale("log") - axes[i].grid(True, alpha=0.3) - -ax1, ax2, ax3 = axes -ax1.set_xlabel("$x$") -ax1.set_ylabel("$\\rho$") - -ax2.set_xlabel("$x$") -ax2.set_ylabel("$v_x$") - -ax3.set_xlabel("$x$") -ax3.set_ylabel("$P$") - -for limiter, solver in cases: - key = f"{limiter}_{solver}" - if key in data: - # Get the last timestep - - ax1.plot(arr_x, data[key][-1]["rho"], label=f"{limiter} {solver} (rho)", linewidth=1) - ax2.plot(arr_x, data[key][-1]["vx"], label=f"{limiter} {solver} (vx)", linewidth=1) - ax3.plot(arr_x, data[key][-1]["P"], label=f"{limiter} {solver} (P)", linewidth=1) - - -ax1.legend() -ax2.legend() -ax3.legend() - -plt.tight_layout() -plt.show() - -# %% -# Plot 2: Animation of vanleer_sym_hllc configuration - -# sphinx_gallery_thumbnail_number = 2 - -from matplotlib.animation import FuncAnimation - -fig2, ax2 = plt.subplots() -ax2.set_xlabel("$x$") -ax2.set_ylabel("$\\rho$") -ax2.set_ylim(-2, 2) -ax2.set_xlim(0.5, 1.5) -ax2.grid(True, alpha=0.3) - -(line_rho,) = ax2.plot(arr_x, data[case_anim][0]["rho"], label="rho", linewidth=2) -(line_vx,) = ax2.plot(arr_x, data[case_anim][0]["vx"], label="vx", linewidth=2) -(line_P,) = ax2.plot(arr_x, data[case_anim][0]["P"], label="P", linewidth=2) - - -ax2.legend() -ax2.set_title(f"{case_anim} - t = {0.0:.3f} s") - - -def animate(frame): - t = tmax * frame / timestamps - line_rho.set_ydata(data[case_anim][frame]["rho"]) - line_vx.set_ydata(data[case_anim][frame]["vx"]) - line_P.set_ydata(data[case_anim][frame]["P"]) - ax2.set_title(f"{case_anim} - t = {t:.3f} s") - return (line_rho, line_vx, line_P) - - -anim = FuncAnimation(fig2, animate, frames=timestamps + 1, interval=50, blit=False, repeat=True) -plt.tight_layout() -plt.show() From 68ee82799b87f2b979b5645bd00ebe6afbcb8636 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sun, 22 Feb 2026 23:52:38 +0100 Subject: [PATCH 6/7] update --- doc/sphinx/examples/ramses/run_sod.py | 125 ++++++++++++++------------ src/shamphys/src/SodTube.cpp | 1 - 2 files changed, 66 insertions(+), 60 deletions(-) diff --git a/doc/sphinx/examples/ramses/run_sod.py b/doc/sphinx/examples/ramses/run_sod.py index b5299aef66..1d6d1f949b 100644 --- a/doc/sphinx/examples/ramses/run_sod.py +++ b/doc/sphinx/examples/ramses/run_sod.py @@ -5,6 +5,7 @@ Compare Sod tube with all slope limiters & Riemann solvers """ +import json import os import matplotlib.pyplot as plt @@ -33,33 +34,19 @@ def run_advect(slope_limiter: str, riemann_solver: str, only_last_step: bool = T model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3") cfg = model.gen_default_config() - scale_fact = 2 / (sz * base * multx) - cfg.set_scale_factor(scale_fact) - cfg.set_eos_gamma(gamma) - - if slope_limiter == "none": - cfg.set_slope_lim_none() - elif slope_limiter == "vanleer": - cfg.set_slope_lim_vanleer_f() - elif slope_limiter == "vanleer_std": - cfg.set_slope_lim_vanleer_std() - elif slope_limiter == "vanleer_sym": - cfg.set_slope_lim_vanleer_sym() - elif slope_limiter == "minmod": - cfg.set_slope_lim_minmod() - else: - raise ValueError(f"Invalid slope limiter: {slope_limiter}") - - if riemann_solver == "rusanov": - cfg.set_riemann_solver_rusanov() - elif riemann_solver == "hll": - cfg.set_riemann_solver_hll() - elif riemann_solver == "hllc": - cfg.set_riemann_solver_hllc() - else: - raise ValueError(f"Invalid Riemann solver: {riemann_solver}") - cfg.set_face_time_interpolation(True) + cfg.from_json( + { + "hydro_riemann_solver": riemann_solver, + "slope_limiter": slope_limiter, + "eos_gamma": gamma, + "face_half_time_interpolation": True, + "grid_coord_to_pos_fact": 2 / (sz * base * multx), + } + ) + + print("used config:\n" + json.dumps(cfg.to_json(), indent=4)) + model.set_solver_config(cfg) model.init_scheduler(int(1e7), 1) @@ -96,16 +83,16 @@ def rhovel_map(rmin, rmax): results = [] def analysis(iplot: int): - rho_vals = model.render_slice("rho", "f64", positions) - rhov_vals = model.render_slice("rhovel", "f64_3", positions) - rhoetot_vals = model.render_slice("rhoetot", "f64", positions) + rho_vals = np.array(model.render_slice("rho", "f64", positions)) + rhov_vals = np.array(model.render_slice("rhovel", "f64_3", positions)) + rhoetot_vals = np.array(model.render_slice("rhoetot", "f64", positions)) - vx = np.array(rhov_vals)[:, 0] / np.array(rho_vals) - P = (np.array(rhoetot_vals) - 0.5 * np.array(rho_vals) * vx**2) * (gamma - 1) + vx = rhov_vals[:, 0] / rho_vals + P = (rhoetot_vals - 0.5 * rho_vals * vx**2) * (gamma - 1) results_dic = { - "rho": np.array(rho_vals), - "vx": np.array(vx), - "P": np.array(P), + "rho": rho_vals, + "vx": vx, + "P": P, } results.append(results_dic) @@ -132,7 +119,7 @@ def analysis(iplot: int): data["minmod_hllc"] = run_advect("minmod", "hllc", only_last_step=False) # %% -# Plot 1: Comparison against analytical solution +# Prepare data to be plotted riemann_solvers = ["rusanov", "hll", "hllc"] slope_limiters = ["none", "minmod"] @@ -140,26 +127,36 @@ def analysis(iplot: int): xrange = 0.5 sod = shamrock.phys.SodTube(gamma=gamma, rho_1=1, P_1=1, rho_5=0.125, P_5=0.1) -#### add analytical soluce arr_x = [x[0] for x in positions] -arr_rho = [] -arr_P = [] -arr_vx = [] -for i in range(len(arr_x)): - x_ = arr_x[i] - xref +def get_data_analytic(frame: int): + t = tmax * frame / timestamps + + arr_rho = [] + arr_P = [] + arr_vx = [] + + for i in range(len(arr_x)): + x_ = arr_x[i] - xref + + _rho, _vx, _P = sod.get_value(t, x_) + arr_rho.append(_rho) + arr_vx.append(_vx) + arr_P.append(_P) + + return { + "rho": arr_rho, + "vx": arr_vx, + "P": arr_P, + } - _rho, _vx, _P = sod.get_value(tmax, x_) - # print(x_,_rho, _vx, _P) - arr_rho.append(_rho) - arr_vx.append(_vx) - arr_P.append(_P) -arr_rho = np.array(arr_rho) -arr_vx = np.array(arr_vx) -arr_P = np.array(arr_P) +data_analytic = [get_data_analytic(i) for i in range(timestamps + 1)] + +# %% +# Plot 1: Comparison against analytical solution fig, axes = plt.subplots(3, 1, figsize=(6, 15)) fig.suptitle(f"t={tmax} (Last Step)", fontsize=14) @@ -170,18 +167,16 @@ def analysis(iplot: int): ax1, ax2, ax3 = axes ax1.set_xlabel("$x$") -ax1.set_ylabel("$\\rho$") +ax1.set_ylabel("$\\delta \\rho$") ax2.set_xlabel("$x$") -ax2.set_ylabel("$v_x$") +ax2.set_ylabel("$\\delta v_x$") ax3.set_xlabel("$x$") -ax3.set_ylabel("$P$") +ax3.set_ylabel("$\\delta P$") -# ax1.plot(arr_x, arr_rho, color="black", label="analytic") -# ax2.plot(arr_x, arr_vx, color="black", label="analytic") -# ax3.plot(arr_x, arr_P, color="black", label="analytic") +last_analytic = data_analytic[-1] for limiter in slope_limiters: @@ -189,9 +184,9 @@ def analysis(iplot: int): key = f"{limiter}_{solver}" if key in data: # Get the last timestep - delta_rho = np.abs(data[key][-1]["rho"] - arr_rho) - delta_vx = np.abs(data[key][-1]["vx"] - arr_vx) - delta_P = np.abs(data[key][-1]["P"] - arr_P) + delta_rho = np.abs(data[key][-1]["rho"] - last_analytic["rho"]) + delta_vx = np.abs(data[key][-1]["vx"] - last_analytic["vx"]) + delta_P = np.abs(data[key][-1]["P"] - last_analytic["P"]) ax1.plot(arr_x, delta_rho, label=f"{limiter} {solver} (rho)", linewidth=1) ax2.plot(arr_x, delta_vx, label=f"{limiter} {solver} (vx)", linewidth=1) @@ -219,6 +214,15 @@ def analysis(iplot: int): ax2.set_xlim(0.5, 1.5) ax2.grid(True, alpha=0.3) +(line_analytic_rho,) = ax2.plot( + arr_x, data_analytic[0]["rho"], color="black", label="analytic", linewidth=2 +) +(line_analytic_vx,) = ax2.plot( + arr_x, data_analytic[0]["vx"], color="black", label="analytic", linewidth=2 +) +(line_analytic_P,) = ax2.plot( + arr_x, data_analytic[0]["P"], color="black", label="analytic", linewidth=2 +) (line_rho,) = ax2.plot(arr_x, data["minmod_hllc"][0]["rho"], label="rho", linewidth=2) (line_vx,) = ax2.plot(arr_x, data["minmod_hllc"][0]["vx"], label="vx", linewidth=2) (line_P,) = ax2.plot(arr_x, data["minmod_hllc"][0]["P"], label="P", linewidth=2) @@ -232,8 +236,11 @@ def animate(frame): line_rho.set_ydata(data["minmod_hllc"][frame]["rho"]) line_vx.set_ydata(data["minmod_hllc"][frame]["vx"]) line_P.set_ydata(data["minmod_hllc"][frame]["P"]) + line_analytic_rho.set_ydata(data_analytic[frame]["rho"]) + line_analytic_vx.set_ydata(data_analytic[frame]["vx"]) + line_analytic_P.set_ydata(data_analytic[frame]["P"]) ax2.set_title(f"minmod_hllc - t = {t:.3f} s") - return (line_rho, line_vx, line_P) + return (line_rho, line_vx, line_P, line_analytic_rho, line_analytic_vx, line_analytic_P) anim = FuncAnimation(fig2, animate, frames=timestamps + 1, interval=50, blit=False, repeat=True) diff --git a/src/shamphys/src/SodTube.cpp b/src/shamphys/src/SodTube.cpp index 9977ca2d9f..607be9ae5c 100644 --- a/src/shamphys/src/SodTube.cpp +++ b/src/shamphys/src/SodTube.cpp @@ -14,7 +14,6 @@ */ #include "shamphys/SodTube.hpp" -#include "shamcomm/logs.hpp" #include "shammath/derivatives.hpp" #include "shammath/solve.hpp" #include From 401d4c5cf2c3149404f34dcebd67bc4562393c7e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20David--Cl=C3=A9ris?= Date: Sun, 22 Feb 2026 23:55:09 +0100 Subject: [PATCH 7/7] split sod test away from main job --- .github/workflows/shamrock-acpp-phys-test.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/shamrock-acpp-phys-test.yml b/.github/workflows/shamrock-acpp-phys-test.yml index efb6cf725d..3054bf9506 100644 --- a/.github/workflows/shamrock-acpp-phys-test.yml +++ b/.github/workflows/shamrock-acpp-phys-test.yml @@ -143,6 +143,8 @@ jobs: artifact_name: gallery_ramses_run_linear_wave_with_bc - testfile: "examples/ramses/run_advect.py" artifact_name: gallery_ramses_run_advect + - testfile: "examples/ramses/run_sod.py" + artifact_name: gallery_ramses_run_sod - testfile: "examples/ramses/run_toro_shocks.py" artifact_name: gallery_ramses_run_toro_shocks