Address code review feedback

- iterative.hpp: Add epsilon to Z_L/Z_R denominators to prevent division by zero
- Solver.cpp: Add explicit error for unsupported BC types (ShearingPeriodic)

Author: Guo-astro

src/shammodels/gsph/include/shammodels/gsph/math/riemann/iterative.hpp

@@ -122,12 +122,13 @@ namespace shammodels::gsph::riemann {
 
             // Derivatives dW/dp for Newton-Raphson
             // Z_L = -dW_L/dp * W_L (note the sign convention)
+            // Add smallp to denominator to prevent division by zero near vacuum/shock
             Tscal Z_L = Tscal{4} * V_L * W_L * W_L;
-            Z_L       = -Z_L * W_L / (Z_L - gp1 * (p_star - p_L));
+            Z_L       = -Z_L * W_L / (Z_L - gp1 * (p_star - p_L) + smallp);
 
             // Z_R = dW_R/dp * W_R
             Tscal Z_R = Tscal{4} * V_R * W_R * W_R;
-            Z_R       = Z_R * W_R / (Z_R - gp1 * (p_star - p_R));
+            Z_R       = Z_R * W_R / (Z_R - gp1 * (p_star - p_R) + smallp);
 
             // Intermediate velocities from each side
             // u*_L = u_L - (p* - p_L) / W_L

src/shammodels/gsph/src/Solver.cpp

@@ -132,6 +132,10 @@ void shammodels::gsph::Solver<Tvec, Kern>::gen_ghost_handler(Tscal time_val) {
         SolverBCPeriodic *c
         = std::get_if<SolverBCPeriodic>(&solver_config.boundary_config.config)) {
         storage.ghost_handler.set(GhostHandle{scheduler(), BCPeriodic{}, storage.patch_rank_owner});
+    } else {
+        // ShearingPeriodic and other BC types not yet supported in GSPH
+        shambase::throw_with_loc<std::runtime_error>(
+            "GSPH: Unsupported boundary condition type. Only Free and Periodic are supported.");
     }
 }
 
@@ -481,6 +485,10 @@ void shammodels::gsph::Solver<Tvec, Kern>::apply_position_boundary(Tscal time_va
         SolverBCPeriodic *c
         = std::get_if<SolverBCPeriodic>(&solver_config.boundary_config.config)) {
         integrators.fields_apply_periodicity(ixyz, std::pair{bmin, bmax});
+    } else {
+        // ShearingPeriodic and other BC types not yet supported in GSPH
+        shambase::throw_with_loc<std::runtime_error>(
+            "GSPH: Unsupported boundary condition type. Only Free and Periodic are supported.");
     }
 
     reatrib.reatribute_patch_objects(storage.serial_patch_tree.get(), "xyz");