Note

Go to the end to download the full example code.

Testing Sod tube with GSPH#

CI test for Sod tube with GSPH using M4 kernel and HLLC Riemann solver. Uses piecewise constant reconstruction (first-order, stable). Computes L2 error against analytical solution and checks for regression.

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

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery