added LAST_VALID_SOLUTION convergence enum

docs/source/examples/shock/example7.rst

@@ -112,6 +112,9 @@ Initialize arrays for storing the solution variables:
     mole_fractions_refl = {}
 
 Loop over the initial shock velocities and solve the shock problem, noting again that we set `reflected=True` to compute the equilibrium reflected shock solution.
+This example still gates downstream use on ``solution.converged``. If a solve
+returns ``solution.last_error == cea.LAST_VALID_SOLUTION``, the retained state
+is still not a converged shock point.
 
 .. code-block:: python
     :emphasize-lines: 4

docs/source/interfaces/c_api.rst

@@ -11,6 +11,9 @@ Notes:
 - ``cea_rocket_solver_get_eqsolver`` is no longer exposed; use ``cea_rocket_solver_get_size`` and other RocketSolver APIs.
 - Rocket, shock, and detonation solution property getters now require a ``len`` argument and will return
   ``CEA_INVALID_SIZE`` if ``len`` is smaller than the internal number of points.
+- ``cea_shock_solver_solve`` may return ``CEA_LAST_VALID_SOLUTION`` when the incident-equilibrium solver retains
+  a last valid state for inspection without converging the shock iteration. In that case the returned properties
+  remain readable, but ``cea_shock_solution_get_converged`` still reports false.
 
 .. doxygenfile:: cea.h
    :project: cea

docs/source/interfaces/python_api.rst

@@ -49,6 +49,12 @@ RocketSolution
 ShockSolver
 -----------
 
+Shock solves can emit a ``RuntimeWarning`` with ``LAST_VALID_SOLUTION`` when the
+incident-equilibrium path retains a last valid state without converging the
+shock iteration. In that case ``ShockSolution.last_error`` records the soft
+failure, while ``ShockSolution.converged`` remains the authoritative strict
+convergence signal.
+
 .. autoclass:: cea.ShockSolver
    :members:
 

source/bind/c/CMakeLists.txt

@@ -97,4 +97,12 @@ if (CEA_BUILD_TESTING)
         WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
     )
 
+    add_executable(cea_bindc_shock_last_valid samples/shock_last_valid.c)
+    target_link_libraries(cea_bindc_shock_last_valid PRIVATE cea::bindc)
+    add_test(
+        NAME cea_bindc_shock_last_valid
+        COMMAND cea_bindc_shock_last_valid
+        WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+    )
+
 endif()

source/bind/c/bindc.F90

@@ -7,6 +7,7 @@ module cea_bindc
     use cea_param, only: empty_dp, gas_constant, get_data_search_dirs
     use cea_input, only: ReactantInput
     use cea_mixture, only: names_match
+    use cea_shock, only: SHOCK_SOLVE_STATUS_LAST_VALID
     use iso_c_binding
     use fb_logging
     use fb_utils, only: assert, locate, is_empty, to_str
@@ -203,6 +204,7 @@ module cea_bindc
         enumerator :: CEA_INVALID_INDEX = 6
         enumerator :: CEA_INVALID_SIZE = 7
         enumerator :: CEA_NOT_CONVERGED = 8
+        enumerator :: CEA_LAST_VALID_SOLUTION = 9
     end enum
 
     !-----------------------------------------------------------------
@@ -2495,7 +2497,16 @@ function cea_shock_solver_solve(sptr, slptr, weights, T0, p0, mach1_or_u1, use_m
             solution = solver%solve(weights(:nr), T0, p0, u1=mach1_or_u1, &
                 reflected=reflected, incident_frozen=incident_frozen, reflected_frozen=reflected_frozen)
         end if
-        if (.not. solution%converged) ierr = CEA_NOT_CONVERGED
+        if (.not. solution%converged) then
+            ! A retained last-valid incident state is intentionally exposed as a
+            ! descriptive soft failure, but it remains a non-converged shock
+            ! solution and callers must check the converged flag separately.
+            if (solution%solve_status == SHOCK_SOLVE_STATUS_LAST_VALID) then
+                ierr = CEA_LAST_VALID_SOLUTION
+            else
+                ierr = CEA_NOT_CONVERGED
+            end if
+        end if
     end function
 
 

source/bind/c/cea.h

@@ -639,6 +639,9 @@ extern "C"
       cea_int *value);
 
   // Solve
