Refactor to eliminate code duplication (SSOT/DRY)

Python scripts:
- Delete duplicate run_gsph_sod_tube.py (same as gsph_sod_shock_tube.py)
- Update animate_sod_tube.py to use shamrock.phys.SodTube instead of custom riemann.py
- Replace riemann.py with sedov.py (keep only SedovAnalytical, use shamphys.SodTube for Sod)

C++ VTK dump utilities:
- Create shared VTKDumpUtils.hpp in shammodels/common/io/
- Refactor both SPH and GSPH VTKDump.cpp to use shared utilities
- Eliminates 5 duplicated helper functions (start_dump, vtk_dump_add_patch_id, etc.)

Net reduction: ~528 lines of duplicate code removed

Author: Guo-astro

exemples/common/analytical/riemann.py

@@ -0,0 +1,131 @@
+#!/usr/bin/env python3
+"""
+Sedov-Taylor Blast Wave Analytical Solution
+
+Computes the self-similar solution for the Sedov-Taylor blast wave problem.
+
+References:
+- Sedov, L.I. (1959) "Similarity and Dimensional Methods in Mechanics"
+- Taylor, G.I. (1950) "The Formation of a Blast Wave by a Very Intense Explosion"
+
+Note: For Sod shock tube, use shamrock.phys.SodTube instead.
+"""
+
+import numpy as np
+
+
+class SedovAnalytical:
+    """
+    Sedov-Taylor blast wave analytical solution.
+
+    Parameters:
+    -----------
+    gamma : float
+        Adiabatic index
+    E_blast : float
+        Total blast energy
+    rho_0 : float
+        Background density
+    ndim : int
+        Dimensionality (1=planar, 2=cylindrical, 3=spherical)
+    """
+
+    def __init__(self, gamma=5.0 / 3.0, E_blast=1.0, rho_0=1.0, ndim=3):
+        self.gamma = gamma
+        self.E_blast = E_blast
+        self.rho_0 = rho_0
+        self.ndim = ndim
+
+        # Similarity exponent
+        self.alpha = 2.0 / (ndim + 2)
+
+        # Sedov constant (depends on gamma and ndim)
+        self.xi_0 = self._compute_xi0()
+
+    def _compute_xi0(self):
+        """Compute the Sedov constant xi_0."""
+        gamma = self.gamma
+        ndim = self.ndim
+
+        # Approximate values for common cases
+        if ndim == 3:  # Spherical
+            if abs(gamma - 5.0 / 3.0) < 0.01:
+                return 1.15167
+            elif abs(gamma - 1.4) < 0.01:
+                return 1.03275
+        elif ndim == 2:  # Cylindrical
+            if abs(gamma - 1.4) < 0.01:
+                return 1.033
+        elif ndim == 1:  # Planar
+            if abs(gamma - 1.4) < 0.01:
+                return 0.911
+
+        # General approximation
+        return 1.0
+
+    def shock_radius(self, t):
+        """Compute shock radius at time t."""
+        if t <= 0:
+            return 0.0
+        return self.xi_0 * (self.E_blast * t**2 / self.rho_0) ** (1.0 / (self.ndim + 2))
+
+    def shock_velocity(self, t):
+        """Compute shock velocity at time t."""
+        if t <= 0:
+            return 0.0
+        R_s = self.shock_radius(t)
+        return 2.0 / (self.ndim + 2) * R_s / t
+
+    def post_shock_density(self):
+        """Compute post-shock density."""
+        return self.rho_0 * (self.gamma + 1) / (self.gamma - 1)
+
+    def solution_at_time(self, t, r_max=None, n_points=500):
+        """
+        Compute the radial profile at time t.
+
+        Returns:
+        --------
+        r, rho, v, p : arrays
+            Radius and primitive variables
+        """
+        if t <= 1e-10:
+            r = np.linspace(0, r_max or 0.01, n_points)
+            return (
+                r,
+                np.ones(n_points) * self.rho_0,
+                np.zeros(n_points),
+                np.ones(n_points) * 1e-10,
+            )
+
+        gamma = self.gamma
+        ndim = self.ndim
+
+        R_s = self.shock_radius(t)
+        v_s = self.shock_velocity(t)
+
+        if r_max is None:
+            r_max = R_s * 1.5
+
+        # Similarity variable lambda = r/R_s
+        lam = np.linspace(0, min(1.0, r_max / R_s), n_points)
+        r = lam * R_s
+
+        # Post-shock values
+        rho_s = self.post_shock_density()
+        v_shock = 2.0 / (gamma + 1) * v_s
+        p_s = 2.0 / (gamma + 1) * self.rho_0 * v_s**2
+
+        # Approximate profiles (self-similar structure)
+        # Velocity: linear profile
+        v = v_shock * lam
+
+        # Density: peaks near shock, low at center
+        omega = (ndim + 2) * gamma / (2 + ndim * (gamma - 1))
+        rho = rho_s * lam ** (omega - 1) * np.maximum(0.1, 1 - 0.8 * (1 - lam) ** 2)
+        rho[0] = rho[1] if n_points > 1 else rho_s * 0.1
+
+        # Pressure: higher at center
+        p = p_s * (0.5 + 0.5 * lam**2)
+
+        return r, rho, v, p

