Fix/reflected shock w5 (nasa#56)

* initial fix for w5 velocity error

* added transport property output for shock probles

* aligned shock values for failed equilibrium cases

* updated shock problem sequencing

* added cached transport solution for failed solves

* resolved reflected-frozen header difference

* fixed transport property values

* quick cleanup and updated changelog

* fixed v_sonic

* updated reflected velocities; removed pandas dependancy

* added conservation equations test for shock problem

* updated changelog

Author: markleader

CHANGELOG.md

@@ -9,9 +9,15 @@ All notable user-visible changes to this project are documented here.
 - Added SI-focused custom-reactant handling at the Python API layer: `Reactant.temperature` is specified in K and `Reactant.enthalpy` in J/kg (converted internally for core input) (`#53`).
 - Legacy input parsing now supports repeated `outp` dataset keywords (including multiline forms) by merging successive `outp` entries during dataset assembly (`#52`).
 - FAC rocket chamber-closure iteration logic in `RocketSolver_solve_fac` was updated toward CEA2 parity: Option-1 pressure correction direction now follows legacy semantics, the Option-1 convergence check is normalized to assigned injector pressure, the fixed 4-pass outer loop was replaced with tolerance-driven iteration plus a bounded safety guard, and FAC combustor-end reseeding now refreshes from the current infinity state each chamber iteration (`#54`).
+- Legacy SHCK branch sequencing and output selection were reworked toward CEA2 parity: explicit incident requests now report only the requested incident branch (with reflected state attached when requested), while reflected-only requests continue exploring the mixed equilibrium/frozen permutations.
+- Shock reports now include incident/reflected transport tables (`Visc`, `Cp`, `Conductivity`, `Prandtl Number`) when `outp tran` is enabled and transport data are available.
+- Python test/main-interface helper workflows no longer require `pandas`; CSV result output is now written with the standard library and the installation guidance/environment were updated accordingly.
 
 ### Fixed
 - Legacy CLI equilibrium/rocket/shock workflows now propagate `include_ions` into generated product mixtures so ionized products are retained when requested (`#52`).
+- Reflected-shock velocity outputs were corrected for reflected equilibrium/frozen workflows: the legacy CEA2 `u5`/`u5+v2` relations effectively used `rho5/rho2` where the reflected-shock mass balance requires `rho5/rho2 - 1`, so the old equations underpredicted reflected velocities and did not satisfy the reflected conservation equations. CEA now reports conservation-consistent `u5`, `u5+v2`, sonic speed, reflected ratios, and reflected-frozen header thermodynamic values.
+- Shock failure handling now stops at the last valid state instead of emitting partially invalid reflected results; failed downstream shock states are explicitly zeroed for consistent output.
+- Shock transport-property calculations now preserve and restore stable transport seeds across singular-recovery paths, fixing incorrect or missing transport values in difficult reflected-shock cases.
 
 ### Added
 - Added C and Python support for custom reactant data (including species not present in `thermo.lib`) in parity with the main interface workflow used by RP-1311 Example 5 (`#53`).
@@ -20,6 +26,8 @@ All notable user-visible changes to this project are documented here.
   - `cea_mixture_create_products_from_input_reactants_w_ions` (`#53`).
 - Added a shared bindc parser path for `cea_reactant_input -> ReactantInput` conversion to reduce duplicated C-binding logic (`#53`).
 - Added Python `cea.Reactant` and mixed-input `Mixture(...)` support in the Cython binding (`#53`).
+- Added shock regression coverage for transport output population across equilibrium/frozen branches and for singular-recovery reflected-shock cases.
+- Added a Python `pytest` reflected-shock conservation check that verifies incident/reflected mass, momentum, and energy balances and confirms the corrected `u5+v2` relation.
 
 ## [3.1.0] - 2026-03-02
 

docs/source/installation.rst

@@ -238,7 +238,7 @@ Build and install:
 
 3. Prepare Python packages used by the optional binding/tests/docs::
 
-      pip install cython numpy setuptools pandas pytest
+      pip install cython numpy setuptools pytest
 
 4. Configure with the Intel preset (this is the corrected preset usage) and
    build::

environment.yml

@@ -10,7 +10,6 @@ dependencies:
   - gxx
   - ninja
   - numpy
-  - pandas
   - pytest
   - python
   - setuptools

source/bind/python/tests/test_shock_conservation.py

@@ -0,0 +1,94 @@
+import numpy as np
+import pytest
+
+
+_CONSERVATION_RTOL = 5.0e-5
+
+
+def _example7_case(cea):
+    reac_names = ["H2", "O2", "Ar"]
+    moles = np.array([0.05, 0.05, 0.9])
+    p0 = cea.units.mmhg_to_bar(10.0)
+    T0 = 300.0
+
+    reac = cea.Mixture(reac_names)
+    prod = cea.Mixture(reac_names, products_from_reactants=True)
+    solver = cea.ShockSolver(prod, reactants=reac)
+    soln = cea.ShockSolution(solver, reflected=True)
+    weights = reac.moles_to_weights(moles)
+    return solver, soln, weights, T0, p0
+
+
+def _scaled_tolerance(lhs, rhs):
+    return _CONSERVATION_RTOL * max(abs(lhs), abs(rhs), 1.0)
+
+
+def _assert_conserved(lhs, rhs):
+    residual = lhs - rhs
+    assert abs(residual) <= _scaled_tolerance(lhs, rhs)
+
+
+@pytest.mark.parametrize(
+    ("mode", "u1", "kwargs"),
+    [
+        ("equilibrium", 1100.0, {"reflected": True}),
+        ("equilibrium", 1400.0, {"reflected": True}),
+        ("frozen", 1100.0, {"reflected": True, "incident_frozen": True, "reflected_frozen": True}),
+        ("frozen", 1400.0, {"reflected": True, "incident_frozen": True, "reflected_frozen": True}),
+    ],
+)
+def test_shock_conservation_equations(cea_module, mode, u1, kwargs):
+    cea = cea_module
+    solver, soln, weights, T0, p0 = _example7_case(cea)
+
+    solver.solve(soln, weights, T0, p0, u1=u1, **kwargs)
+
+    assert soln.converged, f"{mode} shock solve failed to converge for u1={u1}"
+
+    w1 = soln.velocity[0]
+    w2 = soln.velocity[1]
+    ur = soln.velocity[2]
+    v2 = soln.v2
+
+    rho1 = soln.density[0]
+    rho2 = soln.density[1]
+    rho5 = soln.density[2]
+
+    p1 = soln.P[0] * 1.0e5
+    p2 = soln.P[1] * 1.0e5
+    p5 = soln.P[2] * 1.0e5
+
+    h1 = soln.enthalpy[0]
+    h2 = soln.enthalpy[1]
+    h5 = soln.enthalpy[2]
+
+    # Incident shock conservation in the incident shock frame.
+    incident_mass_lhs = rho1 * w1
+    incident_mass_rhs = rho2 * w2
+    _assert_conserved(incident_mass_lhs, incident_mass_rhs)
+
+    incident_momentum_lhs = p1 + rho1 * w1**2
+    incident_momentum_rhs = p2 + rho2 * w2**2
+    _assert_conserved(incident_momentum_lhs, incident_momentum_rhs)
+
+    incident_energy_lhs = h1 + w1**2 / 2000.0
+    incident_energy_rhs = h2 + w2**2 / 2000.0
+    _assert_conserved(incident_energy_lhs, incident_energy_rhs)
+
+    # Reflected shock conservation in the reflected-shock frame.
+    reflected_upstream_speed = v2 + ur
+    reflected_downstream_speed = ur
+
+    reflected_mass_lhs = rho2 * reflected_upstream_speed
+    reflected_mass_rhs = rho5 * reflected_downstream_speed
+    _assert_conserved(reflected_mass_lhs, reflected_mass_rhs)
+
+    reflected_momentum_lhs = p2 + rho2 * reflected_upstream_speed**2
+    reflected_momentum_rhs = p5 + rho5 * reflected_downstream_speed**2
+    _assert_conserved(reflected_momentum_lhs, reflected_momentum_rhs)
+
+    reflected_energy_lhs = h2 + reflected_upstream_speed**2 / 2000.0
+    reflected_energy_rhs = h5 + reflected_downstream_speed**2 / 2000.0
+    _assert_conserved(reflected_energy_lhs, reflected_energy_rhs)
+
+    assert soln.u5_p_v2 == pytest.approx(soln.v2 + soln.velocity[2], rel=1.0e-10, abs=1.0e-10)