+  // Returns CEA_LAST_VALID_SOLUTION when the incident-equilibrium shock path
+  // retains a last valid state for inspection without converging the shock
+  // iteration. In that case, cea_shock_solution_get_converged still reports 0.
   cea_err cea_shock_solver_solve(
       const cea_shock_solver solver,
       cea_shock_solution soln,

source/bind/c/cea_enum.h

@@ -189,4 +189,5 @@
     CEA_INVALID_EQUILIBRIUM_SIZE_TYPE = 5, \
     CEA_INVALID_INDEX            = 6, \
     CEA_INVALID_SIZE             = 7, \
-    CEA_NOT_CONVERGED            = 8
+    CEA_NOT_CONVERGED            = 8, \
+    CEA_LAST_VALID_SOLUTION      = 9

source/bind/c/samples/shock_last_valid.c

@@ -0,0 +1,143 @@
+#include "stdlib.h"
+#include "stdio.h"
+#include "cea.h"
+
+#define LEN(x) (sizeof(x) / sizeof((x)[0]))
+#define NUM_PTS 2
+#define NUM_CANDIDATES 12
+
+int main(void) {
+    const cea_string reactants[] = {"N2 ", "CH4"};
+    const cea_real moles[] = {0.985, 0.015};
+    const cea_string omitted_products[] = {""};
+    const cea_int nr = LEN(reactants);
+    const cea_real p0 = 0.0002986421052631579;
+    const cea_real T0 = 291.0;
+    const cea_real u1_candidates[NUM_CANDIDATES] = {
+        1200.0, 1300.0, 1400.0, 1500.0, 1600.0, 1700.0,
+        1800.0, 1900.0, 2000.0, 2100.0, 2200.0, 2300.0
+    };
+
+    int status = EXIT_FAILURE;
+    int found_last_valid = 0;
+    cea_err ierr;
+    cea_mixture reac = NULL;
+    cea_mixture prod = NULL;
+    cea_shock_solver solver = NULL;
+    cea_shock_solution soln = NULL;
+    cea_real *weights = NULL;
+    cea_real pressure[NUM_PTS];
+    cea_real temperature[NUM_PTS];
+    int converged = 1;
+    size_t i;
+    cea_solver_opts opts;
+
+    cea_set_log_level(CEA_LOG_NONE);
+    ierr = cea_init();
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "cea_init failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    ierr = cea_mixture_create_w_ions(&reac, nr, reactants);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "reactant create failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    ierr = cea_mixture_create_from_reactants_w_ions(&prod, nr, reactants, 0, omitted_products);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "product create failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    ierr = cea_solver_opts_init(&opts);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "opts init failed: %d\n", ierr);
+        goto cleanup;
+    }
+    opts.trace = 1.0e-10;
+    opts.ions = true;
+    opts.transport = true;
+    opts.reactants = reac;
+
+    ierr = cea_shock_solver_create_with_options(&solver, prod, opts);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "solver create failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    ierr = cea_shock_solution_create(&soln, NUM_PTS);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "solution create failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    weights = (cea_real *) calloc(nr, sizeof(cea_real));
+    if (weights == NULL) {
+        fprintf(stderr, "weights allocation failed\n");
+        goto cleanup;
+    }
+
+    ierr = cea_mixture_moles_to_weights(reac, nr, moles, weights);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "weights conversion failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    for (i = 0; i < NUM_CANDIDATES; ++i) {
+        ierr = cea_shock_solver_solve(
+            solver, soln, weights, T0, p0, u1_candidates[i], false, false, false, false
+        );
+        if (ierr == CEA_LAST_VALID_SOLUTION) {
+            found_last_valid = 1;
+            break;
+        }
+        if (ierr != CEA_SUCCESS && ierr != CEA_NOT_CONVERGED) {
+            fprintf(stderr, "unexpected solve status for u1=%g: %d\n", u1_candidates[i], ierr);
+            goto cleanup;
+        }
+    }
+
+    if (!found_last_valid) {
+        fprintf(stderr, "no last-valid incident shock case found across candidate velocities\n");
+        goto cleanup;
+    }
+
+    ierr = cea_shock_solution_get_converged(soln, &converged);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "get_converged failed: %d\n", ierr);
+        goto cleanup;
+    }
+    if (converged != 0) {
+        fprintf(stderr, "expected non-converged shock status\n");
+        goto cleanup;
+    }
+
+    ierr = cea_shock_solution_get_property(soln, CEA_SHOCK_PRESSURE, NUM_PTS, pressure);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "pressure query failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    ierr = cea_shock_solution_get_property(soln, CEA_SHOCK_TEMPERATURE, NUM_PTS, temperature);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "temperature query failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    if (pressure[1] <= 0.0 || temperature[1] <= 0.0) {
+        fprintf(stderr, "expected retained incident state to remain populated\n");
+        goto cleanup;
+    }
+
+    status = EXIT_SUCCESS;
+
+cleanup:
+    if (soln != NULL) cea_shock_solution_destroy(&soln);
+    if (solver != NULL) cea_shock_solver_destroy(&solver);
+    if (prod != NULL) cea_mixture_destroy(&prod);
+    if (reac != NULL) cea_mixture_destroy(&reac);
+    free(weights);
+    return status;
+}

source/bind/python/CEA.pyx

@@ -62,10 +62,19 @@ def _err_name(cea_err ierr):
         return "CEA_INVALID_SIZE"
     if ierr == CEA_NOT_CONVERGED:
         return "CEA_NOT_CONVERGED"
+    if ierr == CEA_LAST_VALID_SOLUTION:
+        return "CEA_LAST_VALID_SOLUTION"
     return f"CEA_ERR_{int(ierr)}"
 
 
 def _check_ierr(cea_err ierr, str context, bint allow_not_converged=True):
