Add GSPH unit tests and CI physics validation tests

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

@@ -51,6 +51,10 @@ jobs:
         include:
           - worldsize: 1
             testfile: sod_tube_sph.py
+          - worldsize: 1
+            testfile: gsph/sod_tube_gsph.py
+          - worldsize: 2
+            testfile: gsph/sod_tube_gsph.py
           - worldsize: 1
             testfile: comp_phantom_sedov_1patch.py
           - worldsize: 2

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

@@ -0,0 +1,151 @@
+"""
+Testing Extreme Blast Wave with GSPH
+====================================
+
+CI test for the extreme blast wave problem from Inutsuka 2002 (Section 4.3).
+This is a severe test with Mach number ~10^5.
+
+Initial conditions (from Inutsuka 2002):
+    rho_L = 1,      rho_R = 1
+    P_L   = 3000,   P_R   = 1e-7
+    v_L   = 0,      v_R   = 0
+"""
+
+import numpy as np
+
+import shamrock
+
+gamma = 1.4
+rho_L, rho_R = 1.0, 1.0
+P_L, P_R = 3000.0, 1e-7
+u_L = P_L / ((gamma - 1) * rho_L)
+u_R = P_R / ((gamma - 1) * rho_R)
+resol = 100
+
+ctx = shamrock.Context()
+ctx.pdata_layout_new()
+
+model = shamrock.get_Model_GSPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
+cfg = model.gen_default_config()
+cfg.set_riemann_hllc()
+cfg.set_reconstruct_piecewise_constant()
+cfg.set_boundary_periodic()
+cfg.set_eos_adiabatic(gamma)
+cfg.print_status()
+model.set_solver_config(cfg)
+model.init_scheduler(int(1e8), 1)
+
+(xs, ys, zs) = model.get_box_dim_fcc_3d(1, resol, 24, 24)
+dr = 1 / xs
+(xs, ys, zs) = model.get_box_dim_fcc_3d(dr, resol, 24, 24)
+model.resize_simulation_box((-xs, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
+
+model.add_cube_hcp_3d(dr, (-xs, -ys / 2, -zs / 2), (0, ys / 2, zs / 2))
+model.add_cube_hcp_3d(dr, (0, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
+model.set_field_in_box("uint", "f64", u_L, (-xs, -ys / 2, -zs / 2), (0, ys / 2, zs / 2))
+model.set_field_in_box("uint", "f64", u_R, (0, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
+
+vol_b = xs * ys * zs
+totmass = (rho_R * vol_b) + (rho_L * vol_b)
+pmass = model.total_mass_to_part_mass(totmass)
+model.set_particle_mass(pmass)
+hfact = model.get_hfact()
+
+model.set_cfl_cour(0.3)
+model.set_cfl_force(0.3)
+
+t_target = 0.015
+print(f"GSPH Extreme Blast Wave Test (M4, HLLC, t={t_target})")
+model.evolve_until(t_target)
+
+sod = shamrock.phys.SodTube(gamma=gamma, rho_1=rho_L, P_1=P_L, rho_5=rho_R, P_5=P_R)
+
+data = ctx.collect_data()
+
+
+def compute_L2_errors(data, sod, t, x_min, x_max):
+    """Compute L2 errors using ctx.collect_data() (no pyvista dependency)."""
+    points = np.array(data["xyz"])
+    velocities = np.array(data["vxyz"])
+    hpart = np.array(data["hpart"])
+    uint = np.array(data["uint"])
+
+    rho_sim = pmass * (hfact / hpart) ** 3
+    P_sim = (gamma - 1) * rho_sim * uint
+
+    x, vx, vy, vz = points[:, 0], velocities[:, 0], velocities[:, 1], velocities[:, 2]
+    mask = (x >= x_min) & (x <= x_max)
+    x_f, rho_f, vx_f, vy_f, vz_f, P_f = (
+        x[mask],
+        rho_sim[mask],
+        vx[mask],
+        vy[mask],
+        vz[mask],
+        P_sim[mask],
+    )
+
+    if len(x_f) == 0:
+        raise RuntimeError("No particles in analysis region")
+
+    rho_ana, vx_ana, P_ana = np.zeros(len(x_f)), np.zeros(len(x_f)), np.zeros(len(x_f))
+    for i, xi in enumerate(x_f):
+        rho_ana[i], vx_ana[i], P_ana[i] = sod.get_value(t, xi)
+
+    rho_norm = max(np.mean(rho_ana), 1e-10)
+    vx_norm = max(np.mean(np.abs(vx_ana)), 0.1)
+    P_norm = max(np.mean(P_ana), 1e-10)
+
+    err_rho = np.sqrt(np.mean((rho_f - rho_ana) ** 2)) / rho_norm
+    err_vx = np.sqrt(np.mean((vx_f - vx_ana) ** 2)) / vx_norm
+    err_vy = np.sqrt(np.mean(vy_f**2))
+    err_vz = np.sqrt(np.mean(vz_f**2))
+    err_P = np.sqrt(np.mean((P_f - P_ana) ** 2)) / P_norm
+    return err_rho, (err_vx, err_vy, err_vz), err_P
+
+
+if shamrock.sys.world_rank() == 0:
+    rho, v, P = compute_L2_errors(ctx, sod, t_target, -0.5, 0.5)
+    vx, vy, vz = v
+
+    print("current errors :")
+    print(f"err_rho = {rho}")
+    print(f"err_vx = {vx}")
+    print(f"err_vy = {vy}")
+    print(f"err_vz = {vz}")
+    print(f"err_P = {P}")
+
+    # Expected L2 error values (calibrated from CI run with M4 kernel)
+    expect_rho = 10.688658207003348
+    expect_vx = 1.0420471749025182
+    expect_vy = 0.11766417324542999
+    expect_vz = 0.0027436730451881886
+    expect_P = 1.6660643954434153
+
+    tol = 1e-8
+
+    test_pass = True
+    err_log = ""
+
+    error_checks = {
+        "rho": (rho, expect_rho),
+        "vx": (vx, expect_vx),
+        "vy": (vy, expect_vy),
+        "vz": (vz, expect_vz),
+        "P": (P, expect_P),
+    }
+
+    for name, (value, expected) in error_checks.items():
+        if abs(value - expected) > tol * expected:
+            err_log += f"error on {name} is outside of tolerances:\n"
+            err_log += f"  expected error = {expected} +- {tol*expected}\n"
+            err_log += (
+                f"  obtained error = {value} (relative error = {(value - expected)/expected})\n"
+            )
+            test_pass = False
+
+    if test_pass:
+        print("\n" + "=" * 50)
+        print("GSPH Extreme Blast Wave Test: PASSED")
+        print("=" * 50)
+    else:
+        exit("Test did not pass L2 margins : \n" + err_log)

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

@@ -0,0 +1,143 @@
+"""
+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.
+"""
+
+import numpy as np
+
+import shamrock
+
+gamma = 1.4
+rho_L, rho_R = 1.0, 0.125
+P_L, P_R = 1.0, 0.1
+fact = (rho_L / rho_R) ** (1.0 / 3.0)
+u_L = P_L / ((gamma - 1) * rho_L)
+u_R = P_R / ((gamma - 1) * rho_R)
+resol = 128
+
+ctx = shamrock.Context()
+ctx.pdata_layout_new()
+
+model = shamrock.get_Model_GSPH(context=ctx, vector_type="f64_3", sph_kernel="M4")
+cfg = model.gen_default_config()
+cfg.set_riemann_hllc()
+cfg.set_reconstruct_piecewise_constant()
+cfg.set_boundary_periodic()
+cfg.set_eos_adiabatic(gamma)
+cfg.print_status()
+model.set_solver_config(cfg)
+model.init_scheduler(int(1e8), 1)
+
+(xs, ys, zs) = model.get_box_dim_fcc_3d(1, resol, 24, 24)
+dr = 1 / xs
+(xs, ys, zs) = model.get_box_dim_fcc_3d(dr, resol, 24, 24)
+model.resize_simulation_box((-xs, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
+
+model.add_cube_hcp_3d(dr, (-xs, -ys / 2, -zs / 2), (0, ys / 2, zs / 2))
+model.add_cube_hcp_3d(dr * fact, (0, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
+model.set_field_in_box("uint", "f64", u_L, (-xs, -ys / 2, -zs / 2), (0, ys / 2, zs / 2))
+model.set_field_in_box("uint", "f64", u_R, (0, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
+
+vol_b = xs * ys * zs
+totmass = (rho_R * vol_b) + (rho_L * vol_b)
+pmass = model.total_mass_to_part_mass(totmass)
+model.set_particle_mass(pmass)
+hfact = model.get_hfact()
+
+model.set_cfl_cour(0.1)
+model.set_cfl_force(0.1)
+
+t_target = 0.245
+print(f"GSPH Sod Shock Tube Test (M4, HLLC, t={t_target})")
+model.evolve_until(t_target)
+
+sod = shamrock.phys.SodTube(gamma=gamma, rho_1=rho_L, P_1=P_L, rho_5=rho_R, P_5=P_R)
+
+data = ctx.collect_data()
+
+
+def compute_L2_errors(data, sod, t, x_min, x_max):
+    """Compute L2 errors using ctx.collect_data() (no pyvista dependency)."""
+    points = np.array(data["xyz"])
+    velocities = np.array(data["vxyz"])
+    hpart = np.array(data["hpart"])
+    uint = np.array(data["uint"])
+
+    rho_sim = pmass * (hfact / hpart) ** 3
+    P_sim = (gamma - 1) * rho_sim * uint
+
+    x, vx, vy, vz = points[:, 0], velocities[:, 0], velocities[:, 1], velocities[:, 2]
+    mask = (x >= x_min) & (x <= x_max)
+    x_f, rho_f, vx_f, vy_f, vz_f, P_f = (
+        x[mask],
+        rho_sim[mask],
+        vx[mask],
+        vy[mask],
+        vz[mask],
+        P_sim[mask],
+    )
+
+    if len(x_f) == 0:
+        raise RuntimeError("No particles in analysis region")
+
+    rho_ana, vx_ana, P_ana = np.zeros(len(x_f)), np.zeros(len(x_f)), np.zeros(len(x_f))
+    for i, xi in enumerate(x_f):
+        rho_ana[i], vx_ana[i], P_ana[i] = sod.get_value(t, xi)
+
+    err_rho = np.sqrt(np.mean((rho_f - rho_ana) ** 2)) / np.mean(rho_ana)
+    err_vx = np.sqrt(np.mean((vx_f - vx_ana) ** 2)) / (np.mean(np.abs(vx_ana)) + 0.1)
+    err_vy = np.sqrt(np.mean(vy_f**2))
+    err_vz = np.sqrt(np.mean(vz_f**2))
+    err_P = np.sqrt(np.mean((P_f - P_ana) ** 2)) / np.mean(P_ana)
+    return err_rho, (err_vx, err_vy, err_vz), err_P
+
+
+if shamrock.sys.world_rank() == 0:
+    rho, v, P = compute_L2_errors(data, sod, t_target, -0.5, 0.5)
+    vx, vy, vz = v
+
+    print("current errors :")
+    print(f"err_rho = {rho}")
+    print(f"err_vx = {vx}")
+    print(f"err_vy = {vy}")
+    print(f"err_vz = {vz}")
+    print(f"err_P = {P}")
+
+    # Expected L2 error values (calibrated from local run with M4 kernel)
+    # Tolerance set very strict for regression testing (like sod_tube_sph.py)
+    expect_rho = 0.029892771160040497
+    expect_vx = 0.10118608617971991
+    expect_vy = 0.006382105147197806
+    expect_vz = 3.118241304703099e-05
+    expect_P = 0.038072557056294656
+
+    tol = 1e-8
+
+    test_pass = True
+    err_log = ""
+
+    error_checks = {
+        "rho": (rho, expect_rho),
+        "vx": (vx, expect_vx),
+        "vy": (vy, expect_vy),
+        "vz": (vz, expect_vz),
+        "P": (P, expect_P),
+    }
+
+    for name, (value, expected) in error_checks.items():
+        if abs(value - expected) > tol * expected:
+            err_log += f"error on {name} is outside of tolerances:\n"
+            err_log += f"  expected error = {expected} +- {tol*expected}\n"
+            err_log += (
+                f"  obtained error = {value} (relative error = {(value - expected)/expected})\n"
+            )
+            test_pass = False
+
+    if test_pass:
+        print("\n" + "=" * 50)
+        print("GSPH Sod Shock Tube Test: PASSED")
+        print("=" * 50)

src/shambase/include/shambase/constants.hpp

@@ -11,6 +11,7 @@
 
 /**
  * @file constants.hpp
+ * @author Guo Yansong (guo.yansong.ngy@gmail.com)
  * @author Timothée David--Cléris (tim.shamrock@proton.me)
  * @brief Class holding the value of numerous constants
  * generated from the following source
@@ -40,6 +41,7 @@
         add("gamma_5_6", gamma(5/6))
         add("gamma_1", gamma(1))
         add("sqrt_2", 2**(1/2))
+        add("sqrt_pi", pi**(1/2))
         add("e", e)
   */
 
@@ -65,6 +67,7 @@ namespace shambase::constants {
     template<class T> constexpr T gamma_5_6 = 1.1287870299081257386;
     template<class T> constexpr T gamma_1 = 1;
     template<class T> constexpr T sqrt_2 = 1.4142135623730951455;
+    template<class T> constexpr T sqrt_pi = 1.7724538509055158819; // same as gamma_1_2
     template<class T> constexpr T e = 2.7182818284590450908;
     // clang-format on
 

src/shammodels/gsph/CMakeLists.txt

@@ -11,11 +11,12 @@ cmake_minimum_required(VERSION 3.9)
 
 project(Shammodels_gsph CXX C)
 
-# Sources: GSPH solver, Model, and VTK I/O
+# Sources: GSPH solver, Model, VTK I/O, and Python bindings
 set(Sources
     src/SolverConfig.cpp
     src/Solver.cpp
     src/Model.cpp
+    src/pyGSPHModel.cpp
     src/modules/UpdateDerivs.cpp
     src/modules/io/VTKDump.cpp
 )

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

@@ -130,6 +130,17 @@ namespace shammodels::gsph {
         void compute_eos_fields();
         void reset_eos_fields();
 
+        /**
+         * @brief Compute gradients for MUSCL reconstruction
+         *
+         * Computes density, pressure, and velocity gradients for each particle
+         * using SPH kernel gradient summation. Only called when MUSCL reconstruction
+         * is enabled (reconstruct_config.is_muscl() == true).
+         *
+         * Reference: Cha & Whitworth (2003)
+         */
+        void compute_gradients();
+
         void prepare_corrector();
 
         /**