chamber_solver.py

Go to the documentation of this file.

"""Chamber pressure solver: solve supply(Pc) = demand(Pc).
 
This is the authoritative, pure-Python chamber solver used by every general
evaluation path (frontend solves, time-varying, flight, Layer-1 finalization). The
native C kernel does not hook in here: it accelerates only the Layer-1 inner loop,
through a single ``ed_evaluate`` call in
``engine.accel.evaluate``. Keeping the accelerator out of this solver
is deliberate — it leaves the general path simple, authoritative, and free of any
runtime native/Python switch.
"""
 
import numpy as np
import logging
from scipy.optimize import brentq, newton
from typing import Tuple, Dict, Any, List, Optional
 
from engine.pipeline.config_schemas import PintleEngineConfig, ensure_chamber_geometry
from engine.pipeline.combustion_eff import eta_cstar, calculate_Lstar
from engine.pipeline.cea_cache import CEACache
from engine.pipeline.thermal.film_cooling import compute_film_cooling
from engine.pipeline.thermal.regen_cooling import compute_regen_heat_transfer
from engine.pipeline.thermal.ablative_cooling import liner_response
from engine.pipeline.numerical_robustness import (
    PhysicalConstraints,
    NumericalStability,
    PhysicsValidator,
    validate_engine_state,
)
from engine.pipeline.constants import (
    DEFAULT_CHAMBER_TEMP_K,
    DEFAULT_CSTAR_IDEAL_M_S,
    DEFAULT_GAMMA_ND,
    DEFAULT_GAS_CONST_J_KG_K,
    DEFAULT_TURBULENCE_INTENSITY_ND,
)
from engine.core.closure import flows
from engine.core.nozzle import contraction_ratio_of, nozzle_stagnation_loss
 
 
# Chamber-pressure search window. Mirrored VERBATIM in engine/accel/kernels.evaluate_core -- the
# parity gate requires the two root-finds to agree, so change both or neither.
#
# Floor: the throat must be sonic for the demand model mdot = Pc*At/c* to hold, which needs
# Pc > P_amb / (2/(gamma+1))^(gamma/(gamma-1)), i.e. ~1.8 atm across gamma 1.1-1.3. Below that the
# efficiency model collapses, demand blows up, and the residual has a spurious zero crossing near
# 20 psi: with the old 1 bar floor Brent converged to it and reported NEGATIVE thrust for the
# shipped default config at its own configured tank pressures.
PC_CHOKE_FLOOR_PA = 2.0 * 101325.0
# Ceiling: Pc sits below the lower tank pressure by at least this fraction of it -- a floor on the
# total feed + injector drop. It is a search bound, not a design rule; the previous 15% "feed loss
# margin" was a guess that excluded the real operating point of any soft-injector case.
PC_MIN_TOTAL_DROP_FRAC = 0.02
# Equal-step scan resolution used to locate the highest-Pc sign change before Brent.
PC_BRACKET_SCAN_POINTS = 32
 
 
def _cea_table_range(cache) -> dict:
    """The CEA table's extent per axis, {"Pc": [Pa, Pa], "MR": [..], "eps": [..]}; {} if unknown.
 
    ``CEACache.eval`` clamps an off-table input to the edge and says so (``clamped``); a clamped
    point is solved on the edge's thermochemistry, which can be far from its own (c* falls
    ~116 m/s per unit O/F just inside the ethalox table's 2.5 edge and is flat past it).
    """
    out = {}
    for axis in ("Pc", "MR", "eps"):
        try:
            out[axis] = [float(getattr(cache, f"{axis}_min")), float(getattr(cache, f"{axis}_max"))]
        except (AttributeError, TypeError, ValueError):
            continue
    return out
 
 