+    if ierr == CEA_LAST_VALID_SOLUTION and allow_not_converged:
+        warnings.warn(
+            f"{context}: CEA_LAST_VALID_SOLUTION "
+            f"(last valid shock state retained; solution.converged is False)",
+            RuntimeWarning,
+        )
+        return
     if ierr == CEA_NOT_CONVERGED and allow_not_converged:
         warnings.warn(f"{context}: { _err_name(ierr) }", RuntimeWarning)
         return
@@ -77,6 +86,8 @@ SUCCESS                  = CEA_SUCCESS
 INVALID_FILENAME         = CEA_INVALID_FILENAME
 INVALID_PROPERTY_TYPE    = CEA_INVALID_PROPERTY_TYPE
 INVALID_EQUILIBRIUM_TYPE = CEA_INVALID_EQUILIBRIUM_TYPE
+NOT_CONVERGED           = CEA_NOT_CONVERGED
+LAST_VALID_SOLUTION     = CEA_LAST_VALID_SOLUTION
 
 # Alias the log levels
 LOG_CRITICAL = CEA_LOG_CRITICAL
@@ -3335,6 +3346,12 @@ cdef class ShockSolver:
         reflected_frozen : bool, default False
             Use frozen chemistry for reflected shock
 
+        Notes
+        -----
+        A retained incident last-valid state emits ``RuntimeWarning`` and sets
+        ``soln.last_error`` to ``LAST_VALID_SOLUTION`` while leaving
+        ``soln.converged`` as ``False``.
+
         Raises
         ------
         ValueError
@@ -3921,7 +3938,7 @@ cdef class ShockSolution:
         Returns
         -------
         bool
-            True if solution converged, False otherwise
+            True if the shock iteration converged, False otherwise
         """
         def __get__(self):
             cdef cea_err ierr

source/bind/python/cea/lib/libcea.pyi

@@ -36,6 +36,8 @@ SUCCESS: int
 INVALID_FILENAME: int
 INVALID_PROPERTY_TYPE: int
 INVALID_EQUILIBRIUM_TYPE: int
+NOT_CONVERGED: int
+LAST_VALID_SOLUTION: int
 
 LOG_CRITICAL: int
 LOG_ERROR: int

source/bind/python/cea_def.pxd

@@ -13,6 +13,7 @@ cdef extern from "cea.h":
         CEA_INVALID_INDEX
         CEA_INVALID_SIZE
         CEA_NOT_CONVERGED
+        CEA_LAST_VALID_SOLUTION
 
     ctypedef enum cea_log_level:
         CEA_LOG_CRITICAL

source/bind/python/tests/test_shock_last_valid.py

@@ -0,0 +1,67 @@
+import re
+import warnings
+
+import numpy as np
+import pytest
+
+
+_LAST_VALID_U1_CANDIDATES = (
+    1200.0,
+    1300.0,
+    1400.0,
+    1500.0,
+    1600.0,
+    1700.0,
+    1800.0,
+    1900.0,
+    2000.0,
+    2100.0,
+    2200.0,
+    2300.0,
+)
+
+
+def _find_last_valid_incident_case(cea):
+    reac = cea.Mixture(["N2 ", "CH4"], ions=True)
+    prod = cea.Mixture(["N2 ", "CH4"], products_from_reactants=True, ions=True)
+    solver = cea.ShockSolver(prod, reactants=reac, ions=True, transport=True, trace=1.0e-10)
+
+    weights = reac.moles_to_weights(np.array([0.985, 0.015]))
+    p0 = cea.units.mmhg_to_bar(0.224)
+    T0 = 291.0
+    outcomes = []
+
+    for u1 in _LAST_VALID_U1_CANDIDATES:
+        soln = cea.ShockSolution(solver, reflected=False)
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always")
+            solver.solve(soln, weights, T0, p0, u1=u1, reflected=False)
+        outcomes.append((u1, soln.last_error, soln.converged))
+        if soln.last_error == cea.LAST_VALID_SOLUTION:
+            return soln, caught
+
+    details = ", ".join(
+        f"u1={u1:.0f}: last_error={last_error}, converged={converged}"
+        for u1, last_error, converged in outcomes
+    )
+    pytest.fail(f"No LAST_VALID_SOLUTION case found across candidate velocities. {details}")
+
+
+def test_shock_last_valid_solution_warns_and_preserves_state(cea_module):
+    cea = cea_module
+    soln, caught = _find_last_valid_incident_case(cea)
+
+    pattern = re.compile(
+        r"CEA_LAST_VALID_SOLUTION.*last valid shock state retained.*converged is False"
+    )
+    assert any(
+        isinstance(w.message, RuntimeWarning) and pattern.search(str(w.message))
+        for w in caught
+    )
+
+    assert soln.last_error == cea.LAST_VALID_SOLUTION
+    assert not soln.converged
+    assert soln.P[1] > 0.0
+    assert soln.T[1] > 0.0
+    assert soln.velocity[1] > 0.0
+    assert soln.sonic_velocity[1] > 0.0