exemples/common/analytical/sedov.py

@@ -24,7 +24,7 @@
 # Import from shared analytical module
 exemples_dir = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
 sys.path.insert(0, exemples_dir)
-from common.analytical.riemann import SedovAnalytical
+from common.analytical.sedov import SedovAnalytical
 
 # Try to import animation tools
 try:
@@ -184,7 +184,7 @@ def compute_radial_profiles(data, n_bins=100):
 print()
 
 # Create analytical solution object
-sedov_analytical = SedovAnalytical(gamma=gamma, E_blast=E_blast, rho_0=rho_0, ndim=3)
+sedov_analytical = SedovAnalytical(gamma=gamma, E_blast=E_blast, rho_0=rho_0)
 
 # Determine frame skip for reasonable animation size
 n_frames = len(files)

exemples/gsph/run_gsph_sod_tube.py

@@ -21,10 +21,8 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-# Import from shared analytical module
-exemples_dir = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
-sys.path.insert(0, exemples_dir)
-from common.analytical.riemann import SodAnalytical
+# Use shamrock's built-in SodTube analytical solution
+import shamrock
 
 # Try to import animation tools
 try:
@@ -135,8 +133,22 @@ def find_snapshots(data_dir):
 print(f"  Right: rho = {rho_R}, P = {p_R}")
 print()
 
-# Create analytical solution object
-sod_analytical = SodAnalytical(gamma=gamma, rho_L=rho_L, rho_R=rho_R, p_L=p_L, p_R=p_R, x0=0.0)
+# Create analytical solution object using shamrock's built-in SodTube
+# Note: SodTube uses rho_1/P_1 for left state and rho_5/P_5 for right state
+sod_analytical = shamrock.phys.SodTube(
+    gamma=gamma, rho_1=rho_L, P_1=p_L, rho_5=rho_R, P_5=p_R
+)
+
+
+def get_analytical_solution(sod, t, x_array):
+    """Get analytical solution at multiple x positions."""
+    rho = np.zeros(len(x_array))
+    vel = np.zeros(len(x_array))
+    pres = np.zeros(len(x_array))
+    for i, x in enumerate(x_array):
+        rho[i], vel[i], pres[i] = sod.get_value(t, x)
+    ene = pres / ((gamma - 1) * np.maximum(rho, 1e-10))
+    return rho, vel, pres, ene
 
 # Determine frame skip for reasonable animation size
 n_frames = len(files)
@@ -203,9 +215,8 @@ def update(frame_num):
     ene_sim = data["ene"][sort_idx]
 
     # Get analytical solution
-    x_ana, rho_ana, vel_ana, pres_ana, ene_ana = sod_analytical.solution_at_time(
-        time, x_min=x_sim.min(), x_max=x_sim.max(), n_points=500
-    )
+    x_ana = np.linspace(x_sim.min(), x_sim.max(), 500)
+    rho_ana, vel_ana, pres_ana, ene_ana = get_analytical_solution(sod_analytical, time, x_ana)
 
     # Clear axes
     ax1.clear()
@@ -309,9 +320,8 @@ def update(frame_num):
     pres_sim = data["pres"][sort_idx]
     ene_sim = data["ene"][sort_idx]
 
-    x_ana, rho_ana, vel_ana, pres_ana, ene_ana = sod_analytical.solution_at_time(
-        time, x_min=x_sim.min(), x_max=x_sim.max(), n_points=500
-    )
+    x_ana = np.linspace(x_sim.min(), x_sim.max(), 500)
+    rho_ana, vel_ana, pres_ana, ene_ana = get_analytical_solution(sod_analytical, time, x_ana)
 
     fig2, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 12))
     fig2.suptitle(