class ChamberSolver:
    """Solves for chamber pressure by balancing supply and demand"""
    
    def __init__(self, config: PintleEngineConfig, cea_cache: CEACache):
        self.config = config
        self.cea_cache = cea_cache
        
        # Ensure chamber_geometry exists
        ensure_chamber_geometry(config)
        
        # Cache for spray quality (updated during solve)
        self.spray_quality_good = True
    
    def residual(self, Pc: float, P_tank_O: float, P_tank_F: float) -> float:
        """
        Calculate residual: supply(Pc) - demand(Pc)
        
        Residual = 0 when supply equals demand.
        
        Parameters:
        -----------
        Pc : float
            Chamber pressure guess [Pa]
        P_tank_O : float
            Oxidizer tank pressure [Pa]
        P_tank_F : float
            Fuel tank pressure [Pa]
        
        Returns:
        --------
        residual : float [kg/s]
        """
        # Validate inputs
        Pc_val = float(Pc)
        P_tank_O_val = float(P_tank_O)
        P_tank_F_val = float(P_tank_F)
        
        # Physical constraint checks
        if not np.isfinite(Pc_val) or Pc_val <= 0:
            return np.nan  # Invalid pressure
        if not np.isfinite(P_tank_O_val) or P_tank_O_val <= 0:
            return np.nan
        if not np.isfinite(P_tank_F_val) or P_tank_F_val <= 0:
            return np.nan
        
        # Supply side: mass flow from injector (via closure)
        # flows() takes TANK PRESSURES and solves for mdot
        # It internally calculates:
        #   1. Feed losses: P_tank → P_injector
        #   2. Injector flow: P_injector - Pc → mdot
        #   3. Spray constraints: validates and adjusts if needed
        try:
            mdot_O, mdot_F, diagnostics = flows(
                P_tank_O_val,  # Tank pressure (INPUT)
                P_tank_F_val,  # Tank pressure (INPUT)
                Pc_val,        # Chamber pressure (GUESS - being solved for)
                self.config
            )
        except Exception as e:
            # If flows() fails, return NaN to signal invalid point
            import warnings
            warnings.warn(f"flows() failed in residual at Pc={Pc_val/1e6:.2f} MPa: {e}")
            return np.nan
        
        # Validate mass flows
        if not (np.isfinite(mdot_O) and np.isfinite(mdot_F)):
            return np.nan
        
        mdot_supply = mdot_O + mdot_F
        
        # Update spray quality for efficiency calculation
        self.spray_quality_good = diagnostics.get("constraints_satisfied", True)
        
        # Demand side: mass flow required by combustion
        MR, mr_valid = NumericalStability.safe_divide(mdot_O, mdot_F, 2.5, "MR")
        if not mr_valid.passed:
            return np.nan
        
        # Get CEA properties (IDEAL - infinite area equilibrium)
        # For 3D cache, use default expansion ratio from config (doesn't affect chamber properties much)
        try:
            cg = ensure_chamber_geometry(self.config)
            eps_default = cg.expansion_ratio
            cea_props = self.cea_cache.eval(MR, Pc_val, 101325.0, eps_default)
        except Exception as e:
            # Log the error for debugging but return NaN to signal failure
            import warnings
            warnings.warn(f"CEA cache eval failed in residual: {e}")
            return np.nan
        
        cstar_ideal = cea_props.get("cstar_ideal", 0.0)
        if not np.isfinite(cstar_ideal) or cstar_ideal <= 0:
            return np.nan
        
        # Apply combustion efficiency for FINITE CHAMBER
        # This corrects CEA's infinite-area assumption
 
 
        # Apply combustion efficiency for FINITE CHAMBER
        # This corrects CEA's infinite-area assumption
        
        # Calculate current Lstar and injector diameter dynamically
        # This is critical for optimization where config is mutated in-place
        cg = ensure_chamber_geometry(self.config)
        current_Lstar = calculate_Lstar(
            cg.volume,
            cg.A_throat,
            Lstar_override=cg.Lstar
        )
        current_injector_diameter = self._infer_injector_diameter()
 
 
        cooling_results, cooling_eff, _ = self._evaluate_cooling_models(
            Pc_val,
            mdot_O,
            mdot_F,
            cea_props,
            diagnostics,
            with_profile=False,  # root-find: profile is display-only, never needed here
        )
 
        geometry = self._get_chamber_geometry()
        advanced_params = {
            "Pc": Pc_val,
            "Tc": cea_props.get("Tc", DEFAULT_CHAMBER_TEMP_K),
            "cstar_ideal": cea_props.get("cstar_ideal", DEFAULT_CSTAR_IDEAL_M_S),
            "gamma": cea_props.get("gamma", DEFAULT_GAMMA_ND),
            "R": cea_props.get("R", DEFAULT_GAS_CONST_J_KG_K),
            "MR": MR,
            "Ac": geometry["area_cross"],
            "At": cg.A_throat,  # Added At for residence time calculation
            "chamber_length": geometry["length"],  # Added for mixing models
            "Dinj": current_injector_diameter,
            "m_dot_total": mdot_supply,
            "spray_diagnostics": diagnostics,
            "turbulence_intensity": diagnostics.get("turbulence_intensity_mix", DEFAULT_TURBULENCE_INTENSITY_ND),
            "momentum_ratio_R": diagnostics.get("momentum_ratio_R"),
            "R_opt": self._rupe_R_opt(),
            "fuel_props": self._get_fuel_props(),
            "ox_props": self._get_ox_props(),
            "cea_cache": self.cea_cache,  # c*(O/F) for the mixing / vaporized-gas c* ratios
        }
 
        # Add injection velocities to advanced_params if available in diagnostics
        if "u_F" in diagnostics:
            advanced_params["u_fuel"] = diagnostics["u_F"]
        if "u_O" in diagnostics:
            advanced_params["u_lox"] = diagnostics["u_O"]
        
        if hasattr(self, '_debug') and self._debug:
            logging.getLogger("evaluate").info(f"[SOLVER_DEBUG] Pc Guess: {Pc_val/1e6:.4f} MPa | Supply mdot: {mdot_supply:.4f} kg/s | MR: {MR:.3f}")
        
        # Calculate efficiency using advanced physics-based model
        eta = eta_cstar(
            current_Lstar,
            self.config.combustion.efficiency,
            cooling_eff,
            advanced_params,
            debug=self._debug if hasattr(self, '_debug') else False,
        )
        
        # Validate efficiency. No upper bound: eta_c* is against c*_ideal at the BULK O/F, and off
        # the c* peak a spread of stream tubes, or a vaporized gas nearer the peak, can beat it
        # (combustion_physics.stream_tube_mixing_efficiency / calculate_vaporization_efficiency).
        if not np.isfinite(eta) or eta <= 0:
            return np.nan
 
        
        # Actual c* accounting for finite chamber volume
        cstar_actual = eta * cstar_ideal
        
        # Demand: mdot = P0 * At / c*_actual, P0 = Pc / kappa the nozzle stagnation pressure
        # (Pc is the injector-end pressure; kappa the finite-area combustor's Rayleigh loss).
        cg = ensure_chamber_geometry(self.config)
        kappa = nozzle_stagnation_loss(contraction_ratio_of(cg), float(cea_props.get("gamma", DEFAULT_GAMMA_ND)))
        mdot_demand, demand_valid = NumericalStability.safe_divide(
            Pc_val / kappa * cg.A_throat,
            cstar_actual,
            0.0,
            "mdot_demand"
        )
        if not demand_valid.passed:
            return np.nan
        
        residual = mdot_supply - mdot_demand
        
        # Validate residual is finite
        if not np.isfinite(residual):
            return np.nan
        
        return float(residual)
 
    def _accel_chamber_pc(self, P_tank_O: float, P_tank_F: float):
        """Solve Pc with the native C kernel (whole residual loop + Brent) when native
        is enabled and can handle this config; otherwise return None so the Python
        Brent solver runs.
 
        Pure capability routing — no runtime parity self-check or process latch. Parity
        is enforced by the golden test suite and the load-time ABI assert in
        the accelerator. A solve that raises or returns a non-finite Pc falls back to the
        Python solver for that call.
        """
        from engine import accel
        if not accel.enabled():
            return None
        try:
            res = accel.chamber_solve(self.config, self.cea_cache, P_tank_O, P_tank_F)
        except Exception:
            # Strict mode (ED_REQUIRE_ACCEL=1, CI parity job): surface a genuine
            # accelerator failure loudly instead of silently falling back to Python
            # (which would report a false green). Default runs fall back quietly.
            if accel.require():
                raise
            return None
        if res is None:
            return None
        Pc_native = float(res[0])
        if not np.isfinite(Pc_native) or Pc_native <= 0.0:
            return None
        return Pc_native
 
    @staticmethod
    def _highest_sign_change(f, Pc_min: float, Pc_max: float, r_min: float, r_max: float,
                             n: int = PC_BRACKET_SCAN_POINTS):
        """[lo, hi] holding the HIGHEST-Pc sign change of the residual on [Pc_min, Pc_max], or
        (None, None). Equal steps scanned from the top; ``r_min``/``r_max`` are reused at the ends.
        Verbatim mirror of the scan in accel.kernels.evaluate_core (parity)."""
        pb, rb = Pc_max, r_max
        for i in range(n - 1, -1, -1):
            pa = Pc_min + (Pc_max - Pc_min) * (i / n)
            ra = r_min if i == 0 else f(pa)
            if np.isfinite(ra) and np.isfinite(rb) and np.sign(ra) != np.sign(rb):
                return pa, pb
            pb, rb = pa, ra
        return None, None
 
    def _raise_supply_exceeds_demand(self, P_tank_O, P_tank_F, Pc_max, residual_min, residual_max, debug):
        """Supply > demand even at the ceiling: injector oversized / throat undersized for these
        pressures. Raise with the supply/demand numbers at Pc_max so the message is actionable."""
        try:
            mdot_O_test, mdot_F_test, diag_test = flows(P_tank_O, P_tank_F, Pc_max, self.config)
            mdot_supply_test = mdot_O_test + mdot_F_test
            MR_test = mdot_O_test / mdot_F_test if mdot_F_test > 0 else np.inf
            cg = ensure_chamber_geometry(self.config)
            cea_props_test = self.cea_cache.eval(MR_test, Pc_max, 101325.0, cg.expansion_ratio)
            cstar_ideal_test = cea_props_test.get("cstar_ideal", 0.0)
            geometry_test = self._get_chamber_geometry()
            advanced_params_test = {
                "Pc": Pc_max,
                "Tc": cea_props_test.get("Tc", DEFAULT_CHAMBER_TEMP_K),
                "cstar_ideal": cstar_ideal_test,
                "gamma": cea_props_test.get("gamma", DEFAULT_GAMMA_ND),
                "R": cea_props_test.get("R", DEFAULT_GAS_CONST_J_KG_K),
                "MR": MR_test,
                "Ac": geometry_test["area_cross"],
                "At": cg.A_throat,
                "chamber_length": geometry_test["length"],
                "Dinj": self._infer_injector_diameter(),
                "m_dot_total": mdot_supply_test,
                "spray_diagnostics": diag_test,
                "turbulence_intensity": diag_test.get("turbulence_intensity_mix", DEFAULT_TURBULENCE_INTENSITY_ND),
                "fuel_props": self._get_fuel_props(),
            "ox_props": self._get_ox_props(),
            "cea_cache": self.cea_cache,  # c*(O/F) for the mixing / vaporized-gas c* ratios
            }
            eta_test = eta_cstar(
                calculate_Lstar(cg.volume, cg.A_throat, Lstar_override=cg.Lstar),
                self.config.combustion.efficiency,
                diag_test.get("cooling_efficiency", 1.0),
                advanced_params_test,
                debug=debug,
            )
            cstar_actual_test = eta_test * cstar_ideal_test
            kappa_test = nozzle_stagnation_loss(contraction_ratio_of(cg), float(cea_props_test.get("gamma", DEFAULT_GAMMA_ND)))
            mdot_demand_test = (Pc_max / kappa_test * cg.A_throat) / cstar_actual_test if cstar_actual_test > 0 else np.inf
            if mdot_demand_test > 0 and mdot_supply_test > mdot_demand_test:
                Pc_estimate = mdot_supply_test * cstar_actual_test / cg.A_throat * kappa_test
                raise ValueError(
                    f"No solution: Supply > Demand at all Pc. "
                    f"Residual at Pc_min: {residual_min:.4f} kg/s, at Pc_max: {residual_max:.4f} kg/s. "
                    f"\nDiagnostics at Pc_max ({Pc_max/1e6:.2f} MPa):"
                    f"\n  - Supply: {mdot_supply_test:.4f} kg/s (mdot_O={mdot_O_test:.4f}, mdot_F={mdot_F_test:.4f})"
                    f"\n  - Demand: {mdot_demand_test:.4f} kg/s (c*_actual={cstar_actual_test:.1f} m/s, At={cg.A_throat*1e6:.2f} mm²)"
                    f"\n  - Estimated Pc needed: {Pc_estimate/1e6:.2f} MPa (vs Pc_max={Pc_max/1e6:.2f} MPa)"
                    f"\nPossible fixes:"
                    f"\n  1. Reduce injector orifice areas (currently oversized)"
                    f"\n  2. Increase throat area (currently undersized)"
                    f"\n  3. Increase tank pressures to allow higher Pc_max"
                    f"\n  4. Check combustion efficiency (low efficiency reduces demand)"
                )
            raise ValueError(
                f"No solution: Supply > Demand at all Pc. "
                f"Residual: [{residual_min:.4f}, {residual_max:.4f}] kg/s. "
                f"Could not compute detailed diagnostics."
            )
        except ValueError:
            raise
        except Exception as diag_e:
            raise ValueError(
                f"No solution: Supply > Demand at all Pc. "
                f"Residual at bounds: [{residual_min:.4f}, {residual_max:.4f}] kg/s. "
                f"Pc_max ({Pc_max/1e6:.2f} MPa) limited by tank pressure. "
                f"Possible causes: Injector oversized, throat undersized, or combustion efficiency too low. "
                f"Diagnostic error: {diag_e}"
            )
 
    def solve(
        self,
        P_tank_O: float,
        P_tank_F: float,
        Pc_guess: float = None,
        debug: bool = False,
        silent: bool = False,
    ) -> Tuple[float, Dict[str, Any]]:
        """
        Solve for chamber pressure.
        
        Parameters:
        -----------
        P_tank_O : float
            Oxidizer tank pressure [Pa]
        P_tank_F : float
            Fuel tank pressure [Pa]
        Pc_guess : float, optional
            Initial guess for chamber pressure [Pa]
        
        Returns:
        --------
        Pc : float [Pa]
            Solved chamber pressure
        diagnostics : dict
            Solution diagnostics
        """
        # Search window: choked-flow floor, tank-pressure ceiling (see the constants above), then
        # the user's solver.Pc_bounds narrow it further.
        Pc_min = max(PC_CHOKE_FLOOR_PA, self.config.solver.Pc_bounds[0])
        Pc_max = min(min(P_tank_O, P_tank_F) * (1.0 - PC_MIN_TOTAL_DROP_FRAC), self.config.solver.Pc_bounds[1])
        
        if Pc_max <= Pc_min:
            raise ValueError(f"Invalid pressure bounds: Pc_max ({Pc_max}) <= Pc_min ({Pc_min})")
        
        # Initial guess
        if Pc_guess is None:
            Pc_guess = (Pc_min + Pc_max) / 2
        
        # Create residual function with tank pressures bound
        # Create residual function with tank pressures bound
        # Also pass debug to residual if needed in future (currently residual uses instance state, but we can't easily pass debug to it without changing signature broadly or storing state)
        # However, residual calls eta_cstar which needs debug...
        # Wait, residual() calls eta_cstar() inside.
        # I need to update residual() to use debug flag.
        # I'll store `self._debug = debug` temporarily or modify residual validation.
        # Storing on self is easiest for this scope.
        self._debug = debug
        self._silent = silent  # gates the display-only cooling profile in post-processing
 
        def residual_func(Pc):
            return self.residual(Pc, P_tank_O, P_tank_F)
 
        # Accelerated path: the whole residual loop + Brent in the numba kernel when it can
        # handle this config. Importantly, this path is hit by the ~30% of Layer-1 CMA candidates
        # that don't converge in the single-call seam and fall back to runner.evaluate -> here, so
        # keeping it accelerated keeps that fallback fast (a pure-Python Brent solve here is ~100x
        # slower). Any failure -> Python Brent below.
        convergence_history: list = []
        _accel_pc = self._accel_chamber_pc(P_tank_O, P_tank_F)
        if _accel_pc is not None:
            Pc = _accel_pc
            success = True
        else:
            residual_min = residual_func(Pc_min)
            residual_max = residual_func(Pc_max)
 
            if not np.isfinite(residual_min):
                # Try to diagnose the issue
                try:
                    test_res = residual_func((Pc_min + Pc_max) / 2)
                    if not np.isfinite(test_res):
                        raise ValueError(
                            f"Residual function returns non-finite values. "
                            f"Pc_min={Pc_min/1e6:.2f} MPa, Pc_max={Pc_max/1e6:.2f} MPa. "
                            f"Check injector geometry, feed system, or CEA cache."
                        )
                except Exception as e:
                    raise ValueError(
                        f"Residual function evaluation failed at bounds. "
                        f"Pc_min={Pc_min/1e6:.2f} MPa, Pc_max={Pc_max/1e6:.2f} MPa. "
                        f"Error: {e}"
                    )
            if not np.isfinite(residual_max):
                raise ValueError(
                    f"Residual function returns non-finite at Pc_max={Pc_max/1e6:.2f} MPa. "
                    f"Check that tank pressures are sufficient and injector geometry is valid."
                )
 
            # The residual is not monotonic in Pc, so neither an endpoint sign test nor a Brent
            # over the whole window is safe: the operating point is the HIGHEST-Pc root, and the
            # crossing nearest the floor (if any) is the barely-flowing spurious one. Locate the
            # top-most sign change and bracket Brent inside it.
            lo, hi = self._highest_sign_change(residual_func, Pc_min, Pc_max, residual_min, residual_max)
            if lo is None:
                if residual_min > 0 and residual_max > 0:
                    self._raise_supply_exceeds_demand(P_tank_O, P_tank_F, Pc_max, residual_min, residual_max, debug)
                raise ValueError(
                    f"No solution: Supply < Demand at all Pc. "
                    f"Residual at bounds: [{residual_min:.4f}, {residual_max:.4f}] kg/s. "
                    f"Insufficient mass flow. Check tank pressures and injector geometry."
                )
 
            def tracked_residual_func(Pc):
                res = residual_func(Pc)
                convergence_history.append(float(res))
                return res
 
            try:
                if self.config.solver.method == "brentq":
                    Pc, result = brentq(
                        tracked_residual_func,
                        lo,
                        hi,
                        xtol=self.config.solver.tolerance,
                        rtol=self.config.solver.tolerance * 1e-3,  # Relative tolerance
                        maxiter=self.config.solver.max_iterations,
                        full_output=True
                    )
                    success = result.converged
                else:
                    # Newton from the middle of the bracket (less robust; kept for configs that ask)
                    Pc = newton(
                        tracked_residual_func,
                        0.5 * (lo + hi),
                        tol=self.config.solver.tolerance,
                        maxiter=self.config.solver.max_iterations
                    )
                    success = True
                conv_check = NumericalStability.check_convergence(
                    convergence_history,
                    self.config.solver.tolerance,
                    min_iterations=3
                )
                if not conv_check.passed and conv_check.severity == "error":
                    raise RuntimeError(f"Convergence validation failed: {conv_check.message}")
            except ValueError:
                raise
            except Exception as e:
                raise RuntimeError(f"Chamber pressure solver failed: {e}")
 
        # Validate solution
        Pc_val = float(Pc)
        if not np.isfinite(Pc_val):
            raise RuntimeError(f"Solver returned non-finite pressure: {Pc_val}")
        
        Pc_check = PhysicalConstraints.validate_pressure(Pc_val, "Pc_solution")
        if not Pc_check.passed and Pc_check.severity == "error":
            raise RuntimeError(f"Solution validation failed: {Pc_check.message}")
        
        # Get final diagnostics
        mdot_O, mdot_F, closure_diag = flows(P_tank_O, P_tank_F, Pc_val, self.config)
        MR = mdot_O / mdot_F if mdot_F > 0 else 0
 
        # Calculate current Lstar and injector diameter dynamically
        # This is critical for optimization where config is mutated in-place
        cg = ensure_chamber_geometry(self.config)
        current_Lstar = calculate_Lstar(
            cg.volume,
            cg.A_throat,
            Lstar_override=cg.Lstar
        )
        current_injector_diameter = self._infer_injector_diameter()
        
        # Use current expansion ratio for 3D cache
        # FIXED: Add safety check for division by zero
        cg = ensure_chamber_geometry(self.config)
        eps_current = cg.A_exit / cg.A_throat if cg.A_throat and cg.A_exit and cg.A_throat > 0 else cg.expansion_ratio
        cea_props = self.cea_cache.eval(MR, Pc_val, 101325.0, eps_current)
        
        # Calculate total mass flow rate (needed for various calculations below)
        mdot_total = mdot_O + mdot_F
        
        # Calculate cooling effects early (needed for conservative reaction kinetics).
        # The display-only per-segment heat-flux profile is skipped on the silent
        # (optimizer) path — it's the dominant post-processing cost and doesn't
        # affect any returned scalar.
        cooling_results, cooling_eff, effective_Tc = self._evaluate_cooling_models(
            Pc_val,
            mdot_O,
            mdot_F,
            cea_props,
            closure_diag,
            with_profile=not getattr(self, "_silent", False),
        )
 
        # Extract and validate mixture diagnostics (diagnostics-only, no efficiency impact)
        # Enable mixture coupling diagnostics if configured
        eff_cfg = self.config.combustion.efficiency
        if getattr(eff_cfg, 'use_mixture_coupling', False) and isinstance(closure_diag, dict):
            # Extract mixture diagnostics in non-strict mode (warn but don't fail)
            mixture_diag = self._extract_and_validate_mixture_diagnostics(closure_diag, strict=False)
            # Log diagnostics for informational purposes
            if debug:
                self._log_mixture_diagnostics(mixture_diag)
            # Store in closure_diag for later analysis/logging
            closure_diag['mixture_diagnostics'] = mixture_diag
 
 
        # Build advanced parameters for combustion efficiency calculation
        geometry = self._get_chamber_geometry()
        
        advanced_params = {
            "Pc": Pc_val,
            "Tc": cea_props["Tc"],  # Ideal Tc (Conservative Residence Time)
            "Tc_kinetics": effective_Tc, # Actual Tc (Conservative Kinetics)
            "cstar_ideal": cea_props.get("cstar_ideal", DEFAULT_CSTAR_IDEAL_M_S),
            "gamma": cea_props.get("gamma", DEFAULT_GAMMA_ND),
            "R": cea_props.get("R", DEFAULT_GAS_CONST_J_KG_K),
            "MR": MR,
            "Ac": geometry["area_cross"],
            "At": cg.A_throat,
            "chamber_length": geometry["length"],
            "Dinj": current_injector_diameter,
            "m_dot_total": mdot_total,
            "u_fuel": closure_diag.get("u_F"),
            "u_lox": closure_diag.get("u_O"),
            "spray_diagnostics": closure_diag,
            "turbulence_intensity": closure_diag.get("turbulence_intensity_mix", DEFAULT_TURBULENCE_INTENSITY_ND),
            "momentum_ratio_R": closure_diag.get("momentum_ratio_R"),
            "R_opt": self._rupe_R_opt(),
            "fuel_props": self._get_fuel_props(),
            "ox_props": self._get_ox_props(),
            "cea_cache": self.cea_cache,  # c*(O/F) for the mixing / vaporized-gas c* ratios
        }
        
        eta = eta_cstar(
            current_Lstar,
            self.config.combustion.efficiency,
            cooling_eff,
            advanced_params,
            debug=debug,
        )
        
        # Comprehensive validation of final solution
        cstar_actual = eta * cea_props["cstar_ideal"]
        gamma = cea_props["gamma"]
        R = cea_props["R"]
 
        if debug:
            logging.getLogger("evaluate").info(
                f"[CSTAR] eta={eta:.4f} | c*_ideal={cea_props['cstar_ideal']:.1f} m/s -> c*_actual={cstar_actual:.1f} m/s | "
                f"ratio={cstar_actual/cea_props['cstar_ideal']:.4f}"
            )   
        
        # Calculate Isp for validation
        # CRITICAL FIX: Correct Isp formula
        # Isp = F / (mdot * g0) = (Cf * Pc * At) / (mdot * g0)
        # OR equivalently: Isp = cstar * Cf / g0
        # The previous formula with gamma * sqrt(...) was incorrect
        # Sanity-check Isp on the delivered basis at sea level (this solver has no site Pa).
        g0 = 9.80665
        cg = ensure_chamber_geometry(self.config)
        kappa = nozzle_stagnation_loss(contraction_ratio_of(cg), float(cea_props["gamma"]))
        F_check = (cg.nozzle_efficiency * cea_props["Cf_vac"] * Pc_val / kappa * cg.A_throat
                   - 101325.0 * (cg.A_exit or cg.A_throat * cg.expansion_ratio))
        Isp = F_check / (mdot_total * g0) if mdot_total > 0 else 0.0
        
        # Validate engine state (use effective temperature after cooling)
        validation_results = validate_engine_state(
            Pc_val, MR, mdot_total, cstar_actual, gamma, effective_Tc, Isp, 0.0  # F not calculated yet
        )
        
        # Check for critical errors
        critical_errors = [r for r in validation_results if r.severity == "error" and not r.passed]
        if critical_errors:
            error_msgs = [r.message for r in critical_errors]
            raise RuntimeError(f"Solution validation failed:\n" + "\n".join(error_msgs))
        
        diagnostics = {
            "Pc": Pc_val,
            "mdot_O": mdot_O,
            "mdot_F": mdot_F,
            "mdot_total": mdot_total,
            "MR": MR,
            "cstar_ideal": cea_props["cstar_ideal"],
            "Tc_ideal": cea_props["Tc"],  # Store original ideal temperature
            "cstar_actual": cstar_actual,
            "Pc_ns": Pc_val / kappa,  # nozzle stagnation pressure
            "stagnation_loss_kappa": kappa,
            "cea_extrapolated": bool(cea_props.get("extrapolated", False)),  # converged point off the table
            "cea_clamped": tuple(cea_props.get("clamped", ())),
            # The table those were read against, so a report can say how far off it the point sits.
            "cea_table_range": _cea_table_range(self.cea_cache),
            "eta_cstar": eta,
            "cooling_efficiency": cooling_eff,
            "Tc": effective_Tc,  # Use effective temperature after cooling (accounts for energy removal)
            "Tc_ideal": cea_props["Tc"],  # Store original CEA temperature for reference
            "gamma": gamma,
            "R": R,
            "M": cea_props.get("M"),  # Molecular weight [kg/kmol]
            "spray_quality_good": self.spray_quality_good,
            "validation_results": validation_results,  # Include validation results
            "convergence_history": convergence_history,  # Include convergence history
            **closure_diag,
        }
 
        diagnostics["cooling"] = cooling_results
 
        return Pc_val, diagnostics
 
    def _extract_and_validate_mixture_diagnostics(
        self, 
        closure_diag: Dict[str, Any], 
        strict: bool = True
    ) -> Dict[str, float]:
        """
        Extract and validate mixture quality diagnostics from closure diagnostics.
        
        DIAGNOSTICS-ONLY - no efficiency impact. Physics comes from:
        - eta_Lstar: actual SMD via gasification model
        - eta_mixing: actual turbulence/velocities via diffusion
        - eta_kinetics: actual T/P via Damköhler number
        
        Parameters:
        - closure_diag: Injector diagnostics dict
        - strict: If True, raise on missing critical data. If False, warn.
        
        Returns: Dict with actual_smd_microns, x_star_mm, We_O, We_F, etc.
        """
        if not isinstance(closure_diag, dict):
            if strict:
                raise ValueError("closure_diag must be a dictionary")
            else:
                import warnings
                warnings.warn("closure_diag is not a dictionary, cannot extract mixture diagnostics")
                return {}
        
        diagnostics = {}
        
        # === EXTRACT SMD (CRITICAL) ===
        D32_O = closure_diag.get("D32_O")
        D32_F = closure_diag.get("D32_F")
        
        if D32_O is None and D32_F is None:
            msg = "Missing SMD data: both D32_O and D32_F are None in closure diagnostics"
            if strict:
                raise ValueError(msg)
            else:
                import warnings
                warnings.warn(msg)
        else:
            actual_smd_m = max(
                float(D32_O) if D32_O is not None else 0.0,
                float(D32_F) if D32_F is not None else 0.0,
            )
            actual_smd_microns = actual_smd_m * 1e6
            
            if not np.isfinite(actual_smd_microns):
                msg = (
                    f"Non-finite SMD: {actual_smd_microns:.2f} microns. "
                    f"D32_O={D32_O}, D32_F={D32_F}. Check spray breakup calculations."
                )
                if strict:
                    raise ValueError(msg)
                else:
                    import warnings
                    warnings.warn(msg)
            elif actual_smd_microns <= 0:
                msg = (
                    f"Invalid SMD: {actual_smd_microns:.2f} microns (must be > 0). "
                    f"D32_O={D32_O}, D32_F={D32_F}. Check spray breakup calculations."
                )
                if strict:
                    raise ValueError(msg)
                else:
                    import warnings
                    warnings.warn(msg)
            else:
                diagnostics["actual_smd_microns"] = float(actual_smd_microns)
        
        # === EXTRACT IMPINGEMENT DISTANCE (CRITICAL) ===
        x_star_m = closure_diag.get("x_star")
        
        if x_star_m is None:
            msg = "Missing x_star in closure diagnostics. Check spray impingement calculations."
            if strict:
                raise ValueError(msg)
            else:
                import warnings
                warnings.warn(msg)
        else:
            x_star_mm = float(x_star_m) * 1000.0
            
            if not np.isfinite(x_star_mm):
                msg = f"Non-finite x_star: {x_star_mm} mm. Must be finite."
                if strict:
                    raise ValueError(msg)
                else:
                    import warnings
                    warnings.warn(msg)
            else:
                # x_star <= 0 is valid (upstream/no impingement)
                diagnostics["x_star_mm"] = float(x_star_mm)
        
        # === EXTRACT WEBER NUMBERS (NON-CRITICAL) ===
        We_O = closure_diag.get("We_O")
        We_F = closure_diag.get("We_F")
        
        if We_O is not None:
            we_o = float(We_O)
            if np.isfinite(we_o) and we_o >= 0:
                diagnostics["We_O"] = we_o
            elif strict:
                raise ValueError(f"Invalid We_O: {we_o}. Must be finite and non-negative.")
        
        if We_F is not None:
            we_f = float(We_F)
            if np.isfinite(we_f) and we_f >= 0:
                diagnostics["We_F"] = we_f
            elif strict:
                raise ValueError(f"Invalid We_F: {we_f}. Must be finite and non-negative.")
        
        # Compute We_min if both present
        if "We_O" in diagnostics and "We_F" in diagnostics:
            we_o = diagnostics["We_O"]
            we_f = diagnostics["We_F"]
            
            # Use minimum of positive values, or max if one is zero
            if we_o > 0 and we_f > 0:
                diagnostics["We_min"] = min(we_o, we_f)
            else:
                diagnostics["We_min"] = max(we_o, we_f)
        
        # === EXTRACT TURBULENCE INTENSITY (OPTIONAL) ===
        I_mix = closure_diag.get("turbulence_intensity_mix")
        
        if I_mix is not None:
            I_mix_val = float(I_mix)
            
            if not np.isfinite(I_mix_val):
                msg = f"Non-finite turbulence intensity: {I_mix_val}"
                if strict:
                    raise ValueError(msg)
                else:
                    import warnings
                    warnings.warn(msg)
            elif I_mix_val <= 0 or I_mix_val > 1:
                msg = (
                    f"Turbulence intensity out of range: {I_mix_val:.4f} (expected 0-1). "
                    f"Check turbulence calculations."
                )
                if strict:
                    raise ValueError(msg)
                else:
                    import warnings
                    warnings.warn(msg)
            else:
                diagnostics["turbulence_intensity_mix"] = I_mix_val
        
        return diagnostics
    
    def _log_mixture_diagnostics(self, diag: Dict[str, float]) -> None:
        """
        Log mixture quality diagnostics for informational purposes.
        
        DIAGNOSTICS-ONLY - no efficiency impact.
        """
        if not diag:
            return
        
        parts = []
        
        if "actual_smd_microns" in diag:
            parts.append(f"SMD={diag['actual_smd_microns']:.1f} μm")
        
        if "x_star_mm" in diag:
            parts.append(f"x*={diag['x_star_mm']:.2f} mm")
        
        if "We_O" in diag and "We_F" in diag:
            parts.append(f"We_O={diag['We_O']:.1f}, We_F={diag['We_F']:.1f}")
            if "We_min" in diag:
                parts.append(f"We_min={diag['We_min']:.1f}")
        
        if "turbulence_intensity_mix" in diag:
            parts.append(f"I_mix={diag['turbulence_intensity_mix']:.3f}")
        
        if parts:
            import logging
            logging.getLogger("evaluate").info(f"[MIXTURE_DIAG] {', '.join(parts)}")
 
    def _compute_cooling_efficiency(
        self,
        cooling_results: Dict[str, Any],
        mdot_total: float,
        Tc: float,
        gamma: float,
        R: float,
    ) -> float:
        eff_cfg = self.config.combustion.efficiency
        if not eff_cfg.use_cooling_coupling:
            return 1.0
 
        total_heat_removed = 0.0
        for source in cooling_results.values():
            if isinstance(source, dict):
                total_heat_removed += float(source.get("heat_removed", 0.0))
 
        if total_heat_removed <= 0:
            return 1.0
 
        cp = gamma * R / max(gamma - 1.0, 1e-6)
        available_energy = mdot_total * cp * max(Tc, 1.0)
        if available_energy <= 0:
            return 1.0
 
        factor = 1.0 - total_heat_removed / available_energy
        return float(np.clip(factor, eff_cfg.cooling_efficiency_floor, 1.0))
 
    def _evaluate_cooling_models(
        self,
        Pc: float,
        mdot_O: float,
        mdot_F: float,
        cea_props: Dict[str, float],
        closure_diag: Dict[str, Any],
        with_profile: bool = True,
    ) -> Tuple[Dict[str, Any], float, float]:
        config = self.config
        mdot_total = mdot_O + mdot_F
        cooling_results: Dict[str, Any] = {}
        Tc = float(cea_props["Tc"])
 
        if mdot_total <= 0:
            closure_diag["cooling"] = cooling_results
            return cooling_results, 1.0, Tc
 
        fuel_fluid = config.fluids["fuel"]
        geometry = self._get_chamber_geometry()
        Pc_val = float(Pc)
        gamma = float(cea_props["gamma"])
        R = float(cea_props["R"])
 
        rho_g = max(Pc_val / (R * max(Tc, 1.0)), 1e-6)
        area_cross = geometry["area_cross"]
        velocity_g = mdot_total / (rho_g * area_cross)
 
        regen_cfg = config.regen_cooling
        # A disabled regen block describes hardware that is not there: none of its gas-side
        # numbers (bore, k, Pr, emissivity, turbulence, recovery) may reach the other models.
        regen_active = regen_cfg is not None and regen_cfg.enabled
        MR = mdot_O / mdot_F if mdot_F > 0 else float(cea_props.get("MR", 0.0) or 0.0)
        tr = self.cea_cache.aux.transport(MR, Pc_val, "chamber")
        mu_g_config = regen_cfg.hot_gas_viscosity if regen_active else tr["mu"]
        M = cea_props.get("M")  # Molecular weight [kg/kmol]
        if regen_active and M is not None and M > 0 and Tc > 0:
            from engine.pipeline.thermal.regen_cooling import calculate_gas_viscosity_huzel
            mu_g_calculated = calculate_gas_viscosity_huzel(Tc, M)
        else:
            mu_g_calculated = tr["mu"]
        mu_g = mu_g_calculated
        k_g = regen_cfg.hot_gas_thermal_conductivity if regen_active else tr["k"]
        Pr_g = regen_cfg.hot_gas_prandtl if (regen_active and regen_cfg.hot_gas_prandtl > 0) else tr["Pr"]
 
        Re_g = rho_g * velocity_g * geometry["diameter"] / max(mu_g, 1e-8)
        Nu_g = 0.023 * (Re_g ** 0.8) * (Pr_g ** 0.4)
 
        # Core turbulence of fully developed pipe flow; the film model's erosion input.
        turbulence_intensity_calc = float(np.clip(0.16 * Re_g ** -0.125, 0.02, 0.25)) if Re_g > 0 else 0.05
        turbulence_boost = 1.0
        if regen_active:
            turbulence_boost = 1.0 + max(regen_cfg.gas_turbulence_intensity, 0.0)
            turbulence_intensity_calc = max(
                turbulence_intensity_calc,
                float(np.clip(regen_cfg.gas_turbulence_intensity, 0.0, 0.5)),
            )
 
        h_hot_base = Nu_g * k_g / geometry["diameter"] * turbulence_boost
 
        gas_props = {
            "Pc": Pc_val,
            "Tc": Tc,
            "gamma": gamma,
            "R": R,
            "rho": rho_g,
            "velocity": velocity_g,
            "length": geometry["length"],
            "circumference": geometry["circumference"],
            "area": geometry["area"],
            "area_cross": area_cross,
            "h_hot_base": h_hot_base,
            "turbulence_intensity": turbulence_intensity_calc,
        }
 
        film_cfg = config.film_cooling
        film_results = {
            "enabled": False,
            "mdot_available_for_regen": mdot_F,
            "effective_gas_temperature": Tc,
            "heat_removed": 0.0,
        }
 
        if film_cfg is not None and film_cfg.enabled:
            film_results = compute_film_cooling(
                mdot_total,
                mdot_F,
                gas_props,
                film_cfg,
                fuel_fluid,
            )
            cooling_results["film"] = film_results
 
        effective_Tc = float(film_results.get("effective_gas_temperature", Tc))
        comp = self.cea_cache.aux.composition(MR, Pc_val, "chamber")
        gas_props_regen = {
            "Pc": Pc_val,
            "Tc": effective_Tc,
            "gamma": gamma,
            "R": R,
            "M": cea_props.get("M"),  # Molecular weight [kg/kmol] for viscosity calculation
            "chamber_area": geometry["area_cross"],
            "A_throat": ensure_chamber_geometry(config).A_throat,
            "chamber_length": geometry["length"],
            "turbulence_intensity": turbulence_intensity_calc,
            "x_H2O": comp["H2O"],
            "x_CO2": comp["CO2"],
        }
 
        coolant_props = {
            "density": float(fuel_fluid.density),
            "viscosity": float(fuel_fluid.viscosity),
            "cp": float(fuel_fluid.specific_heat),
            "thermal_conductivity": float(fuel_fluid.thermal_conductivity),
            "temperature": float(fuel_fluid.temperature),
        }
 
        mdot_coolant = float(film_results.get("mdot_available_for_regen", mdot_F))
 
        if regen_active:
            regen_results = compute_regen_heat_transfer(
                mdot_coolant,
                coolant_props,
                gas_props_regen,
                regen_cfg,
                mdot_total,
            )
            regen_results["mdot_coolant"] = mdot_coolant
            cooling_results["regen"] = regen_results
 
        abl_cfg = config.ablative_cooling
        if abl_cfg is not None and abl_cfg.enabled:
            from engine.pipeline.thermal.gas_side import HotGasState
            cg_now = ensure_chamber_geometry(config)
            contour = self._wall_contour()
            gas = HotGasState(T0=effective_Tc, P0=Pc_val, gamma=gamma,
                              mass_flux_throat=mdot_total / cg_now.A_throat, mu=tr["mu"],
                              cp=tr["cp"], Pr=tr["Pr"], x_H2O=comp["H2O"], x_CO2=comp["CO2"])
            ablative_results = liner_response(gas, contour, abl_cfg, self._liner_end(contour),
                                              mdot_total, with_profile=with_profile)
            ablative_results["incident_heat_flux"] = (ablative_results["heat_flux_from_gas_convective"]
                                                      + ablative_results["heat_flux_from_gas_radiative"])
 
            # Energy the liner takes out of the gas: dT = Q / (mdot cp).
            abl_heat_removed = ablative_results.get("heat_removed", 0.0)
            if abl_heat_removed > 0 and mdot_total > 0:
                cp = gamma * R / max(gamma - 1.0, 1e-6)
                delta_T_abl = abl_heat_removed / max(mdot_total * cp, 1e-6)
                effective_Tc = max(effective_Tc - delta_T_abl, 1.0)
                ablative_results["temperature_reduction"] = float(delta_T_abl)
            else:
                ablative_results["temperature_reduction"] = 0.0
            ablative_results["effective_gas_temperature"] = float(effective_Tc)
            cooling_results["ablative"] = ablative_results
 
        cooling_eff = self._compute_cooling_efficiency(
            cooling_results,
            mdot_total,
            effective_Tc,
            gamma,
            R,
        )
 
        # Store metadata for diagnostics
        metadata = cooling_results.setdefault("metadata", {})
        metadata["gas_turbulence_intensity"] = turbulence_intensity_calc
        metadata["effective_gas_temperature"] = float(effective_Tc)
        metadata["original_gas_temperature"] = float(Tc)
        metadata["gas_viscosity"] = float(mu_g)  # Viscosity used in calculations (calculated from Huzel if available)
        metadata["gas_viscosity_config"] = float(mu_g_config)  # Viscosity from config (for reference)
        metadata["gas_viscosity_calculated"] = float(mu_g_calculated)  # Viscosity from Huzel formula (for reference)
        closure_diag["cooling"] = cooling_results
 
        return cooling_results, cooling_eff, effective_Tc
    
    def _wall_contour(self):
        """The drawn gas-side contour for this geometry (cached per geometry)."""
        from engine.pipeline.thermal.gas_side import contour_for
        cg = ensure_chamber_geometry(self.config)
        key = (cg.A_throat, cg.chamber_diameter, cg.volume, cg.A_exit)
        cached = getattr(self, "_wall_contour_cache", None)
        if cached is None or cached[0] != key:
            cached = (key, contour_for(cg))
            self._wall_contour_cache = cached
        return cached[1]
 
    def _liner_end(self, contour) -> float:
        """Throat-frame x where the ablative liner stops: the graphite insert's upstream edge,
        or the throat when there is no insert."""
        gr = getattr(self.config, "graphite_insert", None)
        if gr is None or not gr.enabled:
            return 0.0
        half = gr.axial_half_length or gr.axial_half_length_ratio * 2.0 * contour.R_t
        return -float(half)
 
    def _infer_injector_diameter(self) -> float:
        """Estimate a characteristic injector diameter for mixing models."""
        injector_cfg = getattr(self.config, "injector", None)
        diameter = None
        injector_type = getattr(injector_cfg, "type", None) if injector_cfg is not None else None
        try:
            if injector_type == "pintle":
                diameter = injector_cfg.geometry.fuel.d_pintle_tip
            elif injector_type == "coaxial":
                diameter = injector_cfg.geometry.core.d_port
            elif injector_type == "impinging":
                diameter = injector_cfg.geometry.oxidizer.d_jet
        except AttributeError:
            diameter = None
        
        if diameter is None or diameter <= 0:
            geometry = self._get_chamber_geometry()
            diameter = np.sqrt(4.0 * geometry["area_cross"] / np.pi)
        return float(max(diameter, 1e-5))
 
    def _get_chamber_geometry(self) -> Dict[str, float]:
        """
        Extract physical chamber geometry from configuration.
        
        Returns a dictionary with:
        - length: Total physical length [m]
        - diameter: Chamber inner diameter [m]
        - area_cross: Cross-sectional area [m²]
        - circumference: Chamber circumference [m]
        - area: Total wetted surface area [m²] (cylindrical + contraction)
        """
        # cg is guaranteed to exist because __init__ calls ensure_chamber_geometry
        cg = ensure_chamber_geometry(self.config)
        regen_cfg = self.config.regen_cooling
 
        # 1. Physical diameter from unified config
        diameter = cg.chamber_diameter
        
        # Fallback to regen if not in unified (though ensure_chamber_geometry should handle it)
        if (diameter is None or diameter <= 0) and regen_cfg is not None and regen_cfg.chamber_inner_diameter is not None:
            diameter = regen_cfg.chamber_inner_diameter
            
        # Final fallback
        if diameter is None or diameter <= 0:
            diameter = 0.08
            
        diameter = max(diameter, 1e-6)
        area_cross = np.pi * (diameter / 2.0)**2
        circumference = np.pi * diameter
        
        # 2. Physical lengths
        length_total = cg.length
        length_cyl = cg.length_cylindrical
        length_cont = cg.length_contraction
        
        # 3. Wetted Surface Area
        # If we have the breakdown (cylindrical + contraction), calculate accurately
        if length_cyl is not None and length_cont is not None:
            # Wetted area = Cylindrical part + Contraction part (frustum of a cone)
            # Area_cyl = pi * D * L_cyl
            area_cyl = circumference * length_cyl
            
            # Area_cont = lateral area of a frustum = pi * (r1 + r2) * slant_height
            r1 = diameter / 2.0
            # Estimate throat radius from A_throat if available
            A_throat = cg.A_throat if cg.A_throat and cg.A_throat > 0 else (area_cross / 3.0)
            r2 = np.sqrt(A_throat / np.pi)
            
            slant_height = np.sqrt((r1 - r2)**2 + length_cont**2)
            area_cont = np.pi * (r1 + r2) * slant_height
            
            area_wetted = area_cyl + area_cont
        else:
            # Fallback to simple cylinder if breakdown not available
            area_wetted = circumference * length_total
 
        return {
            "length": float(length_total),
            "diameter": float(diameter),
            "area_cross": float(area_cross),
            "circumference": float(circumference),
            "area": float(area_wetted),
        }
 
 
    def _rupe_R_opt(self) -> float:
        """Momentum-ratio optimum for the Rupe mixing model.
 
        Honors an explicit ``combustion.efficiency.R_opt`` override; otherwise derives
        it from the impinging-doublet geometry (``sqrt(sin(theta_F)/sin(theta_O))``),
        falling back to 1.0 (balanced) for pintle/coaxial or when angles are absent.
        """
        eff = self.config.combustion.efficiency
        override = getattr(eff, "R_opt", None)
        if override is not None:
            return float(override)
        try:
            from engine.pipeline.combustion_physics import rupe_R_opt_from_angles
            geom = self.config.injector.geometry
            theta_O = getattr(getattr(geom, "oxidizer", None), "impingement_angle", None)
            theta_F = getattr(getattr(geom, "fuel", None), "impingement_angle", None)
            return rupe_R_opt_from_angles(theta_O, theta_F)
        except (AttributeError, TypeError):
            return 1.0
 
    def _fluid_props(self, which: str) -> Optional[Dict[str, float]]:
        """Liquid properties of ``fluids.<which>`` for the c* vaporization march: exactly what the
        config declares (no substituted RP-1 values). Missing keys are left out, and the march
        records each one it has to assume."""
        fl = (getattr(self.config, "fluids", None) or {}).get(which)
        if fl is None:
            return None
        props: Dict[str, float] = {}
        for key in ("density", "boiling_point", "latent_heat", "molecular_weight", "specific_heat",
                    "temperature", "critical_temperature", "viscosity", "surface_tension"):
            v = getattr(fl, key, None)
            if v is not None:
                props[key] = float(v)
        props["Pc_ref"] = float(getattr(fl, "Pc_ref", 2.5e6) or 2.5e6)
        if getattr(fl, "name", None):
            props["name"] = str(fl.name)   # CoolProp liquid cp(T) / density in the heat-up
        return props
 
    def _get_fuel_props(self) -> Optional[Dict[str, float]]:
        """Fuel liquid properties for the c* model (see ``_fluid_props``), plus the fuel-interface
        temperature cap the legacy gasification helper reads."""
        props = self._fluid_props("fuel")
        if props is None:
            return None
        props["T_star_fuel_cap_K"] = float(getattr(self.config.combustion.efficiency, "T_star_fuel_cap_K", 1000.0))
        return props
 
    def _get_ox_props(self) -> Optional[Dict[str, float]]:
        """Oxidizer liquid properties for the c* model (see ``_fluid_props``)."""
        return self._fluid_props("oxidizer")