Note

Go to the end to download the full example code.

Testing Sod tube with GSPH (exact Riemann solver + Inutsuka V2 force)#

CI test for Sod tube with GSPH using M4 kernel, the exact Riemann solver (Toro 2009), and the Inutsuka (2002) effective volume/face force formulation. Uses piecewise constant reconstruction (first-order, stable). Computes L2 error against analytical solution and checks for regression.

 11 import numpy as np
 12
 13 import shamrock
 14
 15 gamma = 1.4
 16 rho_L, rho_R = 1.0, 0.125
 17 P_L, P_R = 1.0, 0.1
 18 fact = (rho_L / rho_R) ** (1.0 / 3.0)
 19 u_L = P_L / ((gamma - 1) * rho_L)
 20 u_R = P_R / ((gamma - 1) * rho_R)
 21 resol = 64
 22
 23 ctx = shamrock.Context()
 24 ctx.pdata_layout_new()
 25
 26 model = shamrock.get_Model_GSPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
 27 cfg = model.gen_default_config()
 28 cfg.set_riemann_exact()
 29 cfg.set_force_inutsuka_v2()
 30 cfg.set_reconstruct_piecewise_constant()
 31 cfg.set_boundary_periodic()
 32 cfg.set_eos_adiabatic(gamma)
 33 cfg.print_status()
 34 model.set_solver_config(cfg)
 35 model.init_scheduler(int(1e8), 1)
 36
 37 (xs, ys, zs) = model.get_box_dim_fcc_3d(1, resol, 12, 12)
 38 dr = 1 / xs
 39 (xs, ys, zs) = model.get_box_dim_fcc_3d(dr, resol, 12, 12)
 40 model.resize_simulation_box((-xs, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
 41
 42 model.add_cube_hcp_3d(dr, (-xs, -ys / 2, -zs / 2), (0, ys / 2, zs / 2))
 43 model.add_cube_hcp_3d(dr * fact, (0, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
 44 model.set_field_in_box("uint", "f64", u_L, (-xs, -ys / 2, -zs / 2), (0, ys / 2, zs / 2))
 45 model.set_field_in_box("uint", "f64", u_R, (0, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
 46
 47 vol_b = xs * ys * zs
 48 totmass = (rho_R * vol_b) + (rho_L * vol_b)
 49 pmass = model.total_mass_to_part_mass(totmass)
 50 model.set_particle_mass(pmass)
 51 hfact = model.get_hfact()
 52
 53 model.set_cfl_cour(0.1)
 54 model.set_cfl_force(0.1)
 55
 56 t_target = 0.245
 57 print(f"GSPH Sod Shock Tube Test (M4, Exact Riemann, Inutsuka V2, t={t_target})")
 58 model.evolve_until(t_target)
 59
 60 sod = shamrock.phys.SodTube(gamma=gamma, rho_1=rho_L, P_1=P_L, rho_5=rho_R, P_5=P_R)
 61
 62 data = ctx.collect_data()
 63
 64
 65 def compute_L2_errors(data, sod, t, x_min, x_max):
 66     """Compute L2 errors using ctx.collect_data() (no pyvista dependency)."""
 67     points = np.array(data["xyz"])
 68     velocities = np.array(data["vxyz"])
 69     hpart = np.array(data["hpart"])
 70     uint = np.array(data["uint"])
 71
 72     rho_sim = pmass * (hfact / hpart) ** 3
 73     P_sim = (gamma - 1) * rho_sim * uint
 74
 75     x, vx, vy, vz = points[:, 0], velocities[:, 0], velocities[:, 1], velocities[:, 2]
 76     mask = (x >= x_min) & (x <= x_max)
 77     x_f, rho_f, vx_f, vy_f, vz_f, P_f = (
 78         x[mask],
 79         rho_sim[mask],
 80         vx[mask],
 81         vy[mask],
 82         vz[mask],
 83         P_sim[mask],
 84     )
 85
 86     if len(x_f) == 0:
 87         raise RuntimeError("No particles in analysis region")
 88
 89     rho_ana, vx_ana, P_ana = np.zeros(len(x_f)), np.zeros(len(x_f)), np.zeros(len(x_f))
 90     for i, xi in enumerate(x_f):
 91         rho_ana[i], vx_ana[i], P_ana[i] = sod.get_value(t, xi)
 92
 93     err_rho = np.sqrt(np.mean((rho_f - rho_ana) ** 2)) / np.mean(rho_ana)
 94     err_vx = np.sqrt(np.mean((vx_f - vx_ana) ** 2)) / (np.mean(np.abs(vx_ana)) + 0.1)
 95     err_vy = np.sqrt(np.mean(vy_f**2))
 96     err_vz = np.sqrt(np.mean(vz_f**2))
 97     err_P = np.sqrt(np.mean((P_f - P_ana) ** 2)) / np.mean(P_ana)
 98     return err_rho, (err_vx, err_vy, err_vz), err_P
 99
100
101 if shamrock.sys.world_rank() == 0:
102     rho, v, P = compute_L2_errors(data, sod, t_target, -0.5, 0.5)
103     vx, vy, vz = v
104
105     print("current errors :")
106     print(f"err_rho = {rho}")
107     print(f"err_vx = {vx}")
108     print(f"err_vy = {vy}")
109     print(f"err_vz = {vz}")
110     print(f"err_P = {P}")
111
112     # Expected L2 error values (calibrated from CI's own run with M4 kernel,
113     # exact Riemann solver, Inutsuka V2 force formulation)
114     # Tolerance set very strict for regression testing (like sod_tube_sph.py).
115     # NOTE: calibrated on CI hardware (Linux x86_64 docker image), not a local
116     # macOS/arm64 build -- the two disagree at the ULP level because the exact
117     # solver's convergence-based bisection early-exit is sensitive to libm's
118     # sqrt/pow rounding. If CI's actual output ever differs again, update
119     # these expect_* values from CI's own printed "current errors" output.
120     expect_rho = 0.05004380788795549
121     expect_vx = 0.1429980515258284
122     expect_vy = 0.005331795999290027
123     expect_vz = 7.752615518913349e-05
124     expect_P = 0.06285331285974535
125
126     tol = 1e-8
127
128     test_pass = True
129     err_log = ""
130
131     error_checks = {
132         "rho": (rho, expect_rho),
133         "vx": (vx, expect_vx),
134         "vy": (vy, expect_vy),
135         "vz": (vz, expect_vz),
136         "P": (P, expect_P),
137     }
138
139     for name, (value, expected) in error_checks.items():
140         if abs(value - expected) > tol * expected:
141             err_log += f"error on {name} is outside of tolerances:\n"
142             err_log += f"  expected error = {expected} +- {tol * expected}\n"
143             err_log += (
144                 f"  obtained error = {value} (relative error = {(value - expected) / expected})\n"
145             )
146             test_pass = False
147
148     if test_pass:
149         print("\n" + "=" * 50)
150         print("GSPH Sod Shock Tube Test (Exact + InutsukaV2): PASSED")
151         print("=" * 50)
152     else:
153         print("\n" + "=" * 50)
154         print("GSPH Sod Shock Tube Test (Exact + InutsukaV2): FAILED")
155         print("=" * 50)
156         print(err_log)
157         raise RuntimeError("Test failed:\n" + err_log)

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery