Skip to content

Commit 97a965e

Browse files
Guo-astroGuo-astro
committed
Add GSPH unit tests and CI physics validation tests
- Add unit tests for GSPH modules: Riemann solvers, forces, integration - Add CI physics tests for Sod shock tube and extreme blast wave - Use ctx.collect_data() for direct memory access (no pyvista dependency) - Use strict tolerances (1e-8) for regression testing - Fix M4 kernel normalization constants - Remove non-existent MUSCL test from CI workflow
1 parent ce69ca0 commit 97a965e

17 files changed

Lines changed: 2220 additions & 174 deletions

File tree

‎.github/workflows/shamrock-acpp-phys-test.yml‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,10 @@ jobs:
5151
include:
5252
- worldsize: 1
5353
testfile: sod_tube_sph.py
54+
- worldsize: 1
55+
testfile: gsph/sod_tube_gsph.py
56+
- worldsize: 2
57+
testfile: gsph/sod_tube_gsph.py
5458
- worldsize: 1
5559
testfile: comp_phantom_sedov_1patch.py
5660
- worldsize: 2

‎doc/sphinx/examples/tests_ci/gsph/blast_wave_gsph.py‎

Lines changed: 147 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,147 @@
1+
"""
2+
Testing Extreme Blast Wave with GSPH
3+
====================================
4+
5+
CI test for the extreme blast wave problem from Inutsuka 2002 (Section 4.3).
6+
This is a severe test with Mach number ~10^5.
7+
8+
Initial conditions (from Inutsuka 2002):
9+
rho_L = 1, rho_R = 1
10+
P_L = 3000, P_R = 1e-7
11+
v_L = 0, v_R = 0
12+
"""
13+
14+
import numpy as np
15+
16+
import shamrock
17+
18+
gamma = 1.4
19+
rho_L, rho_R = 1.0, 1.0
20+
P_L, P_R = 3000.0, 1e-7
21+
u_L = P_L / ((gamma - 1) * rho_L)
22+
u_R = P_R / ((gamma - 1) * rho_R)
23+
resol = 100
24+
25+
ctx = shamrock.Context()
26+
ctx.pdata_layout_new()
27+
28+
model = shamrock.get_Model_GSPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
29+
cfg = model.gen_default_config()
30+
cfg.set_riemann_hllc()
31+
cfg.set_reconstruct_piecewise_constant()
32+
cfg.set_boundary_periodic()
33+
cfg.set_eos_adiabatic(gamma)
34+
cfg.print_status()
35+
model.set_solver_config(cfg)
36+
model.init_scheduler(int(1e8), 1)
37+
38+
(xs, ys, zs) = model.get_box_dim_fcc_3d(1, resol, 24, 24)
39+
dr = 1 / xs
40+
(xs, ys, zs) = model.get_box_dim_fcc_3d(dr, resol, 24, 24)
41+
model.resize_simulation_box((-xs, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
42+
43+
model.add_cube_hcp_3d(dr, (-xs, -ys / 2, -zs / 2), (0, ys / 2, zs / 2))
44+
model.add_cube_hcp_3d(dr, (0, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
45+
model.set_field_in_box("uint", "f64", u_L, (-xs, -ys / 2, -zs / 2), (0, ys / 2, zs / 2))
46+
model.set_field_in_box("uint", "f64", u_R, (0, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
47+
48+
vol_b = xs * ys * zs
49+
totmass = (rho_R * vol_b) + (rho_L * vol_b)
50+
pmass = model.total_mass_to_part_mass(totmass)
51+
model.set_particle_mass(pmass)
52+
hfact = model.get_hfact()
53+
54+
model.set_cfl_cour(0.3)
55+
model.set_cfl_force(0.3)
56+
57+
t_target = 0.015
58+
print(f"GSPH Extreme Blast Wave Test (M4, HLLC, t={t_target})")
59+
model.evolve_until(t_target)
60+
61+
sod = shamrock.phys.SodTube(gamma=gamma, rho_1=rho_L, P_1=P_L, rho_5=rho_R, P_5=P_R)
62+
63+
64+
def compute_L2_errors(ctx, sod, t, x_min, x_max):
65+
"""Compute L2 errors using ctx.collect_data() (no pyvista dependency)."""
66+
data = ctx.collect_data()
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+
rho_norm = max(np.mean(rho_ana), 1e-10)
94+
vx_norm = max(np.mean(np.abs(vx_ana)), 0.1)
95+
P_norm = max(np.mean(P_ana), 1e-10)
96+
97+
err_rho = np.sqrt(np.mean((rho_f - rho_ana) ** 2)) / rho_norm
98+
err_vx = np.sqrt(np.mean((vx_f - vx_ana) ** 2)) / vx_norm
99+
err_vy = np.sqrt(np.mean(vy_f**2))
100+
err_vz = np.sqrt(np.mean(vz_f**2))
101+
err_P = np.sqrt(np.mean((P_f - P_ana) ** 2)) / P_norm
102+
return err_rho, (err_vx, err_vy, err_vz), err_P
103+
104+
105+
rho, v, P = compute_L2_errors(ctx, sod, t_target, -0.5, 0.5)
106+
vx, vy, vz = v
107+
108+
print("current errors :")
109+
print(f"err_rho = {rho}")
110+
print(f"err_vx = {vx}")
111+
print(f"err_vy = {vy}")
112+
print(f"err_vz = {vz}")
113+
print(f"err_P = {P}")
114+
115+
# Expected L2 error values (calibrated from CI run with M4 kernel)
116+
expect_rho = 10.688658207003348
117+
expect_vx = 1.0420471749025182
118+
expect_vy = 0.11766417324542999
119+
expect_vz = 0.0027436730451881886
120+
expect_P = 1.6660643954434153
121+
122+
tol = 1e-8
123+
124+
test_pass = True
125+
err_log = ""
126+
127+
error_checks = {
128+
"rho": (rho, expect_rho),
129+
"vx": (vx, expect_vx),
130+
"vy": (vy, expect_vy),
131+
"vz": (vz, expect_vz),
132+
"P": (P, expect_P),
133+
}
134+
135+
for name, (value, expected) in error_checks.items():
136+
if abs(value - expected) > tol * expected:
137+
err_log += f"error on {name} is outside of tolerances:\n"
138+
err_log += f" expected error = {expected} +- {tol*expected}\n"
139+
err_log += f" obtained error = {value} (relative error = {(value - expected)/expected})\n"
140+
test_pass = False
141+
142+
if test_pass:
143+
print("\n" + "=" * 50)
144+
print("GSPH Extreme Blast Wave Test: PASSED")
145+
print("=" * 50)
146+
else:
147+
exit("Test did not pass L2 margins : \n" + err_log)

‎doc/sphinx/examples/tests_ci/gsph/sod_tube_gsph.py‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
"""
2+
Testing Sod tube with GSPH
3+
==========================
4+
5+
CI test for Sod tube with GSPH using M4 kernel and HLLC Riemann solver.
6+
Uses piecewise constant reconstruction (first-order, stable).
7+
Computes L2 error against analytical solution and checks for regression.
8+
"""
9+
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+
61+
def compute_L2_errors(ctx, sod, t, x_min, x_max):
62+
"""Compute L2 errors using ctx.collect_data() (no pyvista dependency)."""
63+
data = ctx.collect_data()
64+
points = np.array(data["xyz"])
65+
velocities = np.array(data["vxyz"])
66+
hpart = np.array(data["hpart"])
67+
uint = np.array(data["uint"])
68+
69+
rho_sim = pmass * (hfact / hpart) ** 3
70+
P_sim = (gamma - 1) * rho_sim * uint
71+
72+
x, vx, vy, vz = points[:, 0], velocities[:, 0], velocities[:, 1], velocities[:, 2]
73+
mask = (x >= x_min) & (x <= x_max)
74+
x_f, rho_f, vx_f, vy_f, vz_f, P_f = (
75+
x[mask],
76+
rho_sim[mask],
77+
vx[mask],
78+
vy[mask],
79+
vz[mask],
80+
P_sim[mask],
81+
)
82+
83+
if len(x_f) == 0:
84+
raise RuntimeError("No particles in analysis region")
85+
86+
rho_ana, vx_ana, P_ana = np.zeros(len(x_f)), np.zeros(len(x_f)), np.zeros(len(x_f))
87+
for i, xi in enumerate(x_f):
88+
rho_ana[i], vx_ana[i], P_ana[i] = sod.get_value(t, xi)
89+
90+
err_rho = np.sqrt(np.mean((rho_f - rho_ana) ** 2)) / np.mean(rho_ana)
91+
err_vx = np.sqrt(np.mean((vx_f - vx_ana) ** 2)) / (np.mean(np.abs(vx_ana)) + 0.1)
92+
err_vy = np.sqrt(np.mean(vy_f**2))
93+
err_vz = np.sqrt(np.mean(vz_f**2))
94+
err_P = np.sqrt(np.mean((P_f - P_ana) ** 2)) / np.mean(P_ana)
95+
return err_rho, (err_vx, err_vy, err_vz), err_P
96+
97+
98+
rho, v, P = compute_L2_errors(ctx, sod, t_target, -0.5, 0.5)
99+
vx, vy, vz = v
100+
101+
print("current errors :")
102+
print(f"err_rho = {rho}")
103+
print(f"err_vx = {vx}")
104+
print(f"err_vy = {vy}")
105+
print(f"err_vz = {vz}")
106+
print(f"err_P = {P}")
107+
108+
# Expected L2 error values (calibrated from CI run with M4 kernel)
109+
# Tolerance set very strict for regression testing (like sod_tube_sph.py)
110+
expect_rho = 0.03053641590859224
111+
expect_vx = 0.10520433643658435
112+
expect_vy = 0.0002224532508989796
113+
expect_vz = 5.029288149671629e-05
114+
expect_P = 0.037878823126488034
115+
116+
tol = 1e-8
117+
118+
test_pass = True
119+
err_log = ""
120+
121+
error_checks = {
122+
"rho": (rho, expect_rho),
123+
"vx": (vx, expect_vx),
124+
"vy": (vy, expect_vy),
125+
"vz": (vz, expect_vz),
126+
"P": (P, expect_P),
127+
}
128+
129+
for name, (value, expected) in error_checks.items():
130+
if abs(value - expected) > tol * expected:
131+
err_log += f"error on {name} is outside of tolerances:\n"
132+
err_log += f" expected error = {expected} +- {tol*expected}\n"
133+
err_log += f" obtained error = {value} (relative error = {(value - expected)/expected})\n"
134+
test_pass = False
135+
136+
if test_pass:
137+
print("\n" + "=" * 50)
138+
print("GSPH Sod Shock Tube Test: PASSED")
139+
print("=" * 50)
140+
else:
141+
exit("Test did not pass L2 margins : \n" + err_log)

‎src/shambase/include/shambase/constants.hpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@
1111

1212
/**
1313
* @file constants.hpp
14+
* @author Guo Yansong (guo.yansong.ngy@gmail.com)
1415
* @author Timothée David--Cléris (tim.shamrock@proton.me)
1516
* @brief Class holding the value of numerous constants
1617
* generated from the following source
@@ -40,6 +41,7 @@
4041
add("gamma_5_6", gamma(5/6))
4142
add("gamma_1", gamma(1))
4243
add("sqrt_2", 2**(1/2))
44+
add("sqrt_pi", pi**(1/2))
4345
add("e", e)
4446
*/
4547

@@ -65,6 +67,7 @@ namespace shambase::constants {
6567
template<class T> constexpr T gamma_5_6 = 1.1287870299081257386;
6668
template<class T> constexpr T gamma_1 = 1;
6769
template<class T> constexpr T sqrt_2 = 1.4142135623730951455;
70+
template<class T> constexpr T sqrt_pi = 1.7724538509055158819; // same as gamma_1_2
6871
template<class T> constexpr T e = 2.7182818284590450908;
6972
// clang-format on
7073

‎src/shammath/include/shammath/sphkernels.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@
1111

1212
/**
1313
* @file sphkernels.hpp
14+
* @author Guo Yansong (guo.yansong.ngy@gmail.com)
1415
* @author Timothée David--Cléris (tim.shamrock@proton.me)
1516
* @brief sph kernels
1617
*/

‎src/shammodels/gsph/CMakeLists.txt‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,11 +11,12 @@ cmake_minimum_required(VERSION 3.9)
1111

1212
project(Shammodels_gsph CXX C)
1313

14-
# Sources: GSPH solver, Model, and VTK I/O
14+
# Sources: GSPH solver, Model, VTK I/O, and Python bindings
1515
set(Sources
1616
src/SolverConfig.cpp
1717
src/Solver.cpp
1818
src/Model.cpp
19+
src/pyGSPHModel.cpp
1920
src/modules/UpdateDerivs.cpp
2021
src/modules/io/VTKDump.cpp
2122
)

‎src/shammodels/gsph/include/shammodels/gsph/Solver.hpp‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -130,6 +130,17 @@ namespace shammodels::gsph {
130130
void compute_eos_fields();
131131
void reset_eos_fields();
132132

133+
/**
134+
* @brief Compute gradients for MUSCL reconstruction
135+
*
136+
* Computes density, pressure, and velocity gradients for each particle
137+
* using SPH kernel gradient summation. Only called when MUSCL reconstruction
138+
* is enabled (reconstruct_config.is_muscl() == true).
139+
*
140+
* Reference: Cha & Whitworth (2003)
141+
*/
142+
void compute_gradients();
143+
133144
void prepare_corrector();
134145

135146
/**

0 commit comments

Comments
 (0)