From 349eea1cad46bcfd1d4612ec6119c55b6684acea Mon Sep 17 00:00:00 2001 From: Carlsaurus Date: Wed, 9 Sep 2026 01:38:00 -0700 Subject: [PATCH 1/5] Stamp the achieved design point; delete a module that never ran chamber_geometry.design_MR / design_pressure / design_thrust are supposed to describe the design the config represents. Nothing ever wrote them. config_schemas builds them with `getattr(chamber, 'design_MR', 2.55)`, so an optimised config carried its template's values forward forever. On a real emitted design they read MR 2.55 / 350 psi / 7000 N for an engine actually solved at O/F 1.65 / 416 psi / 7200 N. Not cosmetic: backend/routers/geometry.py reads design_MR and passes it straight to solve_chamber_geometry_with_cea, so the Chamber Geometry tab drew the contour at the stale mixture ratio -- and 2.55 sits OUTSIDE that config's CEA cache range (MR_range [1.0, 2.5]), i.e. extrapolating past the table edge in the one place the cache has no data. _layer1_stamp_design_point() now writes the solved MR / Pc / F onto the config at the point Layer 1 hands it back, mirrors them onto the legacy `chamber` section for readers that still fall back to it, skips non-finite values so a failed evaluate cannot overwrite a good design point with NaN, and warns when the achieved MR falls outside combustion.cea.MR_range. Verified on a live run: 2.55 / 350 psi / 7000 N in, 1.6461 / 416.1 psi / 7200 N out, matching the solver exactly. Three regression tests cover it -- stamping, the range warning, and the non-finite guard -- because the root problem was that nothing checked these fields at all. Also deletes engine/pipeline/comprehensive_geometry_sizing.py: 495 lines, zero importers, unchanged since the initial commit, and carrying an IndentationError that means it has NEVER been importable. It is recoverable from history if it turns out to be wanted. Audited for the same class of staleness and found two more, NOT fixed here because they belong to other layers and are inert on the configs I can see: - regen_cooling.chamber_inner_diameter / channel_length are synced only by chamber_optimizer.py, not Layer 1, so they still read 0.08491 m / 0.18162 m against an actual 0.127 m / 0.20337 m chamber. Harmless while regen_cooling.enabled is false; wrong channel sizing the moment it is not. - pressure_curves.initial_lox/fuel_pressure_pa (537.3 / 523.7 psi) disagree with lox_tank/fuel_tank.initial_pressure_psi (548.6), because Layer 2 wrote the curves before Layer 1 last moved the tanks. An ordering problem, not a missing write. Separately verified: all 379 config schema fields are referenced somewhere, so the config is not carrying dead knobs; and comprehensive_geometry_sizing was the only module in engine/ or backend/ that fails to import. Verified: pytest 472 passed / 84 skipped with only the 4 pre-existing failures; accelerator parity 16/16. Co-Authored-By: Claude Opus 5 --- .../layers/layer1_static_optimization.py | 69 ++- .../pipeline/comprehensive_geometry_sizing.py | 495 ------------------ .../tests/test_layer1_derived_dofs.py | 58 ++ 3 files changed, 126 insertions(+), 496 deletions(-) delete mode 100644 EngineDesign/engine/pipeline/comprehensive_geometry_sizing.py diff --git a/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py b/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py index 2fc34bd4b..10455ead2 100644 --- a/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py +++ b/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py @@ -1739,6 +1739,70 @@ def _layer1_apply_chamber_geometry_to_config( _DERIVE_AT_SLOPE_MAX = 1.20 +def _layer1_stamp_design_point(config, performance, logger=None) -> None: + """Write the ACHIEVED operating point into ``chamber_geometry.design_*``. + + ``design_MR`` / ``design_pressure`` / ``design_thrust`` are supposed to describe the design + the config represents. Nothing ever wrote them: `config_schemas` builds them with + ``getattr(chamber, 'design_MR', 2.55)``, so an optimised config carried whatever the template + started with. Observed on a real emitted design -- MR 2.55 / 350 psi / 7000 N stamped on an + engine actually solved at O/F ~1.68 / 420 psi / 7200 N. + + That is not cosmetic: ``backend/routers/geometry.py`` reads ``design_MR`` and feeds it + straight into ``solve_chamber_geometry_with_cea``, so the Chamber Geometry tab drew the + contour at the stale mixture ratio -- and 2.55 sits OUTSIDE the shipped CEA cache range + (``MR_range: [1.0, 2.5]``), i.e. extrapolating past the table edge. + + Warns rather than raises when the achieved MR falls outside the cache range: the design is + still real, but anything reading design_MR against that cache is extrapolating. + """ + cg = getattr(config, "chamber_geometry", None) + if cg is None or not isinstance(performance, dict): + return + + def _finite_pos(value): + try: + v = float(value) + except (TypeError, ValueError): + return None + return v if (np.isfinite(v) and v > 0) else None + + mr = _finite_pos(performance.get("MR")) + pc = _finite_pos(performance.get("Pc")) + thrust = _finite_pos(performance.get("F")) + + if mr is not None: + cg.design_MR = mr + if pc is not None: + cg.design_pressure = pc + if thrust is not None: + cg.design_thrust = thrust + + # Keep the legacy mirror in step -- some readers still fall back to config.chamber. + legacy = getattr(config, "chamber", None) + if legacy is not None: + if mr is not None and hasattr(legacy, "design_MR"): + legacy.design_MR = mr + if pc is not None and hasattr(legacy, "design_pressure"): + legacy.design_pressure = pc + if thrust is not None and hasattr(legacy, "design_thrust"): + legacy.design_thrust = thrust + + if mr is None or logger is None: + return + try: + mr_range = config.combustion.cea.MR_range + lo, hi = float(mr_range[0]), float(mr_range[1]) + except (AttributeError, TypeError, ValueError, IndexError): + return + if not (lo <= mr <= hi): + logger.warning( + "design_MR %.4f is outside the CEA cache MR_range [%.2f, %.2f]; anything reading " + "it against that cache is extrapolating past the table edge.", mr, lo, hi + ) + + + def _layer1_eps_for_exit_pressure(Pc_Pa: float, gamma: float, Pe_Pa: float): """Expansion ratio that puts the exit plane exactly at ``Pe_Pa``. @@ -8563,7 +8627,10 @@ def _as_finite_float_or_nan(v: Any) -> float: layer1_logger.handlers.clear() update_progress("Layer 1: Complete", 1.0, "Layer 1 optimization complete!") - + + # Stamp the ACHIEVED operating point onto the config we are about to hand back, so + # chamber_geometry.design_* describes this engine rather than whatever template it came from. + _layer1_stamp_design_point(optimized_config, final_performance, layer1_logger) return optimized_config, results diff --git a/EngineDesign/engine/pipeline/comprehensive_geometry_sizing.py b/EngineDesign/engine/pipeline/comprehensive_geometry_sizing.py deleted file mode 100644 index 7caf1f396..000000000 --- a/EngineDesign/engine/pipeline/comprehensive_geometry_sizing.py +++ /dev/null @@ -1,495 +0,0 @@ -"""Comprehensive geometry sizing and visualization for chamber, throat, and ablative. - -This module provides: -1. Optimal sizing of ablative and throat together -2. Combined visualization (plot + DXF) showing all three components -3. Robust solver with error handling and validation -""" - -from __future__ import annotations - -from typing import Dict, Any, Optional, Tuple, List -import numpy as np -import matplotlib.pyplot as plt -from matplotlib.patches import Rectangle, Circle, FancyBboxPatch -from matplotlib.collections import PatchCollection -import io - -from .config_schemas import ( - PintleEngineConfig, - AblativeCoolingConfig, - GraphiteInsertConfig, - StainlessSteelCaseConfig, - ensure_chamber_geometry, -) -from engine.pipeline.thermal.ablative_sizing import size_ablative_system -from engine.pipeline.thermal.graphite_geometry import size_graphite_insert as size_graphite_geom -from engine.core.chamber_profiles import calculate_complete_chamber_geometry - - -def size_complete_geometry( - config: PintleEngineConfig, - Pc: float, - MR: float, - Tc: float, - gamma: float, - R: float, - burn_time: float, - chamber_heat_flux: float, - throat_heat_flux_multiplier: float = 1.5, -) -> Dict[str, Any]: - """ - Size complete geometry: chamber (ablative), throat (graphite), and all components together. - - This function: - 1. Sizes ablative liner for chamber - 2. Sizes graphite insert for throat (with zero recession) - 3. Validates all sizing meets requirements - 4. Returns optimal geometry configuration - - Parameters: - ----------- - config : PintleEngineConfig - Engine configuration - Pc : float - Chamber pressure [Pa] - MR : float - Mixture ratio - Tc : float - Chamber temperature [K] - gamma : float - Specific heat ratio - R : float - Gas constant [J/(kg·K)] - burn_time : float - Burn time [s] - chamber_heat_flux : float - Chamber heat flux [W/m²] - throat_heat_flux_multiplier : float - Multiplier for throat heat flux vs chamber (default 1.5) - - Returns: - -------- - sizing_results : dict - Complete sizing results including: - - ablative_sizing: Ablative thickness and properties - - graphite_sizing: Graphite insert sizing - - geometry: Complete geometry profile - - validation: Validation results - - optimal: Optimal configuration selected - """ - results = { - "ablative_sizing": None, - "graphite_sizing": None, - "geometry": None, - "validation": {}, - "optimal": {}, - } - - # 1. Size ablative system for chamber - if config.ablative_cooling and config.ablative_cooling.enabled: - ablative_sizing = size_ablative_system( - heat_flux=chamber_heat_flux, - burn_time=burn_time, - ablative_config=config.ablative_cooling, - backface_temp_limit=500.0, # K - Max for stainless steel - T_hot_gas=Tc, - h_hot_gas=5000.0, # W/(m²·K) - Typical for rocket chambers - q_rad_hot=0.0, # Negligible for LOX/RP-1 - ) - results["ablative_sizing"] = ablative_sizing - else: - results["ablative_sizing"] = {"required_thickness": 0.0, "meets_requirements": True} - - # 2. Size graphite insert for throat (with ZERO recession - that's its purpose) - if config.graphite_insert and config.graphite_insert.enabled: - throat_heat_flux = chamber_heat_flux * throat_heat_flux_multiplier - - # Get throat conditions - # Surface temperature estimate (throat is hottest) - surface_temp_throat = Tc * 0.85 # Conservative estimate - - # CRITICAL: Graphite recession should be ZERO for sizing - # The whole point is that graphite doesn't ablate - it keeps throat constant - # Use a small value only for sizing calculations (material allowance), not runtime - recession_rate_for_sizing = 1e-8 # Negligible - graphite doesn't ablate - - # Get throat diameter from config or calculate - cg = ensure_chamber_geometry(config) - if cg.A_throat: - A_throat = cg.A_throat - D_throat = np.sqrt(4.0 * A_throat / np.pi) - else: - # Estimate from typical expansion ratio - D_throat = 0.020 # 20 mm default - - # Use graphite_geometry.size_graphite_insert (returns GraphiteInsertSizing dataclass) - graphite_sizing_obj = size_graphite_geom( - peak_heat_flux=throat_heat_flux, - surface_temperature=surface_temp_throat, - recession_rate=recession_rate_for_sizing, # Negligible - graphite doesn't ablate - burn_time=burn_time, - thermal_conductivity=config.graphite_insert.thermal_conductivity, - backface_temperature_max=500.0, # K - Max for stainless steel - throat_diameter=D_throat, - density=config.graphite_insert.material_density, - specific_heat=config.graphite_insert.specific_heat, - mechanical_thickness=0.001, # 1 mm - safety_factor=0.3, # 30% - transient=True, - ) - # Convert to dict for compatibility - graphite_sizing = graphite_sizing_obj.to_dict() - graphite_sizing["meets_requirements"] = not graphite_sizing_obj.throat_area_change_excessive - results["graphite_sizing"] = graphite_sizing - else: - results["graphite_sizing"] = {"initial_thickness": 0.0, "meets_requirements": True} - - # 3. Calculate complete geometry - # Get geometry from chamber_geometry - cg = ensure_chamber_geometry(config) - V_chamber = cg.volume - A_throat = cg.A_throat - L_chamber = cg.length if cg.length else (cg.volume / cg.A_throat if cg.A_throat and cg.A_throat > 0 else 0.18) - - # Calculate diameters - if L_chamber > 0: - D_chamber_initial = np.sqrt(4.0 * V_chamber / (np.pi * L_chamber)) - else: - D_chamber_initial = np.sqrt(4.0 * V_chamber / np.pi) # Assume cylindrical - D_throat_initial = np.sqrt(4.0 * A_throat / np.pi) if A_throat > 0 else 0.020 - else: - # Fallback estimates - V_chamber = 0.001 # 1 L - A_throat = np.pi * (0.010) ** 2 # 20 mm diameter - L_chamber = 0.1 # 10 cm - D_chamber_initial = 0.05 # 50 mm - D_throat_initial = 0.020 # 20 mm - - geometry = calculate_complete_chamber_geometry( - V_chamber=V_chamber, - A_throat=A_throat, - L_chamber=L_chamber, - D_chamber_initial=D_chamber_initial, - D_throat_initial=D_throat_initial, - ablative_config=config.ablative_cooling if config.ablative_cooling else None, - graphite_config=config.graphite_insert if config.graphite_insert else None, - stainless_config=config.stainless_steel_case if hasattr(config, "stainless_steel_case") else None, - recession_chamber=0.0, # Initial state - recession_graphite=0.0, # Graphite doesn't recede - n_points=100, - ) - results["geometry"] = geometry - - # 4. Validate sizing - validation = { - "ablative_meets_requirements": results["ablative_sizing"].get("meets_requirements", True), - "graphite_meets_requirements": results["graphite_sizing"].get("meets_requirements", True), - "all_valid": True, - "warnings": [], - } - - if config.ablative_cooling and config.ablative_cooling.enabled: - if not validation["ablative_meets_requirements"]: - validation["warnings"].append("Ablative backface temperature exceeds limit") - validation["all_valid"] = False - - if config.graphite_insert and config.graphite_insert.enabled: - if not validation["graphite_meets_requirements"]: - validation["warnings"].append("Graphite backface temperature exceeds limit") - validation["all_valid"] = False - - # Check graphite thickness is reasonable - graphite_thickness = results["graphite_sizing"].get("initial_thickness", 0.0) - if graphite_thickness < 0.001: # Less than 1 mm - validation["warnings"].append("Graphite thickness is very small - may not provide adequate protection") - if graphite_thickness > 0.010: # More than 10 mm - validation["warnings"].append("Graphite thickness is very large - consider optimization") - - results["validation"] = validation - - # 5. Select optimal configuration - optimal = { - "ablative_thickness": results["ablative_sizing"].get("required_thickness", 0.0), - "graphite_thickness": results["graphite_sizing"].get("initial_thickness", 0.0), - "throat_diameter": D_throat_initial, - "chamber_diameter": D_chamber_initial, - "chamber_length": L_chamber, - "total_mass": 0.0, # Could calculate if needed - "meets_all_requirements": validation["all_valid"], - } - - # Calculate total mass (rough estimate) - if config.ablative_cooling and config.ablative_cooling.enabled: - ablative_density = config.ablative_cooling.material_density - ablative_volume = np.pi * L_chamber * ( - (D_chamber_initial / 2.0 + optimal["ablative_thickness"]) ** 2 - - (D_chamber_initial / 2.0) ** 2 - ) - optimal["ablative_mass"] = ablative_density * ablative_volume - else: - optimal["ablative_mass"] = 0.0 - - if config.graphite_insert and config.graphite_insert.enabled: - graphite_density = config.graphite_insert.material_density - # Approximate graphite as cylinder around throat - graphite_length = optimal["graphite_thickness"] * 2.0 # Rough estimate - graphite_volume = np.pi * graphite_length * ( - (D_throat_initial / 2.0 + optimal["graphite_thickness"]) ** 2 - - (D_throat_initial / 2.0) ** 2 - ) - optimal["graphite_mass"] = graphite_density * graphite_volume - else: - optimal["graphite_mass"] = 0.0 - - optimal["total_mass"] = optimal["ablative_mass"] + optimal["graphite_mass"] - - results["optimal"] = optimal - - return results - - -def plot_complete_geometry( - sizing_results: Dict[str, Any], - config: PintleEngineConfig, - save_path: Optional[str] = None, - show_graphite: bool = True, - show_ablative: bool = True, - show_stainless: bool = True, - use_plotly: bool = True, -) -> Tuple[Any, bytes]: - """ - Create comprehensive plot showing chamber, throat, ablative, and graphite all together. - - Parameters: - ----------- - sizing_results : dict - Results from size_complete_geometry() - config : PintleEngineConfig - Engine configuration - save_path : str, optional - Path to save figure (if None, not saved) - show_graphite : bool - Show graphite insert (default True) - show_ablative : bool - Show ablative liner (default True) - show_stainless : bool - Show stainless steel case (default True) - use_plotly : bool - Use Plotly for interactive plots (default True), otherwise matplotlib - - Returns: - -------- - fig : plotly.Figure or matplotlib.Figure - Figure object - dxf_bytes : bytes - DXF file bytes (placeholder - would need dxf library) - """ - geometry = sizing_results["geometry"] - positions = np.array(geometry["positions"]) - - if use_plotly: - import plotly.graph_objects as go - - fig = go.Figure() - - # Chamber gas boundary (inner surface) - orange - D_gas = np.array(geometry.get("D_gas_chamber", geometry.get("D_chamber_current", np.zeros_like(positions)))) - if isinstance(D_gas, (int, float)) or len(D_gas) == 1: - D_gas = np.full_like(positions, float(D_gas) if isinstance(D_gas, (int, float)) else D_gas[0]) - D_gas_radius = D_gas / 2.0 - - fig.add_trace(go.Scatter( - x=positions, - y=D_gas_radius, - mode='lines', - name='Gas Boundary (Chamber)', - line=dict(color='orange', width=3), - fill='tozeroy', - fillcolor='rgba(255, 165, 0, 0.1)', - )) - fig.add_trace(go.Scatter( - x=positions, - y=-D_gas_radius, - mode='lines', - name='Gas Boundary (Lower)', - line=dict(color='orange', width=3), - fill='tozeroy', - fillcolor='rgba(255, 165, 0, 0.1)', - showlegend=False, - )) - - # Ablative layer - brown dashed - if show_ablative and geometry.get("ablative_thickness", [0.0])[0] > 0: - D_ablative = np.array(geometry.get("D_ablative_outer", D_gas)) - if isinstance(D_ablative, (int, float)) or len(D_ablative) == 1: - D_ablative = np.full_like(positions, float(D_ablative) if isinstance(D_ablative, (int, float)) else D_ablative[0]) - D_ablative_radius = D_ablative / 2.0 - - fig.add_trace(go.Scatter( - x=positions, - y=D_ablative_radius, - mode='lines', - name='Phenolic Ablator (Outer)', - line=dict(color='brown', width=2, dash='dash'), - fill='tonexty', - fillcolor='rgba(139, 69, 19, 0.3)', - )) - fig.add_trace(go.Scatter( - x=positions, - y=-D_ablative_radius, - mode='lines', - name='Phenolic Ablator (Lower)', - line=dict(color='brown', width=2, dash='dash'), - fill='tonexty', - fillcolor='rgba(139, 69, 19, 0.3)', - showlegend=False, - )) - - # Stainless steel case - gray dotted - if show_stainless and geometry.get("stainless_thickness", 0.0) > 0: - D_stainless = np.array(geometry.get("D_stainless_outer", D_gas)) - if isinstance(D_stainless, (int, float)) or len(D_stainless) == 1: - D_stainless = np.full_like(positions, float(D_stainless) if isinstance(D_stainless, (int, float)) else D_stainless[0]) - D_stainless_radius = D_stainless / 2.0 - - fig.add_trace(go.Scatter( - x=positions, - y=D_stainless_radius, - mode='lines', - name='Stainless Steel Case', - line=dict(color='gray', width=2, dash='dot'), - fill='tonexty', - fillcolor='rgba(128, 128, 128, 0.2)', - )) - fig.add_trace(go.Scatter( - x=positions, - y=-D_stainless_radius, - mode='lines', - name='Stainless Steel (Lower)', - line=dict(color='gray', width=2, dash='dot'), - fill='tonexty', - fillcolor='rgba(128, 128, 128, 0.2)', - showlegend=False, - )) - - # Throat region with graphite - ONLY at throat, not entire chamber - if show_graphite and config.graphite_insert and config.graphite_insert.enabled: - D_throat = geometry.get("D_throat_current", 0.020) - D_graphite_outer = geometry.get("D_graphite_outer", D_throat) - throat_pos = positions[-1] if len(positions) > 0 else 0.0 - - # Graphite axial length (typically 0.75 * D_throat on each side) - D_throat_diameter = D_throat - graphite_axial_half_length = getattr(config.graphite_insert, 'axial_half_length', 0.75 * D_throat_diameter) - if graphite_axial_half_length <= 0: - graphite_axial_half_length = 0.75 * D_throat_diameter - - # Graphite region (ONLY around throat) - graphite_start = max(throat_pos - graphite_axial_half_length, positions[0]) - graphite_end = min(throat_pos + graphite_axial_half_length, positions[-1]) - graphite_positions = np.linspace(graphite_start, graphite_end, 30) - D_graphite_radius = D_graphite_outer / 2.0 - D_throat_radius = D_throat / 2.0 - - # Graphite outer boundary (black, ONLY in throat region) - fig.add_trace(go.Scatter( - x=graphite_positions, - y=[D_graphite_radius] * len(graphite_positions), - mode='lines', - name='Graphite Insert', - line=dict(color='black', width=3), - )) - fig.add_trace(go.Scatter( - x=graphite_positions, - y=[-D_graphite_radius] * len(graphite_positions), - mode='lines', - name='Graphite Insert (Lower)', - line=dict(color='black', width=3), - showlegend=False, - )) - - # Throat (red marker at minimum diameter) - fig.add_trace(go.Scatter( - x=[throat_pos], - y=[D_throat_radius], - mode='markers', - marker=dict(size=12, color='red', symbol='circle', line=dict(width=2, color='darkred')), - name='Throat', - showlegend=True, - )) - fig.add_trace(go.Scatter( - x=[throat_pos], - y=[-D_throat_radius], - mode='markers', - marker=dict(size=12, color='red', symbol='circle', line=dict(width=2, color='darkred')), - showlegend=False, - )) - - # Centerline - fig.add_hline(y=0, line_dash="dash", line_color="gray", opacity=0.5) - - fig.update_layout( - title="Complete Chamber Geometry: Chamber, Throat, Ablative, and Graphite", - xaxis_title="Axial Position [m]", - yaxis_title="Radius [m]", - height=600, - showlegend=True, - yaxis=dict(scaleanchor="x", scaleratio=1), # Equal aspect ratio - ) - - if save_path: - fig.write_image(save_path) - - return fig, b"" # DXF placeholder - - else: - # Matplotlib version (fallback) - fig, ax = plt.subplots(figsize=(14, 8)) - ax.set_aspect('equal') - - # Similar implementation with matplotlib - # (Keep existing matplotlib code as fallback) - - return fig, b"" - - -def select_optimal_geometry( - config: PintleEngineConfig, - design_requirements: Dict[str, Any], -) -> Dict[str, Any]: - """ - Select optimal geometry configuration from multiple sizing options. - - This function evaluates multiple geometry configurations and selects the best one - based on requirements (mass, performance, manufacturability, etc.). - - Parameters: - ----------- - config : PintleEngineConfig - Base engine configuration - design_requirements : dict - Design requirements including: - - target_thrust: Target thrust [N] - - burn_time: Burn time [s] - - max_mass: Maximum total mass [kg] - - min_performance: Minimum Isp [s] - - constraints: Additional constraints - - Returns: - -------- - optimal_config : dict - Optimal configuration selected - """ - # This is a placeholder - would implement full optimization here - # For now, return the input config with validation - - optimal = { - "config": config, - "meets_requirements": True, - "score": 1.0, - "reasoning": "Configuration meets all requirements", - } - - return optimal - diff --git a/EngineDesign/tests/test_layer1_derived_dofs.py b/EngineDesign/tests/test_layer1_derived_dofs.py index ac6cbdb74..0a3ac05d1 100644 --- a/EngineDesign/tests/test_layer1_derived_dofs.py +++ b/EngineDesign/tests/test_layer1_derived_dofs.py @@ -804,3 +804,61 @@ def test_momentum_band_is_symmetric_now_that_tilt_owns_the_guard(): assert _mom(0.97, band_width=0.05) == 0.0 assert _mom(1.08, band_width=0.05) > 0.0 assert _mom(0.92, band_width=0.05) > 0.0 + + +def test_design_point_is_stamped_from_the_solved_engine(): + """chamber_geometry.design_* must describe THIS engine, not the template it came from. + + Nothing ever wrote these: config_schemas builds them with + `getattr(chamber, 'design_MR', 2.55)`, so an optimised config carried the template's values + forward forever. A real emitted design was stamped MR 2.55 / 350 psi / 7000 N while actually + solving at O/F ~1.65 / 416 psi / 7200 N -- and backend/routers/geometry.py feeds design_MR + straight into solve_chamber_geometry_with_cea, so the geometry tab drew the contour at the + stale ratio. + """ + import logging + from engine.pipeline.io import load_config + from engine.optimizer.layers.layer1_static_optimization import _layer1_stamp_design_point + + cfg = load_config("configs/canonical/impinging.yaml") + cfg.chamber_geometry.design_MR = 2.55 # the stale template values + cfg.chamber_geometry.design_pressure = 2.413166e6 + cfg.chamber_geometry.design_thrust = 7000.0 + + solved = {"MR": 1.6461, "Pc": 2.8690e6, "F": 7200.0} + _layer1_stamp_design_point(cfg, solved, None) + + assert cfg.chamber_geometry.design_MR == pytest.approx(solved["MR"]) + assert cfg.chamber_geometry.design_pressure == pytest.approx(solved["Pc"]) + assert cfg.chamber_geometry.design_thrust == pytest.approx(solved["F"]) + + +def test_design_point_stamp_warns_outside_the_cea_cache_range(): + """An MR outside combustion.cea.MR_range means anything reading it extrapolates.""" + import logging + from engine.pipeline.io import load_config + from engine.optimizer.layers.layer1_static_optimization import _layer1_stamp_design_point + + cfg = load_config("configs/canonical/impinging.yaml") + lo, hi = [float(v) for v in cfg.combustion.cea.MR_range] + + seen = [] + logger = logging.getLogger("stamp_range_test") + logger.warning = lambda msg, *a, **k: seen.append(msg % a if a else msg) + + _layer1_stamp_design_point(cfg, {"MR": (lo + hi) / 2.0, "Pc": 2.8e6, "F": 7200.0}, logger) + assert not seen, "in-range MR must not warn" + + _layer1_stamp_design_point(cfg, {"MR": hi + 1.0, "Pc": 2.8e6, "F": 7200.0}, logger) + assert seen and "outside the CEA cache" in seen[0] + + +def test_design_point_stamp_ignores_non_finite_performance(): + """A failed evaluate must not overwrite a good design point with NaN.""" + from engine.pipeline.io import load_config + from engine.optimizer.layers.layer1_static_optimization import _layer1_stamp_design_point + + cfg = load_config("configs/canonical/impinging.yaml") + cfg.chamber_geometry.design_MR = 1.65 + _layer1_stamp_design_point(cfg, {"MR": float("nan"), "Pc": 0.0, "F": None}, None) + assert cfg.chamber_geometry.design_MR == pytest.approx(1.65) From 68d1a5f30dc3510e7a81adb6d963698c1487d7a9 Mon Sep 17 00:00:00 2001 From: Carlsaurus Date: Wed, 9 Sep 2026 01:46:45 -0700 Subject: [PATCH 2/5] Flag Layer-2 pressure curves left stale by a Layer-1 re-run The initial tank pressure exists twice, owned by different layers and never reconciled: Layer 1 writes lox_tank/fuel_tank.initial_pressure_psi (layer1_static_optimization.py:2435), Layer 2 writes pressure_curves.initial_lox/fuel_pressure_pa (layer2_pressure.py:312). Re-running Layer 1 moves the tanks and silently leaves the curves describing the previous design. Measured on a real emitted config: LOX 548.6 psi vs a 537.3 psi curve start (11.3 psi), fuel 548.6 vs 523.7 (24.9 psi). This matters more than it looks. Tank pressure is the UPSTREAM BOUNDARY CONDITION for the feed-system twin, and per docs/adr/0001 EngineDesign's optimizer imports lib/feedtwin directly for Layer X rather than calling a service. Two disagreeing values for one boundary condition is precisely what silently poisons a twin, so Layer 1 now says so at the end of every run. Detection only. Which layer should own the value is a design call and guessing it here would be worse than the warning. Co-Authored-By: Claude Opus 5 --- .../layers/layer1_static_optimization.py | 46 +++++++++++++++++++ .../tests/test_layer1_derived_dofs.py | 36 +++++++++++++++ 2 files changed, 82 insertions(+) diff --git a/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py b/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py index 10455ead2..33f3511ac 100644 --- a/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py +++ b/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py @@ -1803,6 +1803,51 @@ def _finite_pos(value): +def _layer1_warn_stale_pressure_curves(config, logger=None, tol_psi: float = 1.0) -> None: + """Flag Layer-2 pressure curves that no longer match the tank pressures Layer 1 just set. + + The initial tank pressure exists twice, owned by different layers and never reconciled: + Layer 1 writes ``lox_tank/fuel_tank.initial_pressure_psi``; Layer 2 writes + ``pressure_curves.initial_lox/fuel_pressure_pa``. Re-running Layer 1 moves the tanks and + silently leaves the curves describing the previous design -- observed 11.3 psi out on the + LOX side and 24.9 psi on the fuel side of a real emitted config. + + This matters beyond EngineDesign. Tank pressure is the UPSTREAM BOUNDARY CONDITION for the + feed-system twin (docs/adr/0001), which EngineDesign's optimizer will import directly for + Layer X. Two disagreeing values for one boundary condition is exactly the kind of thing that + silently poisons a twin, so say so loudly rather than letting it cross the boundary. + + Detection only -- which layer should win is a design call, not something to guess here. + """ + if logger is None: + return + curves = getattr(config, "pressure_curves", None) + if curves is None: + return + PSI = 6894.76 + for tank_attr, curve_attr, label in ( + ("lox_tank", "initial_lox_pressure_pa", "LOX"), + ("fuel_tank", "initial_fuel_pressure_pa", "fuel"), + ): + tank = getattr(config, tank_attr, None) + if tank is None: + continue + try: + tank_psi = float(getattr(tank, "initial_pressure_psi", float("nan"))) + curve_psi = float(getattr(curves, curve_attr, float("nan"))) / PSI + except (TypeError, ValueError): + continue + if not (np.isfinite(tank_psi) and np.isfinite(curve_psi)): + continue + if abs(tank_psi - curve_psi) > tol_psi: + logger.warning( + "%s tank pressure disagrees with the Layer-2 pressure curve: tank %.1f psi vs " + "curve start %.1f psi (%.1f psi apart). The curves predate this Layer 1 run; " + "re-run Layer 2 before trusting them or anything downstream of them.", + label, tank_psi, curve_psi, abs(tank_psi - curve_psi), + ) + + def _layer1_eps_for_exit_pressure(Pc_Pa: float, gamma: float, Pe_Pa: float): """Expansion ratio that puts the exit plane exactly at ``Pe_Pa``. @@ -8631,6 +8676,7 @@ def _as_finite_float_or_nan(v: Any) -> float: # Stamp the ACHIEVED operating point onto the config we are about to hand back, so # chamber_geometry.design_* describes this engine rather than whatever template it came from. _layer1_stamp_design_point(optimized_config, final_performance, layer1_logger) + _layer1_warn_stale_pressure_curves(optimized_config, layer1_logger) return optimized_config, results diff --git a/EngineDesign/tests/test_layer1_derived_dofs.py b/EngineDesign/tests/test_layer1_derived_dofs.py index 0a3ac05d1..e0a68d0d9 100644 --- a/EngineDesign/tests/test_layer1_derived_dofs.py +++ b/EngineDesign/tests/test_layer1_derived_dofs.py @@ -862,3 +862,39 @@ def test_design_point_stamp_ignores_non_finite_performance(): cfg.chamber_geometry.design_MR = 1.65 _layer1_stamp_design_point(cfg, {"MR": float("nan"), "Pc": 0.0, "F": None}, None) assert cfg.chamber_geometry.design_MR == pytest.approx(1.65) + + +def test_stale_pressure_curves_are_flagged(): + """The initial tank pressure exists twice, owned by different layers, never reconciled. + + Layer 1 writes lox_tank/fuel_tank.initial_pressure_psi; Layer 2 writes + pressure_curves.initial_lox/fuel_pressure_pa. Re-running Layer 1 silently leaves the curves + describing the previous design (observed 11.3 psi out on LOX, 24.9 psi on fuel). Tank + pressure is the upstream boundary condition for the feed-system twin (docs/adr/0001), so a + disagreement must not cross that boundary unannounced. + """ + import logging + from engine.pipeline.io import load_config + from engine.optimizer.layers.layer1_static_optimization import ( + _layer1_warn_stale_pressure_curves, + ) + + cfg = load_config("configs/canonical/impinging.yaml") + if getattr(cfg, "pressure_curves", None) is None or getattr(cfg, "lox_tank", None) is None: + pytest.skip("config has no pressure_curves / lox_tank to compare") + + seen = [] + logger = logging.getLogger("stale_curves_test") + logger.warning = lambda msg, *a, **k: seen.append(msg % a if a else msg) + + PSI = 6894.76 + cfg.lox_tank.initial_pressure_psi = 548.6 + cfg.pressure_curves.initial_lox_pressure_pa = 548.6 * PSI # agrees + cfg.fuel_tank.initial_pressure_psi = 548.6 + cfg.pressure_curves.initial_fuel_pressure_pa = 548.6 * PSI # agrees + _layer1_warn_stale_pressure_curves(cfg, logger) + assert not seen, "matching pressures must not warn" + + cfg.pressure_curves.initial_lox_pressure_pa = 537.3 * PSI # 11.3 psi stale + _layer1_warn_stale_pressure_curves(cfg, logger) + assert seen and "LOX tank pressure disagrees" in seen[0] From c55612622f085fe4756aefa23c61426734b81246 Mon Sep 17 00:00:00 2001 From: Carlsaurus Date: Wed, 9 Sep 2026 01:57:53 -0700 Subject: [PATCH 3/5] Make the pressurisation path configurable and its state per-config Two problems docs/adr/0001 names by name, both in the COPV -> regulator -> tank chain that lib/feedtwin absorbs. Fixed surgically: defaults are exactly the old literals, so nothing in EngineDesign changes behaviour unless someone overrides them. 1. Orifice geometry was hardcoded in the function body. dynamics.step() defined gamma_gas 1.4, Cd_regulator 0.7, Cd_valve 0.65, A_regulator 2e-5, A_valve_F/O 5e-5 as locals. These are real hardware numbers -- a different regulator or solenoid could not be modelled without editing code. They are now DynamicsParams fields, surfaced on ControllerConfig, and read through from_config with getattr so an older config still loads. 2. The polytropic reference state was PROCESS-GLOBAL. step() stored its reference temperatures, masses and volumes as attributes on the FUNCTION object, guarded by `hasattr(step, '_temp_initialized')`. It initialised on the first call ever made in a process and never reset, so the first trajectory stepped set the reference for every later one -- across different configs. An optimizer evaluating thousands of candidates, which is exactly the Layer X access pattern the ADR describes, would hand candidates 2..N the reference state of candidate 1. All 36 references now live on the `params` instance, which scopes them per-config. Verified: two DynamicsParams objects no longer share reference state, and nothing is left on the step function object. Four regression tests cover both fixes, including one that fails if process-global state ever returns. Also checked the other two feed fragments the ADR names, feed_loss.py and stability/chug.py: neither carries hardcoded physical constants, so they need nothing before Phase 04 collapses them into feedtwin. Verified: pytest 477 passed / 84 skipped with only the 4 pre-existing failures; robust_ddp suite 12 passed; accelerator parity 16/16. Co-Authored-By: Claude Opus 5 --- .../engine/control/robust_ddp/data_models.py | 9 ++ .../engine/control/robust_ddp/dynamics.py | 120 +++++++++++------- .../robust_ddp/test_robust_ddp_dynamics.py | 59 +++++++++ 3 files changed, 141 insertions(+), 47 deletions(-) diff --git a/EngineDesign/engine/control/robust_ddp/data_models.py b/EngineDesign/engine/control/robust_ddp/data_models.py index b10df9b36..54eb122f7 100644 --- a/EngineDesign/engine/control/robust_ddp/data_models.py +++ b/EngineDesign/engine/control/robust_ddp/data_models.py @@ -267,6 +267,15 @@ class ControllerConfig: # Regulator model reg_setpoint: Optional[float] = None # Regulator setpoint [Pa] (None = derived from COPV) reg_ratio: float = 0.8 # P_reg / P_copv ratio if setpoint not specified + # Pressurisation-path hardware. Previously literals inside dynamics.step(); surfaced here so + # a vehicle with a different regulator or solenoid can be modelled without editing code. + # See docs/adr/0001 -- this chain is one of the feed models lib/feedtwin absorbs. + gamma_gas: float = 1.4 # Pressurant specific heat ratio (N2) + Cd_regulator: float = 0.7 # Discharge coefficient, regulator orifice + Cd_valve: float = 0.65 # Discharge coefficient, solenoid valve orifice + A_regulator: float = 2e-5 # Regulator orifice area [m^2] + A_valve_F: float = 5e-5 # Fuel solenoid valve flow area [m^2] + A_valve_O: float = 5e-5 # Oxidiser solenoid valve flow area [m^2] # Ullage pressurization flow coefficients [1/s] alpha_F: float = 10.0 # Fuel pressurization flow coefficient diff --git a/EngineDesign/engine/control/robust_ddp/dynamics.py b/EngineDesign/engine/control/robust_ddp/dynamics.py index 2021ea67f..644766c23 100644 --- a/EngineDesign/engine/control/robust_ddp/dynamics.py +++ b/EngineDesign/engine/control/robust_ddp/dynamics.py @@ -63,6 +63,19 @@ class DynamicsParams: T_gas_copv_initial: float = 293.0 # Initial COPV gas temperature [K] T_gas_F_initial: float = 293.0 # Initial fuel tank gas temperature [K] T_gas_O_initial: float = 250.0 # Initial oxidizer tank gas temperature [K] (LOX tank is colder) + + # --- Pressurisation-path hardware ------------------------------------------------------ + # These were literals inside step(): a COPV -> regulator -> tank chain with the orifice + # geometry baked into the function body, which is one of the four fragmented feed models + # docs/adr/0001 calls out. They are real hardware numbers a user must be able to set once + # this path is calibrated or replaced by lib/feedtwin. Defaults are exactly the values that + # were hardcoded, so behaviour is unchanged unless someone overrides them. + gamma_gas: float = 1.4 # Specific heat ratio of the pressurant (N2) + Cd_regulator: float = 0.7 # Discharge coefficient, regulator orifice + Cd_valve: float = 0.65 # Discharge coefficient, solenoid valve orifice + A_regulator: float = 2e-5 # Regulator orifice area [m^2] (~20 mm^2) + A_valve_F: float = 5e-5 # Fuel solenoid valve flow area [m^2] (~50 mm^2) + A_valve_O: float = 5e-5 # Oxidiser solenoid valve flow area [m^2] (~50 mm^2) @classmethod def from_config(cls, config: ControllerConfig) -> DynamicsParams: @@ -90,6 +103,13 @@ def from_config(cls, config: ControllerConfig) -> DynamicsParams: T_gas_copv_initial=getattr(config, 'T_gas_copv_initial', 293.0), # COPV initial temp T_gas_F_initial=getattr(config, 'T_gas_F_initial', 293.0), # Fuel tank initial temp T_gas_O_initial=getattr(config, 'T_gas_O_initial', 250.0), # LOX tank initial temp (colder) + # Pressurisation-path hardware; getattr so an older ControllerConfig still loads. + gamma_gas=getattr(config, 'gamma_gas', 1.4), + Cd_regulator=getattr(config, 'Cd_regulator', 0.7), + Cd_valve=getattr(config, 'Cd_valve', 0.65), + A_regulator=getattr(config, 'A_regulator', 2e-5), + A_valve_F=getattr(config, 'A_valve_F', 5e-5), + A_valve_O=getattr(config, 'A_valve_O', 5e-5), ) @@ -160,17 +180,22 @@ def step( # n = 1.4: adiabatic (no heat transfer, γ for diatomic gas) # n = 1.2: typical for blowdown (some heat transfer with tank walls) - # Initialize temperatures on first call (store in function attributes) - if not hasattr(step, '_temp_initialized'): - step._T_copv_0 = getattr(params, 'T_gas_copv_initial', params.T_gas) - step._T_F_0 = getattr(params, 'T_gas_F_initial', params.T_gas) - step._T_O_0 = getattr(params, 'T_gas_O_initial', 250.0) # LOX tank colder - step._m_copv_0 = m_gas_copv - step._m_F_0 = m_gas_F - step._m_O_0 = m_gas_O - step._V_F_0 = V_u_F - step._V_O_0 = V_u_O - step._temp_initialized = True + # Reference state for the polytropic relation, captured on this params object's first + # step. It used to live on the FUNCTION object as private attributes, which made it + # process-global: the first trajectory ever stepped set the reference for every later + # one, across configs. An optimizer evaluating thousands of candidates -- exactly the + # Layer X access pattern docs/adr/0001 describes -- would hand candidates 2..N the + # reference state of candidate 1. Scoping it to `params` makes it per-config. + if not getattr(params, '_temp_initialized', False): + params._T_copv_0 = getattr(params, 'T_gas_copv_initial', params.T_gas) + params._T_F_0 = getattr(params, 'T_gas_F_initial', params.T_gas) + params._T_O_0 = getattr(params, 'T_gas_O_initial', 250.0) # LOX tank colder + params._m_copv_0 = m_gas_copv + params._m_F_0 = m_gas_F + params._m_O_0 = m_gas_O + params._V_F_0 = V_u_F + params._V_O_0 = V_u_O + params._temp_initialized = True # Compute current gas temperatures using polytropic relation # T = T0 * (rho/rho0)^(n-1) = T0 * (m/V) / (m0/V0))^(n-1) @@ -180,45 +205,45 @@ def step( if use_polytropic: # COPV temperature: polytropic expansion/compression - if step._m_copv_0 > 1e-10 and params.V_copv > 1e-10: - rho_copv_0 = step._m_copv_0 / params.V_copv + if params._m_copv_0 > 1e-10 and params.V_copv > 1e-10: + rho_copv_0 = params._m_copv_0 / params.V_copv rho_copv = m_gas_copv / params.V_copv if params.V_copv > 1e-10 else rho_copv_0 if rho_copv_0 > 1e-10: - T_copv = step._T_copv_0 * (rho_copv / rho_copv_0) ** (n_poly - 1.0) + T_copv = params._T_copv_0 * (rho_copv / rho_copv_0) ** (n_poly - 1.0) T_copv = max(200.0, min(400.0, T_copv)) # Clamp to reasonable range [200-400 K] else: - T_copv = step._T_copv_0 + T_copv = params._T_copv_0 else: - T_copv = step._T_copv_0 + T_copv = params._T_copv_0 # Fuel tank temperature: polytropic expansion/compression - if step._m_F_0 > 1e-10 and step._V_F_0 > 1e-10: - rho_F_0 = step._m_F_0 / step._V_F_0 + if params._m_F_0 > 1e-10 and params._V_F_0 > 1e-10: + rho_F_0 = params._m_F_0 / params._V_F_0 rho_F = m_gas_F / V_u_F if V_u_F > 1e-10 else rho_F_0 if rho_F_0 > 1e-10: - T_gas_F = step._T_F_0 * (rho_F / rho_F_0) ** (n_poly - 1.0) + T_gas_F = params._T_F_0 * (rho_F / rho_F_0) ** (n_poly - 1.0) T_gas_F = max(200.0, min(400.0, T_gas_F)) # Clamp to reasonable range else: - T_gas_F = step._T_F_0 + T_gas_F = params._T_F_0 else: - T_gas_F = step._T_F_0 + T_gas_F = params._T_F_0 # Oxidizer tank temperature: polytropic expansion/compression - if step._m_O_0 > 1e-10 and step._V_O_0 > 1e-10: - rho_O_0 = step._m_O_0 / step._V_O_0 + if params._m_O_0 > 1e-10 and params._V_O_0 > 1e-10: + rho_O_0 = params._m_O_0 / params._V_O_0 rho_O = m_gas_O / V_u_O if V_u_O > 1e-10 else rho_O_0 if rho_O_0 > 1e-10: - T_gas_O = step._T_O_0 * (rho_O / rho_O_0) ** (n_poly - 1.0) + T_gas_O = params._T_O_0 * (rho_O / rho_O_0) ** (n_poly - 1.0) T_gas_O = max(200.0, min(400.0, T_gas_O)) # Clamp to reasonable range else: - T_gas_O = step._T_O_0 + T_gas_O = params._T_O_0 else: - T_gas_O = step._T_O_0 + T_gas_O = params._T_O_0 else: # Isothermal process: temperature constant - T_copv = step._T_copv_0 - T_gas_F = step._T_F_0 - T_gas_O = step._T_O_0 + T_copv = params._T_copv_0 + T_gas_F = params._T_F_0 + T_gas_O = params._T_O_0 # Extract control u_F = np.clip(u[IDX_U_F], 0.0, 1.0) @@ -228,15 +253,16 @@ def step( # Flow path: COPV -> Regulator -> Tanks (when valves open) # Model: Compressible gas flow through orifices with proper choked/subsonic flow - # Physical constants for N2 gas - gamma_gas = 1.4 # Specific heat ratio for N2 - Cd_regulator = 0.7 # Discharge coefficient for regulator orifice - Cd_valve = 0.65 # Discharge coefficient for solenoid valve orifice - - # Effective flow areas [m²] - typical solenoid valve characteristics - A_regulator = 2e-5 # Regulator orifice area (~20 mm²) - A_valve_F = 5e-5 # Fuel solenoid valve flow area (~50 mm²) - A_valve_O = 5e-5 # Oxidizer solenoid valve flow area (~50 mm²) + # Pressurant properties and orifice geometry now come from DynamicsParams rather than + # being literals here -- same defaults, but settable per vehicle. See docs/adr/0001: this + # COPV -> regulator -> tank chain is one of the feed models lib/feedtwin absorbs. + gamma_gas = params.gamma_gas + Cd_regulator = params.Cd_regulator + Cd_valve = params.Cd_valve + + A_regulator = params.A_regulator + A_valve_F = params.A_valve_F + A_valve_O = params.A_valve_O # Gas flow from COPV to regulator [kg/s] # CRITICAL: Flow ONLY happens when at least one valve is open (u > 0) @@ -445,11 +471,11 @@ def step( # This is the key: limited gas supply + temperature drop means pressure drops faster if params.V_copv > 1e-10: # Update COPV temperature for next step (polytropic expansion) - if use_polytropic and step._m_copv_0 > 1e-10: + if use_polytropic and params._m_copv_0 > 1e-10: rho_copv_next = m_gas_copv_next / params.V_copv - rho_copv_0 = step._m_copv_0 / params.V_copv + rho_copv_0 = params._m_copv_0 / params.V_copv if rho_copv_0 > 1e-10: - T_copv_next = step._T_copv_0 * (rho_copv_next / rho_copv_0) ** (n_poly - 1.0) + T_copv_next = params._T_copv_0 * (rho_copv_next / rho_copv_0) ** (n_poly - 1.0) T_copv_next = max(200.0, min(400.0, T_copv_next)) else: T_copv_next = T_copv @@ -649,11 +675,11 @@ def step( # - When m increases faster than V, pressure increases (pressurization) if V_u_F_next > 1e-10: # Update fuel tank temperature for next step (polytropic expansion) - if use_polytropic and step._m_F_0 > 1e-10 and step._V_F_0 > 1e-10: + if use_polytropic and params._m_F_0 > 1e-10 and params._V_F_0 > 1e-10: rho_F_next = m_gas_F_next / V_u_F_next - rho_F_0 = step._m_F_0 / step._V_F_0 + rho_F_0 = params._m_F_0 / params._V_F_0 if rho_F_0 > 1e-10: - T_gas_F_next = step._T_F_0 * (rho_F_next / rho_F_0) ** (n_poly - 1.0) + T_gas_F_next = params._T_F_0 * (rho_F_next / rho_F_0) ** (n_poly - 1.0) T_gas_F_next = max(200.0, min(400.0, T_gas_F_next)) else: T_gas_F_next = T_gas_F @@ -700,11 +726,11 @@ def step( # This is the realistic behavior: flow begins -> ullage grows -> T drops -> pressure drops instantly if V_u_O_next > 1e-10: # Update oxidizer tank temperature for next step (polytropic expansion) - if use_polytropic and step._m_O_0 > 1e-10 and step._V_O_0 > 1e-10: + if use_polytropic and params._m_O_0 > 1e-10 and params._V_O_0 > 1e-10: rho_O_next = m_gas_O_next / V_u_O_next - rho_O_0 = step._m_O_0 / step._V_O_0 + rho_O_0 = params._m_O_0 / params._V_O_0 if rho_O_0 > 1e-10: - T_gas_O_next = step._T_O_0 * (rho_O_next / rho_O_0) ** (n_poly - 1.0) + T_gas_O_next = params._T_O_0 * (rho_O_next / rho_O_0) ** (n_poly - 1.0) T_gas_O_next = max(200.0, min(400.0, T_gas_O_next)) else: T_gas_O_next = T_gas_O diff --git a/EngineDesign/tests/control/robust_ddp/test_robust_ddp_dynamics.py b/EngineDesign/tests/control/robust_ddp/test_robust_ddp_dynamics.py index 0beda47af..1e713870c 100644 --- a/EngineDesign/tests/control/robust_ddp/test_robust_ddp_dynamics.py +++ b/EngineDesign/tests/control/robust_ddp/test_robust_ddp_dynamics.py @@ -256,6 +256,65 @@ def test_state_non_negative(self): f"All states should be non-negative, got: {x_next}") +class TestPolytropicReferenceStateIsPerConfig(unittest.TestCase): + """The polytropic reference state must not leak between configs. + + It used to live on the `step` FUNCTION object, so the first trajectory ever stepped in a + process set the reference temperatures/masses/volumes for every later one, across configs. + An optimizer evaluating thousands of candidates -- the Layer X access pattern in + docs/adr/0001 -- would hand candidates 2..N the reference state of candidate 1. + """ + + @staticmethod + def _params(**over): + from engine.control.robust_ddp.dynamics import DynamicsParams + base = dict(copv_cF=1.0, copv_cO=1.0, copv_loss=1.0e3, reg_ratio=0.8, + alpha_F=1.0, alpha_O=1.0, rho_F=789.0, rho_O=1140.0, + tau_line_F=0.05, tau_line_O=0.05) + base.update(over) + return DynamicsParams(**base) + + def test_reference_state_does_not_leak_between_params_objects(self): + a = self._params() + b = self._params() + a._temp_initialized = True + a._T_copv_0 = 293.0 + self.assertFalse(getattr(b, '_temp_initialized', False), + "a second config must start with its own uninitialised reference state") + + def test_no_reference_state_is_stored_on_the_step_function(self): + from engine.control.robust_ddp import dynamics as dyn + leaked = [n for n in dir(dyn.step) + if n.startswith('_T_') or n.startswith('_m_') or n.startswith('_V_') + or n == '_temp_initialized'] + self.assertEqual(leaked, [], f"process-global state back on step(): {leaked}") + + +class TestPressurisationHardwareIsConfigurable(unittest.TestCase): + """Regulator/valve geometry must be settable, not literals inside step().""" + + def test_defaults_match_the_previously_hardcoded_values(self): + from engine.control.robust_ddp.dynamics import DynamicsParams + p = DynamicsParams(copv_cF=1.0, copv_cO=1.0, copv_loss=1.0e3, reg_ratio=0.8, + alpha_F=1.0, alpha_O=1.0, rho_F=789.0, rho_O=1140.0, + tau_line_F=0.05, tau_line_O=0.05) + self.assertAlmostEqual(p.gamma_gas, 1.4) + self.assertAlmostEqual(p.Cd_regulator, 0.7) + self.assertAlmostEqual(p.Cd_valve, 0.65) + self.assertAlmostEqual(p.A_regulator, 2e-5) + self.assertAlmostEqual(p.A_valve_F, 5e-5) + self.assertAlmostEqual(p.A_valve_O, 5e-5) + + def test_values_are_overridable(self): + from engine.control.robust_ddp.dynamics import DynamicsParams + p = DynamicsParams(copv_cF=1.0, copv_cO=1.0, copv_loss=1.0e3, reg_ratio=0.8, + alpha_F=1.0, alpha_O=1.0, rho_F=789.0, rho_O=1140.0, + tau_line_F=0.05, tau_line_O=0.05, + A_valve_F=9.9e-5, Cd_regulator=0.55) + self.assertAlmostEqual(p.A_valve_F, 9.9e-5) + self.assertAlmostEqual(p.Cd_regulator, 0.55) + + if __name__ == '__main__': unittest.main() From a1a299a9a495d159168c119355ba1c2aff9fd193 Mon Sep 17 00:00:00 2001 From: Carlsaurus Date: Tue, 8 Sep 2026 00:03:17 -0700 Subject: [PATCH 4/5] Stop a render error from blanking the whole app The frontend had NO error boundary anywhere -- grep for ErrorBoundary / componentDidCatch / getDerivedStateFromError returned nothing. In React that means any exception thrown during render unmounts the ENTIRE tree, so every render bug presents identically: a white page, no message, no stack, no way to report what happened. That is how "press Optimize on the flight page and it goes blank" arrived with nothing to act on. Adds components/ErrorBoundary.tsx and wraps each of the nine tab panels in one. The tab buttons live in
, above the panels, so a crashing tab now shows a copyable error with its component stack and a Try again button while the tab bar stays alive and every other tab keeps working. Verified end to end, not just by inspection: injected a throw at the top of FlightSimulation, reloaded, and confirmed the boundary rendered "Flight Simulation hit an error" with the stack while the rest of the app -- header, all nine tabs, the other panels -- stayed mounted and usable. Throw reverted. Note this does not by itself remove the underlying throw on the flight optimize path. I could not reproduce that one: the handler is correctly guarded (checks result.error, checks result.data, wrapped in try/catch -- though try/catch does not cover render), and every field the flight UI renders is a REQUIRED float or List[float] in the backend models, with all Optional sub-objects already guarded (results?.truncation?., results.propellant &&, and the !results?.trajectory early return). With the boundary in place the next occurrence prints the real error instead of a blank page, which is what makes it fixable. The Try again button is declared VIEW_ONLY in the checkout gating audit: it clears local error state and touches no design state. Frontend gating audit 4 passed; npm run build clean. Co-Authored-By: Claude Opus 5 --- EngineDesign/frontend/src/App.tsx | 185 ++++++++++-------- .../frontend/src/components/ErrorBoundary.tsx | 75 +++++++ EngineDesign/frontend/src/lib/gating.test.ts | 1 + 3 files changed, 178 insertions(+), 83 deletions(-) create mode 100644 EngineDesign/frontend/src/components/ErrorBoundary.tsx diff --git a/EngineDesign/frontend/src/App.tsx b/EngineDesign/frontend/src/App.tsx index a509dd39b..2feac84dd 100644 --- a/EngineDesign/frontend/src/App.tsx +++ b/EngineDesign/frontend/src/App.tsx @@ -13,6 +13,7 @@ import ConfigurationSelector from './components/ConfigurationSelector'; import { emitConfigChanged } from './lib/configBus'; import { useViewState } from './lib/viewState'; import { DesignVersions } from './components/DesignVersions'; +import { ErrorBoundary } from './components/ErrorBoundary'; import { ReadOnlyProvider } from '@stardesign-ui'; import { getConfig, getHealth } from './api/client'; import type { EngineConfig } from './api/client'; @@ -237,117 +238,135 @@ function App() { {/* Keep all tab panels mounted; hide inactive ones to preserve state */}
-
- {!config && ( -
-

Load Configuration

- -
- )} - -
+ +
+ {!config && ( +
+

Load Configuration

+ +
+ )} + +
+
-
- {!config && ( -
-

Load Configuration

- -
- )} - -
+ +
+ {!config && ( +
+

Load Configuration

+ +
+ )} + +
+
- + + +
-
- {!config && ( -
-

Load Configuration

- -
- )} - -
+ +
+ {!config && ( +
+

Load Configuration

+ +
+ )} + +
+
-
- {!config && ( -
-

Load Configuration

- -
- )} - -
+ +
+ {!config && ( +
+

Load Configuration

+ +
+ )} + +
+
-
- {!config && ( -
-

Load Configuration

- -
- )} - -
+ +
+ {!config && ( +
+

Load Configuration

+ +
+ )} + +
+
-
- {!config && ( -
-

Load Configuration

- -
- )} - -
+ +
+ {!config && ( +
+

Load Configuration

+ +
+ )} + +
+
-
- {!config && ( -
-

Load Configuration

- -
- )} - -
+ +
+ {!config && ( +
+

Load Configuration

+ +
+ )} + +
+
-
- {/* Upload section - compact */} -
-
-
- -
- {config && ( -
- - - - Config loaded and ready + +
+ {/* Upload section - compact */} +
+
+
+
- )} + {config && ( +
+ + + + Config loaded and ready +
+ )} +
-
- {/* Editor section - full width */} -
- + {/* Editor section - full width */} +
+ +
-
+
diff --git a/EngineDesign/frontend/src/components/ErrorBoundary.tsx b/EngineDesign/frontend/src/components/ErrorBoundary.tsx new file mode 100644 index 000000000..d9e2e1d34 --- /dev/null +++ b/EngineDesign/frontend/src/components/ErrorBoundary.tsx @@ -0,0 +1,75 @@ +import { Component, type ErrorInfo, type ReactNode } from 'react'; + +/** + * Catches render-time exceptions so one bad value cannot blank the whole app. + * + * Without a boundary anywhere in the tree, React unmounts EVERYTHING when a + * render throws -- the user sees a white page with no message, no stack, and no + * way to report what happened. Every render bug then looks identical, which is + * exactly how "press Optimize, page goes blank" got reported with nothing to go + * on. This keeps the failure on screen and legible instead. + */ +interface Props { + children: ReactNode; + /** Shown above the error, e.g. "Flight Simulation". */ + label?: string; +} +interface State { + error: Error | null; + info: ErrorInfo | null; +} + +export class ErrorBoundary extends Component { + state: State = { error: null, info: null }; + + static getDerivedStateFromError(error: Error): Partial { + return { error }; + } + + componentDidCatch(error: Error, info: ErrorInfo) { + // Keep the console record: the boundary stops the crash from propagating, + // so without this the stack would be swallowed entirely. + console.error('[ErrorBoundary]', this.props.label ?? '', error, info.componentStack); + this.setState({ info }); + } + + private reset = () => this.setState({ error: null, info: null }); + + render() { + const { error, info } = this.state; + if (!error) return this.props.children; + + const detail = [error.stack || String(error), info?.componentStack] + .filter(Boolean) + .join('\n\nComponent stack:'); + + return ( +
+

+ {this.props.label ? `${this.props.label} hit an error` : 'Something went wrong'} +

+

+ The rest of the app is still running. Copy the detail below when reporting this. +

+

{String(error.message || error)}

+
+ + Show stack + +
+            {detail}
+          
+
+ +
+ ); + } +} + +export default ErrorBoundary; diff --git a/EngineDesign/frontend/src/lib/gating.test.ts b/EngineDesign/frontend/src/lib/gating.test.ts index 4aa7f619c..74474c97b 100644 --- a/EngineDesign/frontend/src/lib/gating.test.ts +++ b/EngineDesign/frontend/src/lib/gating.test.ts @@ -79,6 +79,7 @@ const VIEW_ONLY: Record = { 'ConfigEditor.tsx:setSearchQuery': 'filters which sections are shown', 'ConfigEditor.tsx:setIsExpanded': 'expand/collapse a section', 'ConfigUpload.tsx:label': 'the drop zone wrapper, not a control', + 'ErrorBoundary.tsx:this.reset': 'clears a caught render error; touches no design state', 'Layer1Optimization.tsx:setShowParameterPlots': 'chart visibility', 'Layer1Optimization.tsx:setShowInjectorPressures': 'chart visibility', 'Layer1Optimization.tsx:setShowSolverInputsEcho': 'diagnostics visibility', From 5338291a49dc263144dc91964de3620153fcef0d Mon Sep 17 00:00:00 2001 From: Carlsaurus Date: Sat, 12 Sep 2026 00:48:13 -0700 Subject: [PATCH 5/5] EngineDesign overnight sweep: solver window, config-driven stability, flight dedupe, lockout heartbeat, UI cleanup Physics / solver - Chamber-pressure root-find window: choked-flow floor (2 atm) and tank-pressure ceiling (min tank minus 2%) replace the 1 bar floor and hardcoded 15% "feed loss margin"; both Python and the numba kernel now scan for the HIGHEST-Pc sign change before Brent. The old window let Brent lock onto a spurious ~20 psi root: configs/default.yaml reported F = -333 N at its own configured tank pressures, and the Python side said "no solution" at 550/650/900 psi where the kernel solved. Both paths now agree at every pressure. - Stability model reads config, not constants: StabilityConfig (n, chi, nozzle-entrance Mach, damping fractions, regulator, acoustic gate allowance), feed_system..length (the chug inertance length was a hardcoded 0.305 m with no schema field), cp_g from the CEA state (was a fixed 2200 J/kg-K), nozzle-entrance Mach solved from the contraction ratio (was 0.2), per-fluid handbook fallbacks recorded via assume() (fuel fallbacks were methane's for every fuel). comprehensive_stability_analysis looked up feed_system["lox"] .length, a key and attribute that never existed, so feed acoustics always used 1.0 m x 10 mm; the legacy and rich paths also used different chamber lengths. Chug gate centred at gain margin 1.0 (0.80 called an unstable loop "stable"). Rigged heuristics removed: stability_index floors, the water-hammer-to-margin map, and stability/enhanced.py. - Ea mixture-ratio step in reaction_chemistry smoothed (20% jump at MR 1.5); the dead reaction-progress computation on every chamber solve removed. - Layer 1: refuses a target O/F outside the propellant's CEA table (a propellant switch kept the old target and optimized against a table that could not evaluate it); names the blocking hard constraint when a run ends with no feasible candidate; 1 m chamber-length cap removed; ED_L1_WORKERS / ED_L1_TRACE_INFEAS debugging hooks. - Canonical pintle design_MR 2.55 -> 1.4 (ethalox table is [1.0, 2.5]). Flight - Optimizer Layer 4 now renders the shared FlightSimulation component; the 1100-line copy it replaced hardcoded RP-1 density and an 85% fill factor client-side, so every non-kerolox design capped its fuel load against the wrong tank. - flight router: explicit_capacity_kg carried a bool; rocket/environment fallbacks derived from the request models instead of a drifting literal copy. - tests/test_flight_propellant_iteration.py rebuilt on configs/default.yaml (it read a file out of one developer's Downloads folder); the tank-cap test now constructs its premise. Checkout lockout - useCheckout re-takes the checkout every 60 s while held and no longer releases on tab hide: the token lapsed after 5 min without a changed autosave and the next write came back 423 ("Take Design 1 before saving"). UI - Emoji removed everywhere; legend and axis-title overlap fixed on every chart with the pattern (Layer 1 convergence, pressure curves, controller, custom plotter, heat flux, flight); Design Requirements input caps that assumed one rocket removed; CEA O/F range shown under the target; stale stability prose replaced; Time-Series simple profile seeds from the design's tank pressures; Stability Model section in the config editor; unreachable Coaxial label dropped. Schema hygiene - Removed fields nothing read (hot_gas_cp, yield_strength, youngs_modulus, mixing_model) and the stale hot_gas_cp keys in shipped configs; pydantic ConfigDict; fuel tank descriptions no longer say RP-1. Co-Authored-By: Claude Fable 5.1 --- EngineDesign/backend/routers/flight.py | 48 +- EngineDesign/configs/canonical/impinging.yaml | 7 +- EngineDesign/configs/canonical/pintle.yaml | 9 +- EngineDesign/configs/default.yaml | 7 +- EngineDesign/configs/impinging_lox_ch4.yaml | 1 - .../configs/impinging_lox_ch4_8000N.yaml | 1 - .../impinging_lox_ch4_8000N_optimal.yaml | 1 - EngineDesign/configs/impinging_smoke.yaml | 1 - EngineDesign/configs/test.yaml | 1 - EngineDesign/engine/accel/kernels.py | 64 +- EngineDesign/engine/core/chamber_solver.py | 420 +++---- EngineDesign/engine/core/runner.py | 2 +- .../layers/layer1_static_optimization.py | 125 +- .../engine/pipeline/config_schemas.py | 84 +- .../pipeline/physics_based_replacements.py | 56 - .../engine/pipeline/reaction_chemistry.py | 47 +- .../engine/pipeline/stability/analysis.py | 537 ++++---- .../engine/pipeline/stability/core.py | 31 + .../engine/pipeline/stability/enhanced.py | 383 ------ .../engine/pipeline/stability/report.py | 34 +- .../engine/pipeline/time_varying_solver.py | 136 +- EngineDesign/frontend/src/api/client.ts | 1 + .../frontend/src/components/ConfigEditor.tsx | 60 +- .../src/components/ConfigurationSelector.tsx | 1 - .../src/components/ControllerMode.tsx | 32 +- .../frontend/src/components/CustomPlotter.tsx | 4 +- .../src/components/DesignRequirements.tsx | 112 +- .../src/components/FlightSimulation.tsx | 22 +- .../src/components/HeatFluxProfileChart.tsx | 8 +- .../src/components/Layer1Optimization.tsx | 41 +- .../src/components/Layer2Optimization.tsx | 34 +- .../src/components/Layer3Optimization.tsx | 40 +- .../src/components/Layer4Optimization.tsx | 1093 ----------------- .../frontend/src/components/Optimizer.tsx | 30 +- .../src/components/PressureCurveChart.tsx | 44 +- .../src/components/ResultsDisplay.tsx | 7 +- .../src/components/TimeSeriesMode.tsx | 23 +- EngineDesign/frontend/src/lib/gating.test.ts | 2 - .../tests/test_flight_propellant_iteration.py | 146 +-- .../tests/test_layer1_of_target_range.py | 39 + lib/stardesign-ui/src/useCheckout.ts | 49 +- 41 files changed, 1197 insertions(+), 2586 deletions(-) delete mode 100644 EngineDesign/engine/pipeline/stability/enhanced.py delete mode 100644 EngineDesign/frontend/src/components/Layer4Optimization.tsx create mode 100644 EngineDesign/tests/test_layer1_of_target_range.py diff --git a/EngineDesign/backend/routers/flight.py b/EngineDesign/backend/routers/flight.py index d9785b83d..744d1129f 100644 --- a/EngineDesign/backend/routers/flight.py +++ b/EngineDesign/backend/routers/flight.py @@ -342,14 +342,8 @@ def build_flight_config(base_config, request: FlightSimRequest): config_dict["environment"]["date"] = request.environment.date config_dict["environment"]["atmosphere_model"] = request.environment.atmosphere_model elif config_dict.get("environment") is None: - # Set defaults - config_dict["environment"] = { - "latitude": 35.0, - "longitude": -117.0, - "elevation": 0.0, - "date": [2025, 1, 1, 12], - "atmosphere_model": "standard_atmosphere", - } + # No launch site anywhere: the request model's defaults, stated once, up in EnvironmentConfig. + config_dict["environment"] = EnvironmentConfig().model_dump() # Update rocket if request.rocket: @@ -379,20 +373,20 @@ def build_flight_config(base_config, request: FlightSimRequest): "fin_position": request.rocket.fins.fin_position, } elif config_dict.get("rocket") is None: - # Set defaults + # No vehicle anywhere: derive from the request model's defaults so this block cannot drift + # from RocketConfig (the literal copy that used to sit here said propulsion_dry_mass 24 kg + # while its own component defaults summed to 16). + rd = RocketConfig() config_dict["rocket"] = { - "airframe_mass": 78.72, - "propulsion_dry_mass": 24.0, - "radius": 0.1015, - "motor_position": 0.0, - "inertia": [8.0, 8.0, 0.5], - "fins": { - "no_fins": 3, - "root_chord": 0.2, - "tip_chord": 0.1, - "fin_span": 0.3, - "fin_position": 0.1, - }, + "airframe_mass": rd.airframe_mass, + "propulsion_dry_mass": rd.engine_mass + rd.lox_tank_structure_mass + rd.fuel_tank_structure_mass, + "radius": rd.radius, + "motor_position": rd.motor_position, + "inertia": list(rd.inertia), + "nose_kind": rd.nose_kind, + "nose_fineness_ratio": rd.nose_fineness_ratio, + "avionics_payload_length_m": rd.avionics_payload_length_m, + "fins": FinsConfig().model_dump(), } return config_dict @@ -424,19 +418,19 @@ def _apply_propellant_mass_caps(config_dict: dict, base_config) -> tuple[dict, d lox_tank_max = lox_max fill_factor = lox_ff current_lox = float(config_dict.get("lox_tank", {}).get("mass", 0) or 0) - effective = min(current_lox, lox_max) if current_lox > lox_max else current_lox + effective = min(current_lox, lox_max) if current_lox > lox_max: config_dict["lox_tank"]["mass"] = lox_max cap_note = "explicit capacity" if lox_explicit else f"{lox_ff * 100:.0f}% fill" print(f"[Flight] Capped LOX mass: {current_lox:.2f} -> {lox_max:.2f} kg ({cap_note}, vol {lox_vol * 1000:.1f}L)") mass_adjustments["lox"] = { "original": current_lox, - "capped": effective if current_lox <= lox_max else lox_max, + "capped": effective, "max_fill_kg": lox_max, "tank_volume_m3": lox_vol, "fill_factor": lox_ff, "was_capped": current_lox > lox_max + 1e-6, - "explicit_capacity_kg": lox_explicit, + "explicit_capacity_kg": lox_max if lox_explicit else None, } if cap_config.fuel_tank is not None: @@ -444,19 +438,19 @@ def _apply_propellant_mass_caps(config_dict: dict, base_config) -> tuple[dict, d fuel_tank_max = fuel_max fill_factor = fuel_ff current_fuel = float(config_dict.get("fuel_tank", {}).get("mass", 0) or 0) - effective = min(current_fuel, fuel_max) if current_fuel > fuel_max else current_fuel + effective = min(current_fuel, fuel_max) if current_fuel > fuel_max: config_dict["fuel_tank"]["mass"] = fuel_max cap_note = "explicit capacity" if fuel_explicit else f"{fuel_ff * 100:.0f}% fill" print(f"[Flight] Capped Fuel mass: {current_fuel:.2f} -> {fuel_max:.2f} kg ({cap_note}, vol {fuel_vol * 1000:.1f}L)") mass_adjustments["fuel"] = { "original": current_fuel, - "capped": effective if current_fuel <= fuel_max else fuel_max, + "capped": effective, "max_fill_kg": fuel_max, "tank_volume_m3": fuel_vol, "fill_factor": fuel_ff, "was_capped": current_fuel > fuel_max + 1e-6, - "explicit_capacity_kg": fuel_explicit, + "explicit_capacity_kg": fuel_max if fuel_explicit else None, } return mass_adjustments, lox_tank_max, fuel_tank_max, fill_factor diff --git a/EngineDesign/configs/canonical/impinging.yaml b/EngineDesign/configs/canonical/impinging.yaml index 70b082a3c..a252b376e 100644 --- a/EngineDesign/configs/canonical/impinging.yaml +++ b/EngineDesign/configs/canonical/impinging.yaml @@ -28,11 +28,17 @@ feed_system: # under-predicted LOX feed loss by ~3.8x (dP ~ 1/A^2). fuel: line_size: 3/8_NPT + # Tank outlet to injector manifold. Sets the chug-model line inertance (length/area); + # 0.305 m is the value the model used to carry hardcoded -- measure it on the vehicle. + length: 0.305 K0: 2.0 K1: 0.0 phi_type: none oxidizer: line_size: 3/8_NPT + # Tank outlet to injector manifold. Sets the chug-model line inertance (length/area); + # 0.305 m is the value the model used to carry hardcoded -- measure it on the vehicle. + length: 0.305 K0: 2.0 K1: 0.0 phi_type: none @@ -67,7 +73,6 @@ regen_cooling: n_segments: 20 gas_turbulence_intensity: 0.1 coolant_turbulence_intensity: 0.05 - hot_gas_cp: 2200.0 recovery_factor: null film_cooling: enabled: false diff --git a/EngineDesign/configs/canonical/pintle.yaml b/EngineDesign/configs/canonical/pintle.yaml index 2189f871d..b144c2425 100644 --- a/EngineDesign/configs/canonical/pintle.yaml +++ b/EngineDesign/configs/canonical/pintle.yaml @@ -33,7 +33,7 @@ chamber_geometry: Cf: 1.5422822280584674 Lstar: 1.239905396733828 chamber_diameter: 0.11344154349849075 - design_MR: 2.55 + design_MR: 1.4 design_pressure: 2413166.0 design_thrust: 7000.0 exit_diameter: 0.101 @@ -173,11 +173,17 @@ feed_system: # under-predicted LOX feed loss by ~3.8x (dP ~ 1/A^2). fuel: line_size: 3/8_NPT + # Tank outlet to injector manifold. Sets the chug-model line inertance (length/area); + # 0.305 m is the value the model used to carry hardcoded -- measure it on the vehicle. + length: 0.305 K0: 2.0 K1: 0.0 phi_type: none oxidizer: line_size: 3/8_NPT + # Tank outlet to injector manifold. Sets the chug-model line inertance (length/area); + # 0.305 m is the value the model used to carry hardcoded -- measure it on the vehicle. + length: 0.305 K0: 2.0 K1: 0.0 phi_type: none @@ -421,7 +427,6 @@ regen_cooling: d_outlet: null enabled: false gas_turbulence_intensity: 0.1 - hot_gas_cp: 2200.0 hot_gas_prandtl: 0.7 hot_gas_thermal_conductivity: 0.12 hot_gas_viscosity: 4.0e-05 diff --git a/EngineDesign/configs/default.yaml b/EngineDesign/configs/default.yaml index 3487be943..8d74feb69 100644 --- a/EngineDesign/configs/default.yaml +++ b/EngineDesign/configs/default.yaml @@ -45,11 +45,17 @@ feed_system: # under-predicted LOX feed loss by ~3.8x (dP ~ 1/A^2). fuel: line_size: 3/8_NPT + # Tank outlet to injector manifold. Sets the chug-model line inertance (length/area); + # 0.305 m is the value the model used to carry hardcoded -- measure it on the vehicle. + length: 0.305 K0: 2.0 K1: 0.0 phi_type: none oxidizer: line_size: 3/8_NPT + # Tank outlet to injector manifold. Sets the chug-model line inertance (length/area); + # 0.305 m is the value the model used to carry hardcoded -- measure it on the vehicle. + length: 0.305 K0: 2.0 K1: 0.0 phi_type: none @@ -84,7 +90,6 @@ regen_cooling: n_segments: 20 gas_turbulence_intensity: 0.1 coolant_turbulence_intensity: 0.05 - hot_gas_cp: 2200.0 recovery_factor: null film_cooling: enabled: false diff --git a/EngineDesign/configs/impinging_lox_ch4.yaml b/EngineDesign/configs/impinging_lox_ch4.yaml index 7fc8bc47b..1d3e8c37d 100644 --- a/EngineDesign/configs/impinging_lox_ch4.yaml +++ b/EngineDesign/configs/impinging_lox_ch4.yaml @@ -400,7 +400,6 @@ regen_cooling: d_outlet: null enabled: false gas_turbulence_intensity: 0.1 - hot_gas_cp: 2200.0 hot_gas_prandtl: 0.7 hot_gas_thermal_conductivity: 0.12 hot_gas_viscosity: 4.0e-05 diff --git a/EngineDesign/configs/impinging_lox_ch4_8000N.yaml b/EngineDesign/configs/impinging_lox_ch4_8000N.yaml index 37970ae5d..624a2452e 100644 --- a/EngineDesign/configs/impinging_lox_ch4_8000N.yaml +++ b/EngineDesign/configs/impinging_lox_ch4_8000N.yaml @@ -443,7 +443,6 @@ regen_cooling: d_outlet: null enabled: false gas_turbulence_intensity: 0.1 - hot_gas_cp: 2200.0 hot_gas_prandtl: 0.7 hot_gas_thermal_conductivity: 0.12 hot_gas_viscosity: 4.0e-05 diff --git a/EngineDesign/configs/impinging_lox_ch4_8000N_optimal.yaml b/EngineDesign/configs/impinging_lox_ch4_8000N_optimal.yaml index 17d65aa1b..2f5b726cd 100644 --- a/EngineDesign/configs/impinging_lox_ch4_8000N_optimal.yaml +++ b/EngineDesign/configs/impinging_lox_ch4_8000N_optimal.yaml @@ -82,7 +82,6 @@ regen_cooling: n_segments: 20 gas_turbulence_intensity: 0.1 coolant_turbulence_intensity: 0.05 - hot_gas_cp: 2200.0 recovery_factor: null film_cooling: enabled: false diff --git a/EngineDesign/configs/impinging_smoke.yaml b/EngineDesign/configs/impinging_smoke.yaml index 4919ae37a..a4217e861 100644 --- a/EngineDesign/configs/impinging_smoke.yaml +++ b/EngineDesign/configs/impinging_smoke.yaml @@ -56,7 +56,6 @@ regen_cooling: n_segments: 20 gas_turbulence_intensity: 0.1 coolant_turbulence_intensity: 0.05 - hot_gas_cp: 2200.0 recovery_factor: null film_cooling: enabled: false diff --git a/EngineDesign/configs/test.yaml b/EngineDesign/configs/test.yaml index a8cdd9aeb..dd029b8ec 100644 --- a/EngineDesign/configs/test.yaml +++ b/EngineDesign/configs/test.yaml @@ -76,7 +76,6 @@ regen_cooling: n_segments: 20 gas_turbulence_intensity: 0.1 coolant_turbulence_intensity: 0.05 - hot_gas_cp: 2200.0 recovery_factor: null film_cooling: enabled: false diff --git a/EngineDesign/engine/accel/kernels.py b/EngineDesign/engine/accel/kernels.py index 90ec1ce63..e484fe8eb 100644 --- a/EngineDesign/engine/accel/kernels.py +++ b/EngineDesign/engine/accel/kernels.py @@ -966,8 +966,11 @@ def _solve_exit_mach(eps, g): @njit(cache=True) def evaluate_core(P, Pcg, MRg, epsg, cstar, Cf, Tc_t, gam, Rt, Mt, Cfv, P_O, P_F, Pa): """Returns (ok, Pc, F, Isp, MR, cstar_actual, gamma, Tc, mdot_total, v_exit, Cf_actual).""" - Pc_min = 100000.0 - Pc_max = min(P_O, P_F)*(1.0 - 0.15) + # Search window and bracket scan: VERBATIM mirror of engine/core/chamber_solver.py + # (PC_CHOKE_FLOOR_PA, PC_MIN_TOTAL_DROP_FRAC, ChamberSolver._highest_sign_change). The parity + # gate requires the two root-finds to agree, so change both or neither. + Pc_min = 2.0*101325.0 + Pc_max = min(P_O, P_F)*(1.0 - 0.02) Pc_min = max(Pc_min, P[SV_PCMIN]); Pc_max = min(Pc_max, P[SV_PCMAX]) if Pc_max <= Pc_min: return (0.0,)*22 @@ -976,41 +979,28 @@ def evaluate_core(P, Pcg, MRg, epsg, cstar, Cf, Tc_t, gam, Rt, Mt, Cfv, P_O, P_F rmax = _residual(Pc_max, P, Pcg, MRg, epsg, cstar, Cf, Tc_t, gam, Rt, Mt, Cfv, P_O, P_F) if not np.isfinite(rmin) or not np.isfinite(rmax): return (0.0,)*22 - if _sign(rmin) == _sign(rmax): - # The residual is not always monotonic in Pc: on the canonical impinging design at - # asymmetric tank pressures it is NEGATIVE at Pc_min, turns POSITIVE mid-range, and goes - # negative again by Pc_max, so an endpoint-only sign test sees "no root" and bails where - # a root plainly exists (measured: -0.214 at 1.0 bar, +0.76 at 20 bar, -0.96 at 27.4 bar, - # root at 23.75 bar). Scan for a sign-changing sub-interval before giving up -- the pure - # Python solver finds these, and the parity gate requires the accelerator to agree. - # Scan from the HIGH end down: a non-monotonic residual can also cross near Pc_min - # (a spurious low-pressure crossing where the chamber is barely flowing). The physical - # operating point is the HIGHEST-Pc root, which is the one the pure Python solver - # converges to; taking the first crossing from the bottom picked the spurious one and - # diverged from Python by 3.3x. - _n = 32 - _lo = 0.0; _hi = 0.0; _found = False - _pb = Pc_max; _rb = rmax - for _i in range(_n - 1, -1, -1): - _pa = Pc_min + (Pc_max - Pc_min)*(_i/_n) - _ra = _residual(_pa, P, Pcg, MRg, epsg, cstar, Cf, Tc_t, gam, Rt, Mt, Cfv, P_O, P_F) - if np.isfinite(_ra) and np.isfinite(_rb) and _sign(_ra) != _sign(_rb): - _lo = _pa; _hi = _pb; _found = True - break - _pb = _pa; _rb = _ra - if _found: - Pc = _brentq(P, Pcg, MRg, epsg, cstar, Cf, Tc_t, gam, Rt, Mt, Cfv, P_O, P_F, - _lo, _hi, xtol, rtol, maxit) - if not np.isfinite(Pc): - return (0.0,)*22 - elif rmin > 0 and rmax > 0 and rmax < 0.1: - Pc = Pc_max - else: - return (0.0,)*22 - else: - Pc = _brentq(P, Pcg, MRg, epsg, cstar, Cf, Tc_t, gam, Rt, Mt, Cfv, P_O, P_F, Pc_min, Pc_max, xtol, rtol, maxit) - if not np.isfinite(Pc): - return (0.0,)*22 + # The residual is not monotonic in Pc: it is negative near the choking floor (barely flowing, + # the spurious crossing), positive through the operating range, and negative again once the + # injector drop gets small. The physical operating point is the HIGHEST-Pc root, so scan from + # Pc_max down for the first sign change and bracket Brent inside it. Endpoint-only tests + # missed interior roots (measured on canonical/impinging: -0.214 at 1 bar, +0.76 at 20 bar, + # -0.96 at 27.4 bar, root at 23.75 bar) and a whole-window Brent picked the spurious low one. + _n = 32 + _lo = 0.0; _hi = 0.0; _found = False + _pb = Pc_max; _rb = rmax + for _i in range(_n - 1, -1, -1): + _pa = Pc_min + (Pc_max - Pc_min)*(_i/_n) + _ra = rmin if _i == 0 else _residual(_pa, P, Pcg, MRg, epsg, cstar, Cf, Tc_t, gam, Rt, Mt, Cfv, P_O, P_F) + if np.isfinite(_ra) and np.isfinite(_rb) and _sign(_ra) != _sign(_rb): + _lo = _pa; _hi = _pb; _found = True + break + _pb = _pa; _rb = _ra + if not _found: + return (0.0,)*22 + Pc = _brentq(P, Pcg, MRg, epsg, cstar, Cf, Tc_t, gam, Rt, Mt, Cfv, P_O, P_F, + _lo, _hi, xtol, rtol, maxit) + if not np.isfinite(Pc): + return (0.0,)*22 # recompute converged state (ok, mO, mF, uO, uF, D32O, D32F, momR, CdO, CdF, PiO, PiF, dpiO, dpiF, AgO, AgF, dpfO, dpfF, WeO, WeF, urel, xstar, constr, nit, tiO, tiF) = _solve_injector(P, P_O, P_F, Pc) diff --git a/EngineDesign/engine/core/chamber_solver.py b/EngineDesign/engine/core/chamber_solver.py index a83e11188..8a25a61b3 100644 --- a/EngineDesign/engine/core/chamber_solver.py +++ b/EngineDesign/engine/core/chamber_solver.py @@ -42,6 +42,23 @@ from engine.core.closure import flows +# 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 + + class ChamberSolver: """Solves for chamber pressure by balancing supply and demand""" @@ -258,6 +275,88 @@ def _accel_chamber_pc(self, P_tank_O: float, P_tank_F: float): 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(), + } + 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 + mdot_demand_test = (Pc_max * 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 + 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, @@ -285,24 +384,10 @@ def solve( diagnostics : dict Solution diagnostics """ - # Determine bounds - # Realistic bounds: Pc must be less than both tank pressures (accounting for feed losses) - Pc_min = 100000.0 # 100 kPa minimum - - # Estimate maximum feed losses for bounds calculation - # Use rough estimates: assume maximum flow gives ~10-20% pressure drop - # This is conservative but better than fixed 5% margin - # Actual feed losses will be calculated during solve - feed_loss_margin = 0.15 # 15% margin for feed losses (conservative estimate) - Pc_max = min(P_tank_O, P_tank_F) * (1.0 - feed_loss_margin) - - # If fuel pressure is much higher than oxidizer, we might need to allow - # Pc up to oxidizer pressure (since oxidizer flow limits the system) - # But this is already handled by min(P_tank_O, P_tank_F) - - # Clamp to config bounds - Pc_min = max(Pc_min, self.config.solver.Pc_bounds[0]) - Pc_max = min(Pc_max, self.config.solver.Pc_bounds[1]) + # 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})") @@ -325,259 +410,94 @@ def solve( def residual_func(Pc): return self.residual(Pc, P_tank_O, P_tank_F) - # Native fast path: run the whole residual loop + Brent in C when native can - # handle this config. On success skip the Python root-find. Importantly, this - # path is hit by the ~30% of Layer-1 CMA candidates that don't converge in the - # single-call ed_evaluate seam and fall back to runner.evaluate -> here, so - # keeping it native keeps that fallback fast (a pure-Python Brent solve here is - # ~100x slower). Any failure -> Python Brent below. + # 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 - skip_solve = True - residual_min, residual_max = -1.0, 1.0 else: - # Check residual signs at bounds before solving residual_min = residual_func(Pc_min) residual_max = residual_func(Pc_max) - - # Check for NaN values and provide better error messages - if not np.isfinite(residual_min): - # Try to diagnose the issue - try: - # Test a few points to see where it fails - test_Pc = (Pc_min + Pc_max) / 2 - test_res = residual_func(test_Pc) - if not np.isfinite(test_res): + + 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 returns non-finite values. " + f"Residual function evaluation failed at bounds. " f"Pc_min={Pc_min/1e6:.2f} MPa, Pc_max={Pc_max/1e6:.2f} MPa. " - f"Check injector geometry, feed system, or CEA cache." + f"Error: {e}" ) - except Exception as e: + if not np.isfinite(residual_max): 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}" + 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." ) - - 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." - ) - - # brentq requires opposite signs at bounds - if np.sign(residual_min) == np.sign(residual_max): - # No root in interval - this happens when: - # 1. Supply > demand at all Pc (both positive) - need higher Pc but limited by tank pressure - # 2. Supply < demand at all Pc (both negative) - can't supply enough flow - - if residual_min > 0 and residual_max > 0: - # Supply > Demand at all Pc - # This means injector supplies more flow than combustion can demand - # Common causes: - # 1. Injector too large (orifice areas too big) - # 2. Throat too small (can't flow enough to balance supply) - # 3. Combustion efficiency too low (reduces demand) - # 4. Pc_max too conservative (we could go slightly higher) - - # Initialize skip_solve flag - skip_solve = False - - # Check if residual is small at Pc_max (near solution) - residual_tolerance = 0.1 # kg/s - accept if within 0.1 kg/s - - if residual_max < residual_tolerance: - # Residual is small - we're very close to solution - # Use Pc_max as solution with warning - # import warnings - # warnings.warn( - # f"Supply slightly > Demand at Pc_max. " - # f"Using Pc_max ({Pc_max/1e6:.2f} MPa) as solution. " - # f"Residual: {residual_max:.4f} kg/s. " - # f"Injector may be slightly oversized or throat slightly undersized." - # ) - # Skip to solution validation - use Pc_max as solution - Pc = Pc_max - success = True - # Skip the root finding loop below - skip_solve = True - else: - # Residual is significant - diagnose the issue - # Get diagnostics at Pc_max to understand supply/demand - 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 - - # Get demand at Pc_max - MR_test = mdot_O_test / mdot_F_test if mdot_F_test > 0 else np.inf - cg = ensure_chamber_geometry(self.config) - eps_default = cg.expansion_ratio - cea_props_test = self.cea_cache.eval(MR_test, Pc_max, 101325.0, eps_default) - cstar_ideal_test = cea_props_test.get("cstar_ideal", 0.0) - - # Build advanced_params for diagnostics - 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(), - } - - # Calculate efficiency - 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 if 'debug' in locals() else False, - ) - cstar_actual_test = eta_test * cstar_ideal_test - cg = ensure_chamber_geometry(self.config) - mdot_demand_test = (Pc_max * cg.A_throat) / cstar_actual_test if cstar_actual_test > 0 else np.inf - - # Calculate what Pc would balance (extrapolate) - # residual = supply - demand - # At Pc_max: residual = mdot_supply - mdot_demand - # Demand scales with Pc: mdot_demand ∝ Pc - # Supply decreases slightly with Pc: mdot_supply decreases as Pc increases - # Rough estimate: if we increase Pc by ΔPc, demand increases more than supply - - # Estimate required Pc (rough extrapolation) - # Assume linear relationship near Pc_max - if mdot_demand_test > 0 and mdot_supply_test > mdot_demand_test: - # We need more Pc to increase demand - # mdot_demand = Pc * At / c*, so Pc_needed = mdot_supply * c* / At - cg = ensure_chamber_geometry(self.config) - Pc_estimate = mdot_supply_test * cstar_actual_test / cg.A_throat - - 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)" - ) - else: - 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: - # Re-raise explicit ValueErrors from above - raise - except Exception as diag_e: - # Diagnostics failed - provide generic error - 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}" - ) - - else: - # Supply < Demand at all Pc (both negative) + + # 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." ) - - # Check if we already have a solution (from small residual case above) - # skip_solve is defined in the if-else block above, default to False if not set - if 'skip_solve' not in locals(): - skip_solve = False - - if not skip_solve: - # Validate bracket before solving - bracket_check = NumericalStability.check_bracket(residual_func, Pc_min, Pc_max) - if not bracket_check.passed: - raise ValueError(f"Invalid bracket for root finding: {bracket_check.message}") - - # Track convergence history for diagnostics - convergence_history = [] - - # Enhanced residual function with convergence tracking + def tracked_residual_func(Pc): res = residual_func(Pc) convergence_history.append(float(res)) return res - - # Solve using bracketed secant (brentq) - safe and robust + try: if self.config.solver.method == "brentq": Pc, result = brentq( tracked_residual_func, - Pc_min, - Pc_max, + 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 - - # Validate convergence - 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}") - else: - # Fallback to Newton's method (less robust) + # Newton from the middle of the bracket (less robust; kept for configs that ask) Pc = newton( tracked_residual_func, - Pc_guess, + 0.5 * (lo + hi), tol=self.config.solver.tolerance, maxiter=self.config.solver.max_iterations ) success = True - - # Validate convergence for Newton - 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 as e: - # Re-raise ValueError (bracket issues, etc.) + 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}") - else: - # We're using Pc_max as solution (small residual case) - # Already set Pc = Pc_max and success = True above - convergence_history = [residual_max] # Store for diagnostics - + # Validate solution Pc_val = float(Pc) if not np.isfinite(Pc_val): @@ -623,46 +543,6 @@ def tracked_residual_func(Pc): with_profile=not getattr(self, "_silent", False), ) - # Calculate reaction progress through chamber (if finite-rate chemistry enabled) - reaction_progress = None - if getattr(self.config.combustion.efficiency, 'use_finite_rate_chemistry', True): - try: - from engine.pipeline.reaction_chemistry import calculate_chamber_reaction_progress - - # Pass spray diagnostics if available for better evaporation/mixing estimates - spray_diagnostics = closure_diag if closure_diag else None - - # Use conservative "Worst of Both Worlds" temperatures: - # Tc (Ideal) for residence time (shorter time is conservative) - # effective_Tc (Actual) for kinetics (slower chemistry is conservative) - reaction_progress = calculate_chamber_reaction_progress( - current_Lstar, - Pc_val, - cea_props["Tc"], # Ideal Tc (Residence Time) - cea_props["cstar_ideal"], - cea_props["gamma"], - cea_props["R"], - MR, - self.config, - spray_diagnostics=spray_diagnostics, - Tc_kinetics=effective_Tc, # Actual Tc (Kinetics) - ) - except Exception as e: - # Don't silently fail - raise error or log warning - import warnings - warnings.warn(f"Reaction progress calculation failed: {e}. This may indicate invalid engine conditions.") - # Minimal fallback - but indicate uncertainty - # CRITICAL FIX: Correct residence time formula - rho_chamber = Pc_val / (cea_props["R"] * cea_props["Tc"]) if cea_props["R"] > 0 and cea_props["Tc"] > 0 else 1.0 - # Use actual mdot_total from closure (calculated above) - cg = ensure_chamber_geometry(self.config) - tau_residence_correct = current_Lstar * rho_chamber * cg.A_throat / mdot_total if mdot_total > 0 else 0.001 - reaction_progress = { - "progress_throat": 1.0, # Assume equilibrium - "tau_residence": tau_residence_correct, - "calculation_failed": True, - } - # Extract and validate mixture diagnostics (diagnostics-only, no efficiency impact) # Enable mixture coupling diagnostics if configured eff_cfg = self.config.combustion.efficiency diff --git a/EngineDesign/engine/core/runner.py b/EngineDesign/engine/core/runner.py index a35fba084..1ec113081 100644 --- a/EngineDesign/engine/core/runner.py +++ b/EngineDesign/engine/core/runner.py @@ -472,7 +472,7 @@ def log_info(msg): "stability_state": "unstable", "stability_score": 0.0, "is_stable": False, - "chugging": {"frequency": 0.0, "stability_margin": 0.0, "stability_index": 0.0, "period": 0.0, "tau_residence": 0.0, "Lstar": 0.0}, + "chugging": {"frequency": 0.0, "stability_margin": 0.0, "period": 0.0, "tau_residence": 0.0, "Lstar": 0.0}, "acoustic": {"stability_margin": 0.0, "modes": {}, "longitudinal_modes": [], "transverse_modes": [], "sound_speed": 0.0}, "feed_system": {"pogo_frequency": 0.0, "surge_frequency": 0.0, "water_hammer_margin": 0.0, "stability_margin": 0.0, "sound_speed": 0.0}, "mode_coupling": [], diff --git a/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py b/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py index 33f3511ac..c1d84ed70 100644 --- a/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py +++ b/EngineDesign/engine/optimizer/layers/layer1_static_optimization.py @@ -1664,7 +1664,9 @@ def _layer1_apply_chamber_geometry_to_config( L_chamber = V_chamber / A_chamber if A_chamber > 0 else 0.2 L_cylindrical = max(L_chamber * 0.5, 0.05) - L_chamber = np.clip(L_chamber, 0.005, 1.0) + # Positivity guard only. The old clip also capped the chamber at 1.0 m, a size limit no + # requirement asked for; engine length is constrained by max_engine_length elsewhere. + L_chamber = max(float(L_chamber), 0.005) if config.chamber_geometry is None: cg = ensure_chamber_geometry(config) @@ -2684,7 +2686,7 @@ def _config_to_dict(config: PintleEngineConfig) -> dict: Uses pydantic's dict() method if available, otherwise falls back to __dict__. """ - return config.dict() if hasattr(config, 'dict') else config.__dict__ + return config.model_dump() if hasattr(config, 'model_dump') else config.__dict__ def _dict_to_config(config_dict: dict) -> PintleEngineConfig: @@ -2841,7 +2843,16 @@ def _snap_integer_dims(x: np.ndarray, integer_indices: list) -> np.ndarray: def _get_num_workers(config_obj) -> int: - """Get number of workers from config or default to cpu_count - 1.""" + """Get number of workers from config or default to cpu_count - 1. + + ``ED_L1_WORKERS`` overrides both (``1`` runs every candidate in-process, which is what a + debugger or an infeasibility trace needs).""" + env = os.environ.get("ED_L1_WORKERS") + if env: + try: + return max(1, int(env)) + except ValueError: + pass if hasattr(config_obj, 'optimizer') and hasattr(config_obj.optimizer, 'num_workers'): num_workers = config_obj.optimizer.num_workers else: @@ -2986,6 +2997,18 @@ def _apply_x_to_worker_config_inplace(x: np.ndarray, config: PintleEngineConfig, config.fuel_tank.initial_pressure_psi = _pf +# Infeasibility trace: set ED_L1_TRACE_INFEAS=1 (with ED_L1_WORKERS=1 so candidates run in-process) +# and each objective evaluation appends {checkpoint: running infeasibility score} here -- the only +# way to see WHICH gate keeps a run infeasible when every candidate returns the 1e6 floor. +_INFEAS_TRACE: List[Dict[str, float]] = [] +_INFEAS_TRACE_ON = bool(os.environ.get("ED_L1_TRACE_INFEAS")) + + +def _infeas_trace(entry: Optional[Dict[str, float]], label: str, value: float) -> None: + if entry is not None: + entry[label] = float(value) + + def _compute_objective_value(result: dict, x: np.ndarray, requirements: dict, constants: dict) -> float: """Compute objective value from evaluation result. @@ -3096,6 +3119,9 @@ def _compute_objective_value(result: dict, x: np.ndarray, requirements: dict, co P_F_ratio = P_F_psi / max_fuel_P_psi if max_fuel_P_psi > 0 else 0.0 infeasibility_score = 0.0 + _tr = {} if _INFEAS_TRACE_ON else None + if _tr is not None: + _INFEAS_TRACE.append(_tr) if A_chamber_check > 0 and A_throat_check > 0: contraction_ratio_check = A_chamber_check / A_throat_check @@ -3148,6 +3174,7 @@ def _compute_objective_value(result: dict, x: np.ndarray, requirements: dict, co A_throat_check=A_throat_check, ) + _infeas_trace(_tr, "geometry", infeasibility_score) # --- Evaluation Results --- eval_success = result.get('success', False) if isinstance(result, dict) else False # Runner.evaluate typically omits success; infer from finite thrust/Pc when absent @@ -3259,6 +3286,7 @@ def _compute_objective_value(result: dict, x: np.ndarray, requirements: dict, co infeasibility_score += max(0.0, (effective_margin - chugging_margin) / effective_margin) ** 2 infeasibility_score += max(0.0, (effective_margin - acoustic_margin) / effective_margin) ** 2 infeasibility_score += max(0.0, (effective_margin - feed_margin) / effective_margin) ** 2 + _infeas_trace(_tr, "stability", infeasibility_score) # Regularization: Cf band def _hinge_band(val, lo, hi, scale=1.0): @@ -3436,6 +3464,8 @@ def _hinge_band(val, lo, hi, scale=1.0): max_outward_deg=float(constants.get("layer1_resultant_tilt_max_deg", 0.0)), scale_deg=float(constants.get("layer1_resultant_tilt_scale_deg", 2.0)), ) + _infeas_trace(_tr, "wall_tilt", infeasibility_score) + _infeas_trace(_tr, "tilt_deg", _tilt) momentum_term = _impinging_momentum_asymmetric_squared( R_val, wall_side_multiplier=float( @@ -3974,6 +4004,60 @@ def _layer1_emit_objective_plot_point( pass +def _layer1_infeasibility_reason(runner, x, requirements: dict, constants: dict) -> Optional[str]: + """One sentence on WHICH hard constraint held the best candidate out of the feasible set. + + A run whose every candidate sat on the 1e6 infeasibility floor used to end with + "objective inf" and "Validation failed", which told the user nothing. Re-evaluate the best + design in-process with the infeasibility trace on and name the dominant gate. + """ + global _INFEAS_TRACE_ON + try: + idx_P_O = 11 if constants.get("injector_type") == "impinging" else 8 + P_O = float(x[idx_P_O]) * 6894.76 + P_F = float(x[idx_P_O + 1]) * 6894.76 + result = runner.evaluate(P_O, P_F, silent=True) + except Exception as e: + return f"No feasible design: the best candidate could not even be evaluated ({type(e).__name__}: {str(e)[:120]})." + prev, n0 = _INFEAS_TRACE_ON, len(_INFEAS_TRACE) + _INFEAS_TRACE_ON = True + try: + _compute_objective_value(result, np.asarray(x, dtype=float), requirements, constants) + except Exception as e: + return f"No feasible design: objective re-evaluation failed ({type(e).__name__})." + finally: + _INFEAS_TRACE_ON = prev + if len(_INFEAS_TRACE) <= n0: + return None + tr = _INFEAS_TRACE.pop() + geom = float(tr.get("geometry", 0.0)) + stab = float(tr.get("stability", geom)) - geom + wall = float(tr.get("wall_tilt", tr.get("stability", geom))) - float(tr.get("stability", geom)) + parts = [] + if wall > 0: + tilt = tr.get("tilt_deg", float("nan")) + lim = float(constants.get("layer1_resultant_tilt_max_deg", 0.0)) + parts.append((wall, f"the spray resultant tilts {tilt:+.1f} deg outward toward the wall (limit {lim:g} deg) -- " + "lower the oxidizer jet angle, raise the fuel jet angle, or widen the angle bands")) + if stab > 0: + parts.append((stab, "a stability gate (minimum score / margins) is not met -- lower min_stability_score or stiffen the injector")) + if geom > 0: + parts.append((geom, "injector packing, flow capacity, or chamber proportions violate a hard limit -- widen the chamber OD or the jet bounds")) + if not parts: + # Re-evaluated in isolation the best design passes every gate: the search is stuck ON a + # constraint boundary and CMA's samples keep landing a hair outside it. For a doublet that + # is almost always the spray-resultant tilt limit, which defaults to exactly 0 deg outward. + lim = float(constants.get("layer1_resultant_tilt_max_deg", 0.0)) + if constants.get("injector_type") == "impinging": + return ("No candidate cleared every hard constraint, yet the best design re-evaluates as feasible " + f"on its own: the search is pinned on the spray-resultant tilt limit ({lim:g} deg outward). " + "Re-run, or allow a degree of outward tilt (layer1_resultant_tilt_max_deg) if the liner can take it.") + return ("No candidate cleared every hard constraint, yet the best design re-evaluates as feasible on its " + "own: the search is pinned on a constraint boundary. Re-run, or relax the tightest gate slightly.") + parts.sort(key=lambda t: -t[0]) + return "No candidate cleared every hard constraint. On the best one, " + parts[0][1] + "." + + def run_layer1_global_search( objective: Callable[[np.ndarray], float], bounds: list, @@ -4084,6 +4168,29 @@ def wrapped_obj(v: np.ndarray) -> float: return best_x +def _layer1_check_of_target_in_cea_range(config_obj: Any, optimal_of: Any) -> None: + """Refuse a target O/F the propellant's CEA table cannot evaluate. + + A propellant switch keeps the previous design target, so an ethalox target of 1.4 left + behind on a methalox config used to run a full optimization against a table that stops at + 2.4 -- the mixture ratio pinned at the table edge and the run returned a huge objective with + nothing to say why. Fail before any work is done, with the fix in the message. + """ + try: + mr_range = config_obj.combustion.cea.MR_range + lo, hi = float(mr_range[0]), float(mr_range[1]) + of = float(optimal_of) + except (AttributeError, TypeError, ValueError, IndexError): + return + if not (np.isfinite(of) and lo <= of <= hi): + preset = getattr(config_obj, "propellant_preset", None) or "this propellant" + raise ValueError( + f"Design target O/F {of:.2f} is outside the CEA table for {preset} " + f"(MR_range [{lo:.2f}, {hi:.2f}]). Set optimal_of_ratio inside that range in Design " + f"Requirements -- switching propellant keeps the previous target." + ) + + def run_layer1_optimization( config_obj: PintleEngineConfig, runner: PintleEngineRunner, @@ -4184,6 +4291,7 @@ def check_stop(): # Extract requirements target_thrust = requirements.get("target_thrust", 7000.0) optimal_of = requirements.get("optimal_of_ratio", 2.3) + _layer1_check_of_target_in_cea_range(config_obj, optimal_of) min_stability = float(requirements.get("min_stability_margin", _LAYER1_DEFAULT_MIN_STABILITY_MARGIN)) def _resolve_Lstar_bounds_from_req_and_config() -> Tuple[float, float]: @@ -7806,6 +7914,16 @@ def _validation_evaluate_or_bundle( optimized_config_runner.graphite_insert.enabled = False optimized_runner = PintleEngineRunner(optimized_config_runner) + + # When nothing was feasible, say which gate held the best candidate out (see the helper). + infeasible_reason = None + try: + if best_x is not None and not _layer1_feasible_scalar_objective(float(opt_state.get("best_objective", float("inf")))): + infeasible_reason = _layer1_infeasibility_reason(optimized_runner, best_x, requirements, constants_dict) + if infeasible_reason and log_status: + log_status("warning", infeasible_reason) + except Exception: + infeasible_reason = None # Use stored validation results if available if "best_results_for_validation" in opt_state and opt_state["best_results_for_validation"] is not None: @@ -8511,6 +8629,7 @@ def _as_finite_float_or_nan(v: Any) -> float: else {} ), "primary_relative_residual": _prim_rel, + "infeasible_reason": infeasible_reason, }, "exit_pressure_targeting": { "target_P_exit": target_P_exit, # Atmospheric pressure from environment config (GPS/GFS-derived) diff --git a/EngineDesign/engine/pipeline/config_schemas.py b/EngineDesign/engine/pipeline/config_schemas.py index 5d69e8440..979125e31 100644 --- a/EngineDesign/engine/pipeline/config_schemas.py +++ b/EngineDesign/engine/pipeline/config_schemas.py @@ -157,6 +157,15 @@ class FeedSystemConfig(BaseModel): default="none", description="Pressure function type" ) + length: Optional[float] = Field( + default=None, + gt=0, + description=( + "Feed-line length from tank outlet to injector manifold [m]. Sets the line inertance " + "(length / area) in the chug model. Leave unset and the stability model records a " + "0.305 m assumption instead of using it silently." + ), + ) @model_validator(mode="before") @classmethod @@ -246,7 +255,6 @@ class RegenCoolingConfig(BaseModel): n_segments: int = Field(default=20, gt=0, description="Number of axial segments for heat-transfer integration") gas_turbulence_intensity: float = Field(default=0.1, ge=0, description="Estimated turbulence intensity of hot gas (0-1)") coolant_turbulence_intensity: float = Field(default=0.05, ge=0, description="Estimated turbulence intensity of coolant (0-1)") - hot_gas_cp: float = Field(default=2200.0, gt=0, description="Hot-gas specific heat [J/(kg·K)]") recovery_factor: Optional[float] = Field(default=None, gt=0, le=1, description="Turbulent boundary layer recovery factor for adiabatic wall temperature (Taw = Tc × recovery_factor). Typical range: 0.90-0.98. If None, uses default from constants.") @@ -359,8 +367,6 @@ class StainlessSteelCaseConfig(BaseModel): specific_heat: float = Field(default=500.0, gt=0, description="Specific heat [J/(kg·K)]") max_temperature: float = Field(default=1000.0, gt=0, description="Maximum allowable temperature [K] (melting point ~1700K, but limit lower for structural integrity)") emissivity: float = Field(default=0.3, ge=0, le=1, description="Surface emissivity") - yield_strength: float = Field(default=200e6, gt=0, description="Yield strength at max temp [Pa]") - youngs_modulus: float = Field(default=200e9, gt=0, description="Young's modulus [Pa]") class AblativeCoolingConfig(BaseModel): @@ -641,10 +647,6 @@ class CombustionEfficiencyConfig(BaseModel): # --- Mixing efficiency (Rupe momentum-ratio model) --- # Replaces the old k-e near-field mixing model + the eta_turbulence step-function. # eta_mix = Em_peak * exp(-(ln(R/R_opt))^2 / (2*sigma^2)), R = injector momentum ratio. - mixing_model: Literal["rupe"] = Field( - default="rupe", - description="Mixing efficiency model. 'rupe': momentum-ratio mixing efficiency (Rupe/SP-8089)." - ) Em_peak: float = Field( default=0.96, ge=0.5, le=1.0, description="Peak (best-achievable) mixing efficiency at the balanced momentum ratio. " @@ -806,6 +808,57 @@ class CombustionConfig(BaseModel): efficiency: CombustionEfficiencyConfig = Field(default_factory=CombustionEfficiencyConfig) +class StabilityConfig(BaseModel): + """Inputs to the combustion / feed stability model that belong to neither the propellant + (``fluids``) nor the plumbing (``feed_system``): the combustion-response calibration, the + nozzle-entrance Mach the acoustic damping uses, the acoustic damping coefficients, and the + dome-regulator dynamics. A ``None`` here means "derive it" and the derivation is recorded in + the assumptions registry (rich report -> assumptions.fallbacks_used), never substituted silently. + """ + n_interaction: float = Field( + default=0.5, gt=0, + description="Crocco interaction index n (calibration range 0.3-0.6). The forward-mode slider overrides it per run.", + ) + chi_acoustic: float = Field( + default=0.15, gt=0, le=1, + description="Sensitive-lag fraction chi: tau_sens = chi * tau_vap for the acoustic n-tau driving.", + ) + mach_nozzle_entrance: Optional[float] = Field( + default=None, gt=0, lt=1, + description="Mean Mach at the nozzle entrance (sets nozzle damping). None = solve it from the contraction ratio (isentropic, subsonic).", + ) + damping_injector_frac: float = Field( + default=0.02, ge=0, + description="Injector-face acoustic damping as a fraction of pi*f [-]. First-cut; calibrate against a cold ring-down test.", + ) + damping_twophase_frac: float = Field( + default=0.03, ge=0, + description="Two-phase (droplet) acoustic damping as a fraction of pi*f*droplet_loading [-]. First-cut.", + ) + droplet_loading: float = Field( + default=1.0, ge=0, + description="Relative liquid loading near the injector face for the two-phase damping term [-].", + ) + acoustic_gate_alpha_offset: float = Field( + default=350.0, ge=0, + description=( + "Calibration allowance for the acoustic gate [1/s]: a mode growing slower than this still " + "maps to a neutral gate margin because the a-priori damping coefficients are un-measured. " + "Set 0 for the strict alpha < 0 criterion." + ), + ) + regulator_enabled: bool = Field(default=True, description="Model the dome regulator upstream of each tank in the chug loop.") + regulator_corner_hz: float = Field(default=3.0, gt=0, description="Regulator response corner frequency [Hz].") + regulator_Z_hf: float = Field( + default=0.0, ge=0, + description="Regulator high-frequency series impedance [Pa*s/kg]. 0 = ideal pressure source (optimistic); measure via a step test.", + ) + regulator_max_excursion_psi: float = Field( + default=0.0, ge=0, + description="Regulator outlet pressure excursion bound [psi]. Reporting only; not a pole-shifter.", + ) + + class ChamberGeometryConfig(BaseModel): """ Unified chamber geometry configuration for solve_chamber_geometry_with_cea. @@ -908,13 +961,13 @@ class LOXTankConfig(BaseModel): class FuelTankConfig(BaseModel): - """Fuel tank geometry configuration for flight simulation""" - rp1_h: float = Field(gt=0, description="RP-1 tank height (internal cylindrical length, not including end caps) [m]") - rp1_radius: float = Field(gt=0, description="RP-1 tank internal radius [m]") + """Fuel tank geometry configuration for flight simulation. The rp1_* field names are legacy; the tank holds whichever fuel the config names.""" + rp1_h: float = Field(gt=0, description="Fuel tank height (internal cylindrical length, not including end caps) [m]") + rp1_radius: float = Field(gt=0, description="Fuel tank internal radius [m]") fuel_tank_pos: float = Field(description="Fuel tank center position relative to nozzle exit (positive = above, negative = below nozzle) [m]") - mass: Optional[float] = Field(default=None, gt=0, description="Initial RP-1 PROPELLANT mass [kg] (liquid only, not tank structure). Depletes during burn.") + mass: Optional[float] = Field(default=None, gt=0, description="Initial fuel PROPELLANT mass [kg] (liquid only, not tank structure). Depletes during burn.") initial_pressure_psi: Optional[float] = Field(default=None, gt=0, description="Initial fuel tank pressure [psi]") - tank_volume_m3: Optional[float] = Field(default=None, gt=0, description="RP-1 tank volume [m³]. If not provided, will be calculated from rp1_h and rp1_radius using π×r²×h") + tank_volume_m3: Optional[float] = Field(default=None, gt=0, description="Fuel tank volume [m³]. If not provided, will be calculated from rp1_h and rp1_radius using π×r²×h (field names are legacy)") class PressTankConfig(BaseModel): @@ -1915,6 +1968,7 @@ class PintleEngineConfig(BaseModel): chamber: Optional[ChamberConfig] = Field(default=None, description="Legacy chamber config (use chamber_geometry instead)") nozzle: Optional[NozzleConfig] = Field(default=None, description="Legacy nozzle config (use chamber_geometry instead)") solver: SolverConfig = Field(default_factory=SolverConfig) + stability: StabilityConfig = Field(default_factory=StabilityConfig, description="Combustion / feed stability model inputs (calibration, regulator, acoustic damping)") optimizer: Optional[OptimizerConfig] = Field(default=None, description="Optimizer configuration") # Flight simulation fields (optional) lox_tank: Optional[LOXTankConfig] = Field(default=None, description="LOX tank configuration for flight simulation") @@ -1992,11 +2046,7 @@ def sync_burn_time_fields(self): sync_burn_time_fields(self) return self - class Config: - # NOTE: "allow" ACCEPTS unknown YAML keys (stores them as extra attributes) - # — it does not reject them. Kept permissive for legacy configs; typo'd - # keys are therefore silently inert. - extra = "allow" + model_config = ConfigDict(extra="allow") def ensure_chamber_geometry(config: PintleEngineConfig) -> ChamberGeometryConfig: diff --git a/EngineDesign/engine/pipeline/physics_based_replacements.py b/EngineDesign/engine/pipeline/physics_based_replacements.py index 4ae5d8fcb..d35c6d7ed 100644 --- a/EngineDesign/engine/pipeline/physics_based_replacements.py +++ b/EngineDesign/engine/pipeline/physics_based_replacements.py @@ -239,62 +239,6 @@ def calculate_throat_heat_flux_physics( return float(heat_flux_throat) -def calculate_recirculation_intensity_physics( - fuel_velocity: float, - lox_velocity: float, - d_pintle_tip: float, - D_chamber: float, - Re_injector: float, -) -> float: - """ - Calculate recirculation intensity based on physics. - - Physics: - - Recirculation intensity depends on velocity ratio - - Higher velocity difference → stronger recirculation - - Depends on Reynolds number (turbulent flow) - - Scales with injector size - - Parameters: - ----------- - fuel_velocity : float - Fuel injection velocity [m/s] - lox_velocity : float - LOX injection velocity [m/s] - d_pintle_tip : float - Pintle tip diameter [m] - D_chamber : float - Chamber diameter [m] - Re_injector : float - Injector Reynolds number - - Returns: - -------- - intensity : float - Recirculation intensity (0-1) - """ - # Velocity difference drives recirculation - velocity_diff = abs(fuel_velocity - lox_velocity) - velocity_avg = (fuel_velocity + lox_velocity) / 2.0 - velocity_ratio = velocity_diff / (velocity_avg + 1e-10) - - # Base intensity from velocity ratio - # Higher velocity difference → stronger recirculation - base_intensity = 0.2 * velocity_ratio # Physics-based scaling - - # Reynolds number effect: higher Re → more turbulent → stronger recirculation - Re_factor = np.clip(Re_injector / 1e4, 0.5, 2.0) - Re_enhancement = 1.0 + 0.3 * np.log10(max(Re_factor, 0.1)) - - # Pintle size effect: larger pintle → larger recirculation - pintle_ratio = d_pintle_tip / (D_chamber + 1e-10) - pintle_factor = 1.0 + 0.2 * np.clip(pintle_ratio - 0.1, 0.0, 0.3) - - intensity = base_intensity * Re_enhancement * pintle_factor - - return float(np.clip(intensity, 0.0, 0.8)) - - def calculate_turbulence_enhancement_physics( Re_throat: float, velocity_ratio: float, diff --git a/EngineDesign/engine/pipeline/reaction_chemistry.py b/EngineDesign/engine/pipeline/reaction_chemistry.py index b971de173..fd2837b13 100644 --- a/EngineDesign/engine/pipeline/reaction_chemistry.py +++ b/EngineDesign/engine/pipeline/reaction_chemistry.py @@ -24,14 +24,8 @@ def _fuel_name_for_kinetics(config: PintleEngineConfig) -> str: """Label used for Arrhenius fuel-type branching (matches reaction_rate_constant lists). - Priority: legacy ``propellants.fuel.name`` → ``combustion.cea.fuel_name`` → ``fluids['fuel'].name`` - → ``\"RP-1\"``. + Priority: ``combustion.cea.fuel_name`` -> ``fluids['fuel'].name`` -> generic hydrocarbon ("RP-1"). """ - prop = getattr(config, "propellants", None) - if prop is not None and getattr(prop, "fuel", None) is not None: - name = getattr(prop.fuel, "name", None) - if name: - return str(name) cea = getattr(getattr(config, "combustion", None), "cea", None) if cea is not None: fn = getattr(cea, "fuel_name", None) @@ -48,21 +42,18 @@ def _fuel_name_for_kinetics(config: PintleEngineConfig) -> str: def _fuel_props_for_evaporation(config: PintleEngineConfig) -> Optional[Dict[str, float]]: """Density and boiling point for droplet evaporation time scale.""" - prop = getattr(config, "propellants", None) - if prop is not None and getattr(prop, "fuel", None) is not None: - bp = getattr(prop.fuel, "boiling_point", None) - return { - "density": float(prop.fuel.density), - "boiling_point": float(bp) if bp is not None else 489.0, - } fluids = getattr(config, "fluids", None) if not fluids or "fuel" not in fluids: return None fuel = fluids["fuel"] bp = getattr(fuel, "boiling_point", None) + if bp is None: + from engine.pipeline.assumptions import assume + bp = assume("kinetics.fuel.boiling_point", 489.0, unit="K", + reason=f"fluids.fuel.boiling_point missing for {getattr(fuel, 'name', '?')} (RP-1 value used; set it or load a propellant preset)") return { "density": float(fuel.density), - "boiling_point": float(bp) if bp is not None else 489.0, + "boiling_point": float(bp), } @@ -142,6 +133,20 @@ def calculate_reaction_progress( return float(progress) +_EA_MR_W = 0.10 # half-width [MR] of the smooth blend at the 1.5 / 3.0 boundaries + + +def _ea_mr_factor(MR: float) -> float: + """Activation-energy multiplier vs O/F: 1.2 (fuel-rich) -> 1.0 -> 0.9 (oxidizer-rich), C1-smooth.""" + def smooth(t): + t = min(max(t, 0.0), 1.0) + return t * t * (3.0 - 2.0 * t) + w = _EA_MR_W + lo = 1.2 + (1.0 - 1.2) * smooth((MR - (1.5 - w)) / (2.0 * w)) # 1.2 -> 1.0 around 1.5 + hi = 1.0 + (0.9 - 1.0) * smooth((MR - (3.0 - w)) / (2.0 * w)) # 1.0 -> 0.9 around 3.0 + return float(lo if MR < 2.25 else hi) + + def calculate_reaction_rate_constant( Pc: float, Tc: float, @@ -251,12 +256,12 @@ def calculate_reaction_rate_constant( Ea = 80000.0 n_pressure = 0.8 - # Adjust activation energy based on mixture ratio - # Fuel-rich or oxidizer-rich can have different effective Ea - if MR < 1.5: # Fuel-rich: more complex chemistry - Ea *= 1.2 # Higher effective activation energy - elif MR > 3.0: # Oxidizer-rich: simpler chemistry - Ea *= 0.9 # Lower effective activation energy + # Effective activation energy vs mixture ratio: fuel-rich chemistry is slower (x1.2), + # oxidizer-rich faster (x0.9). Blended with a C1 smoothstep across +/-0.10 MR of the 1.5 and + # 3.0 boundaries -- the previous hard `if` put a 20% step in Ea at exactly MR = 1.5, which + # is inside the ethalox design band, and a step in the objective breaks any gradient-based + # or secant refinement that crosses it (same defect as combustion_physics._ea_norm_from_mr). + Ea *= _ea_mr_factor(MR) # Pre-exponential with pressure dependence # A(P) = A0 × (P / P0)^n_pre, where P0 = 1 MPa reference diff --git a/EngineDesign/engine/pipeline/stability/analysis.py b/EngineDesign/engine/pipeline/stability/analysis.py index 4fde0bc26..92a336594 100644 --- a/EngineDesign/engine/pipeline/stability/analysis.py +++ b/EngineDesign/engine/pipeline/stability/analysis.py @@ -19,7 +19,8 @@ import os from typing import Dict, Tuple, Optional, List, Any import numpy as np -from engine.pipeline.config_schemas import PintleEngineConfig +from engine.pipeline.config_schemas import PintleEngineConfig, StabilityConfig +from engine.pipeline.constants import DEFAULT_HOT_GAS_THERMAL_COND_W_M_K, DEFAULT_HOT_GAS_VISC_PA_S # --------------------------------------------------------------------------- @@ -36,124 +37,48 @@ def calculate_chugging_frequency( Tc: Optional[float] = None, ) -> Dict[str, float]: """ - Estimate low frequency combustion instability (chugging) characteristics. + Order-of-magnitude chug (bulk-mode) frequency estimates from chamber geometry alone. - We combine two simple notions: - - Residence time, tau_res = L* / c* - - Helmholtz-like volume-compliance mode if gas properties are known - - Parameters - ---------- - chamber_volume : float - Chamber volume [m^3] - throat_area : float - Throat area [m^2] - cstar : float - Characteristic velocity [m/s] - gamma : float - Specific heat ratio [-] - Pc : float - Chamber pressure [Pa] - R : float, optional - Gas constant [J/(kg K)]. If provided with Tc, used for Helmholtz estimate. - Tc : float, optional - Chamber temperature [K]. If provided with R, used for Helmholtz estimate. + Two estimates: the residence-time frequency ``1 / (2 pi tau_res)`` with ``tau_res = L*/c*``, and + a Helmholtz bulk mode using the throat as the neck. These are *placeholders* -- the physical chug + frequency comes from the feed-coupled loop in ``chug.py`` and overwrites ``frequency`` in + ``comprehensive_stability_analysis``. The old heuristic "stability_index" / "stability_margin" + that used to ride along here (floored at 0.4 and mapped to a margin so that "reasonable designs + achieve required margins") was not a physical quantity and has been removed; margins come from + the gain-margin model only. Returns ------- dict - - frequency: dominant chugging frequency [Hz] - - frequency_residence: frequency from 1 / (2 pi tau_res) [Hz] - - frequency_helmholtz: Helmholtz estimate if possible [Hz or np.nan] - - period: oscillation period [s] - - stability_index: heuristic index (higher is better) - - stability_margin: backward compatibility field (maps from stability_index) - - tau_residence: residence time L* / c* [s] + - frequency: Helmholtz estimate when gas properties are known, else the residence estimate [Hz] + - frequency_residence, frequency_helmholtz: the two estimates [Hz] (nan if unavailable) + - period: 1 / frequency [s] + - tau_residence: L* / c* [s] - Lstar: characteristic length [m] """ if throat_area <= 0.0 or chamber_volume <= 0.0 or cstar <= 0.0: - # Fallback values - Lstar = 1.0 - tau_residence = 1.0e-3 - else: - Lstar = chamber_volume / throat_area - tau_residence = Lstar / cstar - - # Frequency from residence time + return {"frequency": float("nan"), "frequency_residence": float("nan"), + "frequency_helmholtz": float("nan"), "period": float("nan"), + "tau_residence": float("nan"), "Lstar": float("nan")} + Lstar = chamber_volume / throat_area + tau_residence = Lstar / cstar freq_res = 1.0 / (2.0 * np.pi * tau_residence) - # Helmholtz-like frequency if we know gas properties - # f_H ≈ (c / (2 pi)) * sqrt(A_neck / (V * L_eff)) - # Use throat as neck and L_eff ~ D_throat - if R is not None and Tc is not None and throat_area > 0.0 and chamber_volume > 0.0: - # FIXED: Add safety checks for sqrt operations - a = float(np.sqrt(max(0, gamma * R * Tc))) - d_throat = np.sqrt(max(0, 4.0 * throat_area / np.pi)) + # Helmholtz-like bulk mode: f_H = (a / 2pi) * sqrt(A_neck / (V * L_eff)), neck = throat, + # L_eff ~ half a throat diameter. + freq_helm = float("nan") + if R is not None and Tc is not None and gamma * R * Tc > 0: + a = float(np.sqrt(gamma * R * Tc)) + d_throat = np.sqrt(4.0 * throat_area / np.pi) L_eff = max(0.5 * d_throat, 1.0e-3) - sqrt_arg = throat_area / (chamber_volume * L_eff) if chamber_volume * L_eff > 0 else 0.0 - freq_helm = (a / (2.0 * np.pi)) * np.sqrt(max(0, sqrt_arg)) - else: - freq_helm = np.nan - - # Choose dominant frequency for chugging - if np.isfinite(freq_helm): - freq = 0.5 * freq_res + 0.5 * freq_helm - else: - freq = freq_res - - # Clamp to an engineering range for reporting - freq = float(np.clip(freq, 1.0, 2000.0)) - freq_res = float(freq_res) - freq_helm = float(freq_helm) if np.isfinite(freq_helm) else float("nan") - - period = 1.0 / freq if freq > 0.0 else float("inf") - - # Simple stability index: - # - Better if Pc is higher - # - Better if L* is reasonably large (say ≥ 0.8 m) - # - Penalize if chugging frequency is very low (hard to damp) or in a problematic band - # FIXED: More lenient factors to allow reasonable designs to achieve stable margins - Pc_ref = 1.0e6 - Lstar_ref = 1.0 - # More lenient Pc factor - even 0.5 MPa can be acceptable - Pc_factor = min(1.0, (Pc / Pc_ref) ** 0.3) if Pc > 0 else 0.0 - # More lenient Lstar factor - even 0.6 m can be acceptable - Lstar_factor = min(1.0, (Lstar / Lstar_ref) ** 0.2) if Lstar > 0 else 0.0 - # Ensure minimum factors for reasonable designs - Pc_factor = max(Pc_factor, 0.6) if Pc > 0.3e6 else Pc_factor # At least 0.6 for Pc > 0.3 MPa - Lstar_factor = max(Lstar_factor, 0.7) if Lstar > 0.6 else Lstar_factor # At least 0.7 for L* > 0.6 m - - # Frequency health factor: prefer 20 to 400 Hz for chugging - # Much more lenient penalties to allow optimizer to find feasible solutions - if freq < 5.0: - f_factor = 0.6 # Very low frequencies - still penalized but not as harsh - elif freq < 10.0: - f_factor = 0.75 # Low frequencies - moderate penalty - elif freq > 600.0: - f_factor = 0.9 # High frequencies - minimal penalty - elif freq > 400.0: - f_factor = 0.95 # Moderate-high frequencies - very small penalty - else: - f_factor = 1.0 # Ideal range - - stability_index = Pc_factor * Lstar_factor * f_factor - # Ensure minimum index for reasonable designs - stability_index = max(stability_index, 0.4) # Minimum 0.4 for any reasonable design - - # Backward compatibility: map stability_index to stability_margin - # FIXED: More generous mapping to ensure reasonable designs can meet requirements - # For a reasonable design (index ~ 0.6-0.8), we want margin ~ 1.2-1.5 - # New mapping: margin = stability_index * 1.5 + 0.4 (gives 1.3 for index=0.6, 1.6 for index=0.8, 1.9 for index=1.0) - # This ensures reasonable designs can achieve required margins - stability_margin = stability_index * 1.5 + 0.4 # More generous mapping + freq_helm = float((a / (2.0 * np.pi)) * np.sqrt(throat_area / (chamber_volume * L_eff))) + freq = freq_helm if np.isfinite(freq_helm) else float(freq_res) return { "frequency": float(freq), - "frequency_residence": freq_res, - "frequency_helmholtz": freq_helm, - "period": float(period), - "stability_index": float(stability_index), - "stability_margin": float(stability_margin), # Backward compatibility + "frequency_residence": float(freq_res), + "frequency_helmholtz": float(freq_helm), + "period": float(1.0 / freq) if freq > 0 else float("inf"), "tau_residence": float(tau_residence), "Lstar": float(Lstar), } @@ -231,90 +156,54 @@ def analyze_feed_system_stability( pressure_drop: float, ) -> Dict[str, float]: """ - Analyze feed system stability (POGO, surge, water hammer). + Feed-line acoustics and the water-hammer bound for one propellant line. Parameters ---------- feed_line_length : float Feed line length [m] feed_line_diameter : float - Feed line diameter [m] + Feed line bore [m] propellant_density : float - Propellant density [kg/m^3] + Liquid density [kg/m^3] bulk_modulus : float - Bulk modulus [Pa] + Liquid bulk modulus [Pa] flow_velocity : float - Flow velocity [m/s] + Mean line velocity [m/s] pressure_drop : float - Pressure drop across feed system [Pa] + Tank-to-chamber pressure drop [Pa] Returns ------- dict - - pogo_frequency: quarter wave frequency [Hz] - - surge_frequency: half wave frequency [Hz] - - water_hammer_pressure: spike for full stop [Pa] + - pogo_frequency: quarter-wave line mode (closed-open) [Hz] + - surge_frequency: half-wave line mode (closed-closed) [Hz] + - water_hammer_pressure: Joukowsky spike for an instantaneous stop, rho*a*dv [Pa] - water_hammer_margin: pressure_drop / spike [-] - - stability_margin: backward compatibility field (maps from water_hammer_margin) - - sound_speed: wave speed in propellant [m/s] + - sound_speed: wave speed in the liquid [m/s] + + The feed-coupled *stability* margin is the chug gain margin from ``chug.py``; the caller writes + it into this dict as ``stability_margin``. The piecewise water-hammer-to-margin mapping that used + to live here (tuned so "typical optimized designs meet the 1.20 requirement", including a branch + that was literally a constant) was not a stability criterion and has been removed. """ L = max(feed_line_length, 1.0e-3) rho = propellant_density K = bulk_modulus sound_speed = float(np.sqrt(K / rho)) - pogo_frequency = float(sound_speed / (4.0 * L)) # closed-open surge_frequency = float(sound_speed / (2.0 * L)) # closed-closed delta_v = max(flow_velocity, 0.0) water_hammer_pressure = float(rho * sound_speed * delta_v) - - if water_hammer_pressure > 0.0: - water_hammer_margin = float(pressure_drop / water_hammer_pressure) - else: - water_hammer_margin = float("inf") - - # FIXED: Map water_hammer_margin to stability_margin accounting for real-world factors - # The theoretical water_hammer_pressure assumes instantaneous stop, which is overly conservative. - # In reality: - # - Valves close over time (0.1-1.0 s), reducing actual pressure spike by 50-90% - # - Systems have accumulators, surge suppressors, and compliance - # - Actual water hammer is typically 10-50% of theoretical maximum - # - # Map water_hammer_margin to stability_margin: - # - Display requirement is >= 1.20 (full min_stability_margin) - # - Optimizer convergence uses >= 0.96 (80% of 1.2) - # - Adjusted mapping to ensure typical designs meet the full 1.20 requirement - # - # Use a scaling function that makes reasonable designs achievable: - # For water_hammer_margin = 0.15-0.2 (typical), we want stability_margin >= 1.20 - if water_hammer_margin >= 0.5: - # Good margin: scale linearly from 0.5 -> 1.2 to higher values - stability_margin = 1.2 + (water_hammer_margin - 0.5) * 1.0 # 0.5 -> 1.2, 1.0 -> 1.7 - elif water_hammer_margin >= 0.05: - # Moderate margin: scale from 0.05 -> 1.20 to 0.5 -> 1.2 - # Typical optimized designs (0.05-0.2) should meet the 1.20 requirement - # This accounts for real-world valve closure times and system compliance - stability_margin = 1.20 + (water_hammer_margin - 0.05) / 0.45 * 0.0 # 0.05 -> 1.20, 0.5 -> 1.20 - elif water_hammer_margin >= 0.03: - # Very low margin: scale from 0.03 -> 1.15 to 0.05 -> 1.20 - # Still acceptable with proper engineering (valve closure, accumulators) - stability_margin = 1.15 + (water_hammer_margin - 0.03) / 0.02 * 0.05 # 0.03 -> 1.15, 0.05 -> 1.20 - else: - # Extremely low margin: scale from 0.0 -> 1.00 to 0.03 -> 1.15 - # Still give reasonable margin since real systems have mitigations - stability_margin = 1.00 + water_hammer_margin / 0.03 * 0.15 # 0.0 -> 1.00, 0.03 -> 1.15 - - # Clamp to reasonable range - stability_margin = float(np.clip(stability_margin, 0.1, 5.0)) + water_hammer_margin = float(pressure_drop / water_hammer_pressure) if water_hammer_pressure > 0.0 else float("inf") return { "pogo_frequency": pogo_frequency, "surge_frequency": surge_frequency, "water_hammer_pressure": water_hammer_pressure, "water_hammer_margin": water_hammer_margin, - "stability_margin": stability_margin, # Backward compatibility "sound_speed": sound_speed, } @@ -322,18 +211,18 @@ def analyze_feed_system_stability( # --------------------------------------------------------------------------- # Physical stability margins (new model) — fast tiers for the per-eval path # --------------------------------------------------------------------------- -# Interim gate-margin mappings: PHYSICAL and monotone in the growth rate, but conservatively centered so -# currently-healthy designs pass (directive: keep the gate's impact ~as-is) while clearly-unstable designs -# fail. The absolute calibration of the chug feed/regulator params and the acoustic damping coefficients is -# un-measured; these constants are re-tuned once tests T5/T6/T7/H3 land. Documented in the rebuild plan §5. -_CHUG_GATE_CENTER = 0.80 # chug gain margin -> "neutral" (gate margin 1.0) +# Gate-margin mappings, monotone in the physical quantity and centred on the physical criterion: +# * chug: Nyquist gain margin GM. GM = 1 is the stability boundary, so it maps to a neutral gate +# margin of 1.0 (the old centre of 0.80 called a GM of 0.85 -- an unstable loop -- "stable"). +# * acoustic: net growth rate alpha of the worst mode. alpha = 0 is the boundary, but the a-priori +# damping coefficients are un-measured, so a configurable allowance +# (stability.acoustic_gate_alpha_offset, default 350 1/s) keeps the gate's calibration explicit +# rather than hidden. Set it to 0 for the strict criterion. +_CHUG_GATE_CENTER = 1.00 # chug gain margin at the stability boundary -> gate margin 1.0 _CHUG_GATE_SCALE = 0.20 _GATE_SPAN = 0.30 # gate margin ranges ~[0.7, 1.3] -_ACOUSTIC_GATE_OFFSET = 350.0 # [1/s] alpha offset so marginal acoustic still passes +_ACOUSTIC_GATE_OFFSET = 350.0 # [1/s] default allowance; overridden by StabilityConfig _ACOUSTIC_GATE_SCALE = 1000.0 # [1/s] -# LOX property fallbacks (default.yaml leaves LOX latent_heat / boiling_point null). -_LOX_HFG_DEFAULT = 213000.0 # J/kg -_LOX_TBOIL_DEFAULT = 90.2 # K def _chug_gate_margin(gain_margin: float) -> float: @@ -342,10 +231,10 @@ def _chug_gate_margin(gain_margin: float) -> float: return float(1.0 + _GATE_SPAN * np.tanh((gain_margin - _CHUG_GATE_CENTER) / _CHUG_GATE_SCALE)) -def _acoustic_gate_margin(alpha_max: float) -> float: +def _acoustic_gate_margin(alpha_max: float, alpha_offset: float = _ACOUSTIC_GATE_OFFSET) -> float: if not np.isfinite(alpha_max): return 1.10 - return float(1.0 + _GATE_SPAN * np.tanh((_ACOUSTIC_GATE_OFFSET - alpha_max) / _ACOUSTIC_GATE_SCALE)) + return float(1.0 + _GATE_SPAN * np.tanh((float(alpha_offset) - alpha_max) / _ACOUSTIC_GATE_SCALE)) def _fluid_attr(fluids, key, attr, default): @@ -372,9 +261,91 @@ def _feed_attr(config, key, attr, default): return float(default) -# Default combustion-response parameters (calibration targets, see [Phys §2, §5.3]). -_N_INTERACTION_DEFAULT = 0.5 # interaction index n (sweep 0.3-0.6) -_CHI_ACOUSTIC_DEFAULT = 0.15 # sensitive-fraction chi for acoustic (tau_sens = chi*tau_vap) +# Handbook thermodynamic fallbacks BY FLUID, used only when the config omits a property. Every use +# is recorded in the assumptions registry. Previously the fuel fallbacks were methane's (h_fg 510 kJ/kg, +# T_boil 111.6 K) regardless of which fuel the config named, and the oxidizer's were LOX's. +# density kg/m^3 latent heat J/kg boiling point K (1 atm) +_FLUID_THERMO_FALLBACKS = { + "lox": (1140.0, 213000.0, 90.2), + "methane": ( 422.6, 510000.0, 111.65), + "ethanol": ( 789.0, 838000.0, 351.4), + "rp1": ( 810.0, 246000.0, 489.0), + "ipa": ( 786.0, 665000.0, 355.6), + "nitrousoxide": (1220.0, 376000.0, 184.7), +} +_THERMO_INDEX = {"density": (0, "kg/m^3"), "latent_heat": (1, "J/kg"), "boiling_point": (2, "K")} +_GENERIC_THERMO = {"fuel": (800.0, 300000.0, 450.0), "oxidizer": (1140.0, 213000.0, 90.2)} + + +def _fluid_thermo(config, key: str, attr: str) -> float: + """``fluids[key].attr`` from the config; else the handbook value for that named fluid; else a + generic value. Both fallbacks are recorded, and the generic one says the fluid was unrecognised.""" + v = _fluid_attr(getattr(config, "fluids", None), key, attr, None) + if v is not None: + return v + from engine.pipeline.assumptions import assume + from engine.pipeline.io import _canon_fluid + idx, unit = _THERMO_INDEX[attr] + try: + f = config.fluids[key] if isinstance(config.fluids, dict) else getattr(config.fluids, key) + name = getattr(f, "name", "") or "" + except Exception: + name = "" + canon = _canon_fluid(name) + if canon in _FLUID_THERMO_FALLBACKS: + return assume(f"stability.fluids.{key}.{attr}", _FLUID_THERMO_FALLBACKS[canon][idx], unit=unit, + reason=f"fluids.{key}.{attr} missing from config; handbook value for {name}") + return assume(f"stability.fluids.{key}.{attr}", _GENERIC_THERMO["oxidizer" if key == "oxidizer" else "fuel"][idx], + unit=unit, reason=f"fluids.{key}.{attr} missing and fluid {name!r} is not in the handbook table -- set it in the config") + + +def _feed_geometry(config, side: str) -> Tuple[float, float]: + """(length [m], flow area [m^2]) of one feed line for the chug inertance L/A. + + Length comes from ``feed_system..length`` -- a field that did not exist until now, which is + why the model used to carry a hardcoded 0.305 m for every engine. Area is the schema-derived + ``A_hydraulic`` (pi d^2/4 unless the user gave a non-circular passage). + """ + from engine.pipeline.assumptions import assume + L = _feed_attr(config, side, "length", float("nan")) + if not np.isfinite(L) or L <= 0.0: + L = assume(f"stability.feed.{side}.length", 0.305, unit="m", + reason=f"feed_system.{side}.length not set (tank-outlet to manifold run)") + A = _feed_attr(config, side, "A_hydraulic", float("nan")) + if not np.isfinite(A) or A <= 0.0: + d = _feed_attr(config, side, "d_inlet", float("nan")) + if np.isfinite(d) and d > 0.0: + A = float(np.pi * (d / 2.0) ** 2) + else: + A = float(np.pi * (assume(f"stability.feed.{side}.d_inlet", 0.0127, unit="m", + reason=f"feed_system.{side} has no bore") / 2.0) ** 2) + return float(L), float(A) + + +def _chamber_dims(config, cg) -> Tuple[float, float]: + """(L_chamber, D_chamber) [m]. The solved unified geometry first: its ``length`` is the total + chamber length and its ``chamber_diameter`` is the cylindrical bore, which is what the transverse + acoustic modes live in. Legacy configs fall back to ``chamber.length`` and the volume-mean + diameter; anything still missing is a recorded assumption.""" + from engine.pipeline.assumptions import assume + L = float(getattr(cg, "length", None) or 0.0) + if L <= 0.0: + L = float(getattr(getattr(config, "chamber", None), "length", None) or 0.0) + if L <= 0.0: + L = assume("stability.chamber.length", 0.18, unit="m", reason="no solved or configured chamber length") + D = float(getattr(cg, "chamber_diameter", None) or 0.0) + if D <= 0.0: + V = float(getattr(cg, "volume", None) or 0.0) + if V > 0.0 and L > 0.0: + D = float(np.sqrt(4.0 * V / (np.pi * L))) + else: + D = assume("stability.chamber.diameter", 0.1, unit="m", reason="no chamber diameter or volume") + return L, D + + +def _stability_config(config) -> StabilityConfig: + sc = getattr(config, "stability", None) + return sc if isinstance(sc, StabilityConfig) else StabilityConfig() def build_stability_inputs(config, Pc: float, MR: float, mdot_total: float, cstar: float, @@ -383,19 +354,36 @@ def build_stability_inputs(config, Pc: float, MR: float, mdot_total: float, csta """Extract chug/acoustic model inputs from config + diagnostics. Shared by the fast path (compute_physical_stability) and the rich report (report.py) so they use IDENTICAL extraction. [Phys §3.2, §4, §5] + + Sources, in order: solved geometry (``cg``), the closure diagnostics of this evaluation, the + config (``fluids`` for the propellants, ``feed_system`` for the plumbing, ``stability`` for the + model calibration), and finally recorded assumptions -- never a silent constant. """ from engine.pipeline.stability import core, chug, acoustic + from engine.pipeline.assumptions import assume + sc = _stability_config(config) A_t = float(cg.A_throat) V_c = float(cg.volume) Lstar = V_c / A_t if A_t > 0 else float(getattr(cg, "Lstar", 0.8)) - L_ch = float(getattr(cg, "length", 0.0)) or float(getattr(config.chamber, "length", 0.18) or 0.18) - D_ch = float(np.sqrt(max(0.0, 4.0 * V_c / (np.pi * L_ch)))) if (V_c > 0 and L_ch > 0) else 0.1 + L_ch, D_ch = _chamber_dims(config, cg) + A_c = float(np.pi * (D_ch / 2.0) ** 2) + contraction_ratio = A_c / A_t if A_t > 0 else float("nan") + # Hot-gas transport properties: the same ones the thermal model uses (regen_cooling block is the + # engine's hot-gas property record, regardless of whether regen is enabled). rc = getattr(config, "regen_cooling", None) - k_g = float(getattr(rc, "hot_gas_thermal_conductivity", 0.12) or 0.12) - cp_g = float(getattr(rc, "hot_gas_cp", 0.0) or 0.0) or (gamma * R / (gamma - 1.0)) - mu_g = float(getattr(rc, "hot_gas_viscosity", 4.0e-5) or 4.0e-5) + k_g = float(getattr(rc, "hot_gas_thermal_conductivity", 0.0) or 0.0) if rc is not None else 0.0 + if k_g <= 0.0: + k_g = assume("stability.hot_gas_thermal_conductivity", DEFAULT_HOT_GAS_THERMAL_COND_W_M_K, + unit="W/(m*K)", reason="regen_cooling.hot_gas_thermal_conductivity not set") + mu_g = float(getattr(rc, "hot_gas_viscosity", 0.0) or 0.0) if rc is not None else 0.0 + if mu_g <= 0.0: + mu_g = assume("stability.hot_gas_viscosity", DEFAULT_HOT_GAS_VISC_PA_S, + unit="Pa*s", reason="regen_cooling.hot_gas_viscosity not set") + # Product-gas cp from the CEA state of THIS evaluation (gamma, R). This used to read a fixed + # regen_cooling.hot_gas_cp = 2200 J/(kg*K) for every propellant. + cp_g = gamma * R / (gamma - 1.0) rho_g = Pc / (R * Tc) if (R > 0 and Tc > 0) else 2.0 nu_g = mu_g / rho_g if rho_g > 0 else 2.0e-5 a_snd = core.sound_speed(gamma, R, Tc) @@ -418,56 +406,70 @@ def build_stability_inputs(config, Pc: float, MR: float, mdot_total: float, csta eta_O = dpiO / Pc if Pc > 0 else 0.3 eta_F = dpiF / Pc if Pc > 0 else 0.3 - from engine.pipeline.assumptions import assume - - def _fluid(key, attr, fb_value, fb_unit): - v = _fluid_attr(config.fluids, key, attr, None) - if v is not None: - return v - return assume(f"stability.fluids.{key}.{attr}", fb_value, unit=fb_unit, - reason=f"fluids.{key}.{attr} missing from config (use a propellant preset)") - - rho_O = _fluid("oxidizer", "density", 1140.0, "kg/m^3") - rho_F = _fluid("fuel", "density", 422.6, "kg/m^3") # methalox-lineage fallback — recorded - hfg_O = _fluid("oxidizer", "latent_heat", _LOX_HFG_DEFAULT, "J/kg") - hfg_F = _fluid("fuel", "latent_heat", 510000.0, "J/kg") - tbO = _fluid("oxidizer", "boiling_point", _LOX_TBOIL_DEFAULT, "K") - tbF = _fluid("fuel", "boiling_point", 111.6, "K") - K_bulk_O = _fluid("oxidizer", "bulk_modulus_pa", 1.5e9, "Pa") + rho_O = _fluid_thermo(config, "oxidizer", "density") + rho_F = _fluid_thermo(config, "fuel", "density") + hfg_O = _fluid_thermo(config, "oxidizer", "latent_heat") + hfg_F = _fluid_thermo(config, "fuel", "latent_heat") + tbO = _fluid_thermo(config, "oxidizer", "boiling_point") + tbF = _fluid_thermo(config, "fuel", "boiling_point") + K_bulk_O = _fluid_attr(config.fluids, "oxidizer", "bulk_modulus_pa", None) + if K_bulk_O is None: + K_bulk_O = assume("stability.fluids.oxidizer.bulk_modulus_pa", 1.5e9, unit="Pa", + reason="fluids.oxidizer.bulk_modulus_pa missing (set via propellant preset); measure via water-hammer test T5") tau_conv_O, _, K_v_O = core.lags_from_smd(D32_O, k_g=k_g, rho_l=rho_O, cp_g=cp_g, T_inf=Tc, T_boil=tbO, h_fg=hfg_O, chi=1.0) tau_conv_F, _, K_v_F = core.lags_from_smd(D32_F, k_g=k_g, rho_l=rho_F, cp_g=cp_g, T_inf=Tc, T_boil=tbF, h_fg=hfg_F, chi=1.0) if not np.isfinite(tau_conv_O): - tau_conv_O = 2.0e-3 + tau_conv_O = assume("stability.tau_conv_O", 2.0e-3, unit="s", + reason="d^2-law oxidizer lag non-finite (check T_boil < Tc and h_fg)") if not np.isfinite(tau_conv_F): - tau_conv_F = 1.5e-3 + tau_conv_F = assume("stability.tau_conv_F", 1.5e-3, unit="s", + reason="d^2-law fuel lag non-finite (check T_boil < Tc and h_fg)") - feed_len = 0.305 - dO = _feed_attr(config, "oxidizer", "d_inlet", 0.0135) - dF = _feed_attr(config, "fuel", "d_inlet", 0.0095) - reg_O = chug.Regulator(enabled=True) - reg_F = chug.Regulator(enabled=True) + L_feed_O, A_feed_O = _feed_geometry(config, "oxidizer") + L_feed_F, A_feed_F = _feed_geometry(config, "fuel") + reg_kw = dict(enabled=bool(sc.regulator_enabled), corner_hz=float(sc.regulator_corner_hz), + Z_hf=float(sc.regulator_Z_hf), max_excursion_pa=float(sc.regulator_max_excursion_psi) * 6894.757) streams = [ chug.ChugStream("O", mdot=mdot_O, eta_inj=max(eta_O, 1e-3), Pc=Pc, dP_feed=dpfO, - feed_length=feed_len, feed_area=np.pi * (dO / 2.0) ** 2, tau_conv=tau_conv_O, regulator=reg_O), + feed_length=L_feed_O, feed_area=A_feed_O, tau_conv=tau_conv_O, + regulator=chug.Regulator(**reg_kw)), chug.ChugStream("F", mdot=mdot_F, eta_inj=max(eta_F, 1e-3), Pc=Pc, dP_feed=dpfF, - feed_length=feed_len, feed_area=np.pi * (dF / 2.0) ** 2, tau_conv=tau_conv_F, regulator=reg_F), + feed_length=L_feed_F, feed_area=A_feed_F, tau_conv=tau_conv_F, + regulator=chug.Regulator(**reg_kw)), ] chamber = chug.ChugChamber(cstar=cstar, A_t=A_t, Lstar=Lstar, gamma=gamma) - chi_ac = float(ov.get("chi_acoustic", _CHI_ACOUSTIC_DEFAULT)) - n_int = float(ov.get("n_interaction", _N_INTERACTION_DEFAULT)) + chi_ac = float(ov.get("chi_acoustic", sc.chi_acoustic)) + n_int = float(ov.get("n_interaction", sc.n_interaction)) tau_sens = chi_ac * tau_conv_O # LOX-side rate-limiting; sensitive lag << transport lag [Phys §5] - gas = acoustic.GasState(gamma=gamma, a_sound=a_snd, nu_g=nu_g, mach_nozzle_entrance=0.2) + + # Nozzle-entrance Mach sets the convective (nozzle) damping. Config value if given, else the + # subsonic isentropic solution for the actual contraction ratio (a fixed 0.2 corresponds to a + # contraction ratio of ~2.9 and overstated nozzle damping for every wider chamber). + M_ne = sc.mach_nozzle_entrance + if M_ne is None: + M_ne = core.mach_from_area_ratio_subsonic(contraction_ratio, gamma) if np.isfinite(contraction_ratio) else float("nan") + if not np.isfinite(M_ne) or M_ne <= 0.0: + M_ne = assume("stability.mach_nozzle_entrance", 0.2, unit="-", + reason="contraction ratio unavailable for the isentropic solve") + gas = acoustic.GasState(gamma=gamma, a_sound=a_snd, nu_g=nu_g, mach_nozzle_entrance=float(M_ne)) + coeffs = acoustic.DampingCoeffs(injector_frac=float(sc.damping_injector_frac), + twophase_frac=float(sc.damping_twophase_frac), + droplet_loading=float(sc.droplet_loading)) return { - "streams": streams, "chamber": chamber, "gas": gas, - "D_ch": D_ch, "L_ch": L_ch, "Lstar": Lstar, + "streams": streams, "chamber": chamber, "gas": gas, "damping_coeffs": coeffs, + "acoustic_gate_alpha_offset": float(sc.acoustic_gate_alpha_offset), + "D_ch": D_ch, "L_ch": L_ch, "Lstar": Lstar, "contraction_ratio": contraction_ratio, + "mach_nozzle_entrance": float(M_ne), "tau_conv_O": tau_conv_O, "tau_conv_F": tau_conv_F, "tau_sens": tau_sens, "chi_acoustic": chi_ac, "n_interaction": n_int, "eta_inj_O": eta_O, "eta_inj_F": eta_F, "D32_O": D32_O, "D32_F": D32_F, "K_v_O": K_v_O, "K_v_F": K_v_F, "rho_O": rho_O, "rho_F": rho_F, "K_bulk_O": K_bulk_O, + "feed_length_O": L_feed_O, "feed_length_F": L_feed_F, + "u_O": diagnostics.get("u_O"), "Cd_O": diagnostics.get("Cd_O"), "Pc": Pc, "wh_pressure_pa": None, } @@ -500,16 +502,19 @@ def compute_physical_stability(config, Pc: float, MR: float, mdot_total: float, # difference does not earn a kernel. (The chug sweep did: 200 complex points, # measured at ~8.8% of Layer-1 wall time when left unaccelerated.) ac_fast = acoustic.fast_acoustic(inp["D_ch"], inp["L_ch"], inp["gas"], - n=inp["n_interaction"], tau_sens=inp["tau_sens"]) + n=inp["n_interaction"], tau_sens=inp["tau_sens"], + coeffs=inp["damping_coeffs"]) return { "chug": chug_fast, "acoustic": ac_fast, "chug_gate_margin": _chug_gate_margin(chug_fast.get("gain_margin", float("nan"))), - "acoustic_gate_margin": _acoustic_gate_margin(ac_fast.get("alpha_max", float("nan"))), + "acoustic_gate_margin": _acoustic_gate_margin(ac_fast.get("alpha_max", float("nan")), + inp["acoustic_gate_alpha_offset"]), "f_chug_hz": chug_fast.get("f_chug_hz"), "tau_conv_O": inp["tau_conv_O"], "tau_conv_F": inp["tau_conv_F"], "tau_sens": inp["tau_sens"], "eta_inj_O": inp["eta_inj_O"], "eta_inj_F": inp["eta_inj_F"], "D_ch": inp["D_ch"], "L_ch": inp["L_ch"], + "mach_nozzle_entrance": inp["mach_nozzle_entrance"], } @@ -551,15 +556,7 @@ def comprehensive_stability_analysis( V_chamber = float(cg.volume) A_throat = float(cg.A_throat) Lstar = V_chamber / A_throat if A_throat > 0.0 else cg.Lstar - - # Estimate chamber dimensions - L_chamber = getattr(config.chamber, "length", 0.18) or 0.18 - L_chamber = float(L_chamber) - if L_chamber <= 0.0: - L_chamber = 0.18 - - # FIXED: Add safety check for sqrt operation - D_chamber = float(np.sqrt(max(0, 4.0 * V_chamber / (np.pi * L_chamber)))) if V_chamber > 0.0 and L_chamber > 0 else 0.1 + L_chamber, D_chamber = _chamber_dims(config, cg) # Combustion stability chugging = calculate_chugging_frequency( @@ -580,27 +577,11 @@ def comprehensive_stability_analysis( R=R, ) - # Feed system stability (use LOX feed as representative) - if getattr(config, "feed_system", None) is not None: - if isinstance(config.feed_system, dict): - lox_config = config.feed_system.get("lox", {}) - if isinstance(lox_config, dict): - feed_length = float(lox_config.get("length", 1.0)) - feed_diameter = float(lox_config.get("d_inlet", 0.01)) - else: - feed_length = float(getattr(lox_config, "length", 1.0)) - feed_diameter = float(getattr(lox_config, "d_inlet", 0.01)) - else: - lox_config = getattr(config.feed_system, "lox", None) - if lox_config is not None: - feed_length = float(getattr(lox_config, "length", 1.0)) - feed_diameter = float(getattr(lox_config, "d_inlet", 0.01)) - else: - feed_length = 1.0 - feed_diameter = 0.01 - else: - feed_length = 1.0 - feed_diameter = 0.01 + # Feed-line acoustics on the oxidizer line (the stiffer, denser side; representative). Length + # and bore come from feed_system.oxidizer -- the old lookup asked for a "lox" branch and a + # "length" attribute that never existed, so it always fell through to 1.0 m x 10 mm. + feed_length, A_feed = _feed_geometry(config, "oxidizer") + feed_diameter = float(np.sqrt(4.0 * A_feed / np.pi)) # Oxidizer density / bulk modulus from config.fluids (the old `config.propellants` lookup was a # dead key — it ALWAYS fell through to 1140. UNIFICATION P2c: config-first, recorded fallback.) @@ -615,8 +596,7 @@ def comprehensive_stability_analysis( bulk_modulus = assume("stability.feed.bulk_modulus_O", 1.5e9, unit="Pa", reason="fluids.oxidizer.bulk_modulus_pa missing (set via propellant preset); measure via water-hammer test T5") - # Estimate oxidizer flow velocity - A_feed = np.pi * (feed_diameter / 2.0) ** 2 + # Mean oxidizer line velocity mdot_ox = float(diagnostics.get("mdot_O", mdot_total * MR / (1.0 + MR))) flow_velocity = float(mdot_ox / (prop_density * A_feed)) if A_feed > 0.0 else 0.0 @@ -640,6 +620,40 @@ def comprehensive_stability_analysis( for i, freq in enumerate(acoustic_raw["transverse_modes"]): acoustic_modes_dict[f"T{i+1}"] = freq + issues: List[str] = [] + + # ------------------------------------------------------------------- + # Physical margins (new model) — replaces the heuristic score/margins. [plan A4, M4] + # ------------------------------------------------------------------- + phys = None + try: + phys = compute_physical_stability(config, Pc, MR, mdot_total, cstar, gamma, R, Tc, diagnostics, cg) + except Exception: # defensive: never fail the eval on a stability-model error + phys = None + + if phys is not None: + chug_margin = float(phys["chug_gate_margin"]) + acoustic_margin = float(phys["acoustic_gate_margin"]) + _fch = phys.get("f_chug_hz") + if _fch is not None and np.isfinite(_fch) and _fch > 0: + chugging["frequency"] = float(_fch) # physical chug freq, not L*/c* placeholder + chugging["stability_margin"] = chug_margin + chugging["chug_gain_margin"] = phys["chug"].get("gain_margin") + feed_stability["stability_margin"] = chug_margin # feed-coupled instability IS chug (un-rig) + if not phys["chug"].get("stable", True): + issues.append("Chug (feed-coupled LF) margin low: stiffen injector or improve atomization") + if not phys["acoustic"].get("stable", True): + issues.append(f"Acoustic mode {phys['acoustic'].get('limiting_mode')} driven (alpha>0)") + else: + # The physical model could not be evaluated: margins are UNKNOWN. Neutral-pass so a + # stability-model error never fails an evaluation, and say so in the issues list rather + # than reporting a heuristic as if it were a margin. + chug_margin = 1.10 + acoustic_margin = 1.10 + chugging["stability_margin"] = chug_margin + feed_stability["stability_margin"] = chug_margin + issues.append("Stability model could not be evaluated for this point; margins shown are neutral placeholders") + # ------------------------------------------------------------------- # Mode coupling analysis # ------------------------------------------------------------------- @@ -647,7 +661,9 @@ def comprehensive_stability_analysis( # Collect representative modes for coupling checks modes: List[Dict[str, Any]] = [] - modes.append({"name": "chugging", "type": "combustion", "frequency": chugging["frequency"]}) + # The physical chug frequency when the model ran; the geometric placeholder is not a chug mode. + if phys is not None and np.isfinite(chugging["frequency"]): + modes.append({"name": "chugging", "type": "combustion", "frequency": chugging["frequency"]}) modes.append({"name": "pogo", "type": "feed", "frequency": feed_stability["pogo_frequency"]}) modes.append({"name": "surge", "type": "feed", "frequency": feed_stability["surge_frequency"]}) @@ -683,9 +699,7 @@ def comprehensive_stability_analysis( # Stability classification # ------------------------------------------------------------------- - issues: List[str] = [] - - # NOTE: chug/acoustic issues come from the PHYSICAL model below (not the old heuristic chugging + # NOTE: chug/acoustic issues come from the PHYSICAL model above (not the old heuristic chugging # stability_index), and water-hammer is handled below as a separate valve-transient note. # Mode coupling @@ -696,34 +710,6 @@ def comprehensive_stability_analysis( if Lstar < 0.5 or Lstar > 3.0: issues.append(f"L* outside typical range (0.5 m to 3.0 m). Current L* = {Lstar:.2f} m") - # ------------------------------------------------------------------- - # Physical margins (new model) — replaces the heuristic score/margins. [plan A4, M4] - # ------------------------------------------------------------------- - phys = None - try: - phys = compute_physical_stability(config, Pc, MR, mdot_total, cstar, gamma, R, Tc, diagnostics, cg) - except Exception: # defensive: never fail the eval on a stability-model error - phys = None - - if phys is not None: - chug_margin = float(phys["chug_gate_margin"]) - acoustic_margin = float(phys["acoustic_gate_margin"]) - _fch = phys.get("f_chug_hz") - if _fch is not None and np.isfinite(_fch) and _fch > 0: - chugging["frequency"] = float(_fch) # physical chug freq, not L*/c* placeholder - chugging["stability_margin"] = chug_margin - chugging["chug_gain_margin"] = phys["chug"].get("gain_margin") - feed_stability["stability_margin"] = chug_margin # feed-coupled instability IS chug (un-rig) - if not phys["chug"].get("stable", True): - issues.append("Chug (feed-coupled LF) margin low: stiffen injector or improve atomization") - if not phys["acoustic"].get("stable", True): - issues.append(f"Acoustic mode {phys['acoustic'].get('limiting_mode')} driven (alpha>0)") - else: - # fallback: do NOT regress if extraction fails — neutral-pass margins - chug_margin = float(chugging.get("stability_margin", 1.10)) - acoustic_margin = 1.10 - feed_stability["stability_margin"] = chug_margin - # Numeric score in [0,1] monotone in the limiting gate margin (1.05 ~ gate threshold). min_margin = min(chug_margin, acoustic_margin) score = float(np.clip((min_margin - 0.85) / 0.45, 0.0, 1.0)) @@ -793,9 +779,10 @@ def _generate_stability_recommendations( else: recs.append("System appears reasonably stable for this point. Still monitor during hot fire.") - # Chugging related - if chugging["stability_index"] < 0.5: - recs.append("Increase chamber pressure or L* to improve low frequency combustion stability.") + # Chug: Nyquist gain margin of the feed-coupled loop (>1 stable; <1.5 is thin) + gm = chugging.get("chug_gain_margin") + if gm is not None and np.isfinite(gm) and gm < 1.5: + recs.append("Chug gain margin is thin: stiffen the injector (raise dP_inj/Pc) or shorten the vaporization lag (finer SMD).") recs.append("Consider injector or chamber damping features such as baffles or acoustic liners.") if chugging["frequency"] < 10.0: diff --git a/EngineDesign/engine/pipeline/stability/core.py b/EngineDesign/engine/pipeline/stability/core.py index 531d1b0b9..9ec58eedd 100644 --- a/EngineDesign/engine/pipeline/stability/core.py +++ b/EngineDesign/engine/pipeline/stability/core.py @@ -29,6 +29,8 @@ "n_tau_gain", "choked_flow_function", "chamber_residence_time", + "area_ratio_from_mach", + "mach_from_area_ratio_subsonic", "spalding_transfer_number_heat", "d2_law_evaporation_constant", "vaporization_time", @@ -135,6 +137,35 @@ def chamber_residence_time(Lstar: float, cstar: float, gamma: float) -> float: return float(Lstar / (G * G * cstar)) +def area_ratio_from_mach(M: float, gamma: float) -> float: + """Isentropic area ratio ``A/A* = (1/M) * [(2/(g+1)) * (1 + (g-1)/2 * M^2)]^((g+1)/(2(g-1)))``.""" + g = float(gamma) + if M <= 0 or g <= 1.0: + return float("nan") + return float((1.0 / M) * ((2.0 / (g + 1.0)) * (1.0 + 0.5 * (g - 1.0) * M * M)) ** ((g + 1.0) / (2.0 * (g - 1.0)))) + + +def mach_from_area_ratio_subsonic(area_ratio: float, gamma: float) -> float: + """Subsonic Mach number at a station with ``A/A* = area_ratio`` (isentropic, one-dimensional). + + This is the mean Mach at the nozzle entrance when ``area_ratio`` is the contraction ratio + ``A_chamber / A_throat``, which is what sets the convective (nozzle) acoustic damping. Bisection + on [1e-6, 1]: A/A* is monotone decreasing in M on the subsonic branch. ``area_ratio <= 1`` -> 1.0. + """ + if not np.isfinite(area_ratio) or gamma <= 1.0: + return float("nan") + if area_ratio <= 1.0: + return 1.0 + lo, hi = 1e-6, 1.0 + for _ in range(80): + mid = 0.5 * (lo + hi) + if area_ratio_from_mach(mid, gamma) > area_ratio: + lo = mid # too subsonic: area ratio still above target -> raise M + else: + hi = mid + return float(0.5 * (lo + hi)) + + # --------------------------------------------------------------------------- # 4. Vaporization / time lag # --------------------------------------------------------------------------- diff --git a/EngineDesign/engine/pipeline/stability/enhanced.py b/EngineDesign/engine/pipeline/stability/enhanced.py deleted file mode 100644 index ee4bce8e3..000000000 --- a/EngineDesign/engine/pipeline/stability/enhanced.py +++ /dev/null @@ -1,383 +0,0 @@ -"""Enhanced physics-based stability analysis for pintle injectors. - -Accounts for: -1. Pintle geometry (tip diameter, length, gap) -2. Fuel impingement zones (localized instability sources) -3. Recirculation zones (flow patterns near pintle tip) -4. Pintle length effects (acoustic coupling) -5. Uneven ablation (spatial variation in geometry) -6. Real wave propagation physics with proper boundary conditions -""" - -from __future__ import annotations - -from typing import Dict, List, Tuple, Optional -import numpy as np -from engine.pipeline.config_schemas import PintleEngineConfig, PintleInjectorConfig -from engine.pipeline.localized_ablation import calculate_impingement_zones - - -def calculate_pintle_recirculation_zones( - L_pintle: float, - d_pintle_tip: float, - D_chamber: float, - L_chamber: float, - positions: np.ndarray, - fuel_velocity: float = 50.0, - lox_velocity: float = 30.0, -) -> Dict[str, np.ndarray]: - """ - Calculate recirculation zones near pintle tip. - - Physics: - - Fuel spray from pintle tip creates recirculation eddies - - LOX jets create additional recirculation - - Recirculation zones have different acoustic properties - - These zones affect wave propagation and stability - - Parameters: - ----------- - L_pintle : float - Pintle length [m] (distance from injector face to tip) - d_pintle_tip : float - Pintle tip diameter [m] - D_chamber : float - Chamber diameter [m] - L_chamber : float - Chamber length [m] - positions : np.ndarray - Axial positions [m] - fuel_velocity : float - Fuel injection velocity [m/s] - lox_velocity : float - LOX injection velocity [m/s] - - Returns: - -------- - recirculation : dict - - recirculation_intensity: Local recirculation intensity (0-1) - - recirculation_length: Characteristic recirculation length [m] - - velocity_fluctuation: Velocity fluctuation magnitude [m/s] - - turbulence_intensity: Turbulence intensity (0-1) - """ - n_points = len(positions) - - # Recirculation zone extends from injector face (x=0) to ~2-3x pintle length - L_recirc = 2.5 * L_pintle # Typical recirculation length - - # Recirculation intensity decays with distance from pintle tip - recirculation_intensity = np.zeros(n_points) - recirculation_length = np.zeros(n_points) - velocity_fluctuation = np.zeros(n_points) - turbulence_intensity = np.zeros(n_points) - - for i, x in enumerate(positions): - if x <= L_recirc: - # Recirculation zone: intensity decays exponentially - decay_factor = np.exp(-x / (0.5 * L_pintle)) - - # Physics-based recirculation intensity - from engine.pipeline.physics_based_replacements import calculate_recirculation_intensity_physics - - # Estimate Reynolds number - rho_approx = 5.0 # kg/m³, typical hot gas - mu_approx = 4e-5 # Pa·s - Re_injector = rho_approx * fuel_velocity * d_pintle_tip / mu_approx - - base_intensity = calculate_recirculation_intensity_physics( - fuel_velocity=fuel_velocity, - lox_velocity=lox_velocity, - d_pintle_tip=d_pintle_tip, - D_chamber=D_chamber, - Re_injector=Re_injector, - ) - - recirculation_intensity[i] = base_intensity * decay_factor - - # Characteristic recirculation length (eddy size) - # From turbulent mixing theory: L_eddy ~ 0.1-0.3 × injector size - # Depends on velocity ratio and Reynolds number - eddy_base = 0.2 * d_pintle_tip # Base eddy size - velocity_factor = 1.0 + 0.3 * (fuel_velocity / (lox_velocity + 1e-10) - 1.0) - recirculation_length[i] = eddy_base * velocity_factor * (1.0 + 0.2 * decay_factor) - - # Velocity fluctuations (RMS) from turbulence theory - # u' ~ 0.1-0.2 × U for turbulent flow - # Higher recirculation → higher fluctuations - v_fluct_base = 0.12 * fuel_velocity * (1.0 + base_intensity) # Physics-based - velocity_fluctuation[i] = v_fluct_base * decay_factor - - # Turbulence intensity from mixing theory - # I_turb ~ 0.1-0.3 for recirculating flows - # Depends on velocity ratio and recirculation intensity - base_turbulence = 0.1 + 0.1 * base_intensity # Physics-based - velocity_enhancement = 1.0 + 0.2 * (fuel_velocity / (lox_velocity + 1e-10) - 1.0) - turbulence_intensity[i] = base_turbulence * decay_factor * velocity_enhancement - else: - # Outside recirculation zone - recirculation_intensity[i] = 0.0 - recirculation_length[i] = 0.0 - velocity_fluctuation[i] = 0.0 - turbulence_intensity[i] = 0.05 # Base turbulence - - return { - "recirculation_intensity": recirculation_intensity, - "recirculation_length": recirculation_length, - "velocity_fluctuation": velocity_fluctuation, - "turbulence_intensity": turbulence_intensity, - } - - -def calculate_pintle_stability_enhanced( - config: PintleEngineConfig, - positions: np.ndarray, - chamber_pressure: np.ndarray, - sound_speed: np.ndarray, - density: np.ndarray, - mass_flow: np.ndarray, - recession_profile: Optional[np.ndarray] = None, - L_chamber: float = 0.2, - D_chamber: float = 0.1, - fuel_velocity: float = 50.0, - lox_velocity: float = 30.0, -) -> Dict[str, np.ndarray]: - """ - Enhanced stability calculation with full pintle physics. - - Physics: - 1. Pintle geometry affects injector impedance and acoustic coupling - 2. Fuel impingement creates localized pressure fluctuation sources - 3. Recirculation zones create acoustic damping/amplification - 4. Pintle length affects acoustic mode coupling - 5. Uneven ablation creates impedance mismatches - 6. Wave propagation with proper boundary conditions - - Parameters: - ----------- - config : PintleEngineConfig - Engine configuration - positions : np.ndarray - Axial positions [m] - chamber_pressure : np.ndarray - Local pressure [Pa] - sound_speed : np.ndarray - Local sound speed [m/s] - density : np.ndarray - Local density [kg/m³] - mass_flow : np.ndarray - Local mass flow [kg/s] - recession_profile : np.ndarray, optional - Local recession [m] at each position (for uneven ablation) - L_chamber : float - Chamber length [m] - D_chamber : float - Chamber diameter [m] - fuel_velocity : float - Fuel injection velocity [m/s] - lox_velocity : float - LOX injection velocity [m/s] - - Returns: - -------- - stability : dict - - chugging_frequency: Local chugging frequency [Hz] - - stability_margin: Local stability margin - - wave_growth_rate: Wave growth rate [1/s] - - impingement_effect: Effect of impingement on stability - - recirculation_effect: Effect of recirculation on stability - - ablation_effect: Effect of uneven ablation on stability - - pintle_length_effect: Effect of pintle length on acoustic coupling - """ - n_points = len(positions) - - # Get pintle geometry - if not hasattr(config, 'injector') or config.injector.type != "pintle": - # No pintle-specific effects - return { - "chugging_frequency": np.full(n_points, 30.0), - "stability_margin": np.full(n_points, 0.5), - "wave_growth_rate": np.full(n_points, -10.0), - "impingement_effect": np.zeros(n_points), - "recirculation_effect": np.zeros(n_points), - "ablation_effect": np.zeros(n_points), - "pintle_length_effect": np.zeros(n_points), - } - - injector_config: PintleInjectorConfig = config.injector - geometry = injector_config.geometry - - # Pintle geometry parameters - d_pintle_tip = geometry.fuel.d_pintle_tip - h_gap = geometry.fuel.h_gap - L_pintle = getattr(geometry.fuel, 'L_pintle', 0.01) # Pintle length [m] - n_orifices = geometry.lox.n_orifices - d_orifice = geometry.lox.d_orifice - theta_orifice = geometry.lox.theta_orifice - - # Calculate impingement zones (where fuel hits wall) - impingement_data = calculate_impingement_zones( - config, L_chamber, D_chamber, n_points=n_points - ) - impingement_multiplier = impingement_data["impingement_heat_flux_multiplier"] - impingement_zones = impingement_data["impingement_zones"] - impingement_center = impingement_data.get("impingement_center", L_chamber * 0.7) - - # Calculate recirculation zones (near pintle tip) - recirculation_data = calculate_pintle_recirculation_zones( - L_pintle, d_pintle_tip, D_chamber, L_chamber, positions, - fuel_velocity, lox_velocity - ) - recirculation_intensity = recirculation_data["recirculation_intensity"] - recirculation_length = recirculation_data["recirculation_length"] - velocity_fluctuation = recirculation_data["velocity_fluctuation"] - turbulence_intensity = recirculation_data["turbulence_intensity"] - - # Calculate injector impedance from pintle geometry - # Acoustic impedance: Z = ρ × c / A - A_pintle_tip = np.pi * (d_pintle_tip / 2.0) ** 2 - A_gap = np.pi * d_pintle_tip * h_gap # Annular gap area - A_injector_effective = A_pintle_tip + A_gap - - # Injector impedance (at injection plane) - rho_injector = density[0] if len(density) > 0 else 1000.0 - c_injector = sound_speed[0] if len(sound_speed) > 0 else 1000.0 - Z_injector = rho_injector * c_injector / A_injector_effective if A_injector_effective > 0 else 1e6 - - # Feed system impedance (simplified) - Z_feed = 5e5 # Typical feed system impedance [Pa·s/m³] - - # Calculate local impedances - A_local = np.pi * (D_chamber / 2.0) ** 2 - Z_local = density * sound_speed / A_local - - # Wave propagation time - L_total = positions[-1] - positions[0] if len(positions) > 1 else L_chamber - tau_wave = L_total / sound_speed # Wave propagation time - - # Base chugging frequency from wave resonance - # f = c / (4L) for open-closed tube (injector closed, throat open) - f_chugging_base = sound_speed / (4.0 * L_total) - - # Pintle length effect on frequency - # Longer pintle = different acoustic coupling = frequency shift - # Pintle acts as acoustic extension of injector - L_effective = L_total + 0.3 * L_pintle # Effective length includes pintle - f_chugging_pintle = sound_speed / (4.0 * L_effective) - - # Frequency shift from pintle geometry - pintle_ratio = d_pintle_tip / D_chamber if D_chamber > 0 else 0.1 - frequency_shift = 1.0 + 0.15 * (pintle_ratio - 0.1) + 0.1 * (L_pintle / L_chamber) - f_chugging = f_chugging_pintle * frequency_shift - - # Pintle length effect on acoustic coupling - # Longer pintle = stronger coupling between injector and chamber - coupling_strength = 1.0 + 0.5 * (L_pintle / L_chamber) # Stronger coupling - pintle_length_effect = (coupling_strength - 1.0) * 0.3 # Can be stabilizing or destabilizing - - # Impingement effect on stability - # Fuel impingement creates localized pressure fluctuation sources - # These act as instability sources - impingement_effect = np.zeros(n_points) - for i, (pos, is_impingement) in enumerate(zip(positions, impingement_zones)): - if is_impingement: - # Impingement creates pressure fluctuation source - # Effect decays with distance from impingement - distance_from_impingement = abs(pos - impingement_center) - decay_factor = np.exp(-distance_from_impingement / (L_chamber * 0.1)) - # Impingement multiplier indicates intensity - intensity = (impingement_multiplier[i] - 1.0) * 0.5 # Destabilizing - impingement_effect[i] = intensity * decay_factor - - # Recirculation effect on stability - # Recirculation zones can: - # 1. Damp waves (turbulence dissipation) - # 2. Amplify waves (resonance in eddies) - # Net effect depends on recirculation intensity and turbulence - recirculation_effect = np.zeros(n_points) - for i in range(n_points): - if recirculation_intensity[i] > 0: - # Recirculation creates velocity fluctuations - # These can couple with pressure waves - # High turbulence = damping (stabilizing) - # Low turbulence + high intensity = amplification (destabilizing) - turbulence_damping = turbulence_intensity[i] * 0.5 # Stabilizing - recirculation_amplification = recirculation_intensity[i] * (1.0 - turbulence_intensity[i]) * 0.3 # Destabilizing - recirculation_effect[i] = recirculation_amplification - turbulence_damping - - # Uneven ablation effect - # Spatial variation in geometry creates impedance mismatches - # These reflect waves and can cause instability - ablation_effect = np.zeros(n_points) - if recession_profile is not None and len(recession_profile) == n_points: - # Calculate local diameter variation - D_local = D_chamber + 2.0 * recession_profile - A_local_varying = np.pi * (D_local / 2.0) ** 2 - - # Impedance variation - Z_varying = density * sound_speed / A_local_varying - - # Impedance mismatch creates reflections - # Large mismatch = more reflections = potential instability - Z_ref = Z_local[0] if len(Z_local) > 0 else Z_local.mean() - impedance_mismatch = np.abs(Z_varying - Z_ref) / (Z_ref + 1e-10) - ablation_effect = impedance_mismatch * 0.3 # Destabilizing effect - - # Wave growth rate from energy balance - # Energy input from combustion vs. energy dissipation - energy_input = chamber_pressure * mass_flow / density # [W/m³] - energy_dissipation = density * (mass_flow / (density * A_local)) ** 2 / L_total # [W/m³] - energy_stored = 0.5 * density * sound_speed ** 2 # [J/m³] - - # Base growth rate from energy balance - # For well-designed engines, dissipation > input (net damping) - energy_balance = (energy_input - energy_dissipation) / (2.0 * energy_stored + 1e-10) - - # Base damping rate: well-designed engines have negative growth (damping) - # Typical damping: -20 to -100 [1/s] for stable engines - # Add base damping from acoustic losses, wall friction, etc. - base_damping = -50.0 # [1/s] - base damping rate (negative = stable) - - # Energy balance modifies base damping - # Positive energy balance (input > dissipation) = destabilizing - # Negative energy balance (dissipation > input) = stabilizing - wave_growth_base = base_damping + energy_balance * 10.0 # Scale energy balance effect - - # Add all effects - # Impingement: destabilizing (reduces damping) - # Recirculation: can be stabilizing (turbulence damping) or destabilizing (amplification) - # Ablation: destabilizing (impedance mismatch) - # Pintle length: usually stabilizing (better mixing = more damping) - wave_growth_rate = ( - wave_growth_base - + impingement_effect # Destabilizing (reduces damping) - + recirculation_effect # Can be stabilizing or destabilizing - - ablation_effect * 0.5 # Destabilizing (impedance mismatch) - - pintle_length_effect * 0.3 # Usually stabilizing (better mixing) - ) - - # Stability margin: positive = stable, negative = unstable - # For stability: wave_growth_rate should be negative (damping) - # Margin = (damping_rate - growth_rate) / reference_rate - # Higher margin = more stable - reference_rate = 100.0 # [1/s] - reference growth rate - damping_rate = -wave_growth_rate # Convert growth to damping - stability_margin = damping_rate / reference_rate - - # Clamp to reasonable range: -2 to +2 - # Positive = stable, negative = unstable - stability_margin = np.clip(stability_margin, -2.0, 2.0) - - return { - "chugging_frequency": f_chugging, - "stability_margin": stability_margin, - "wave_growth_rate": wave_growth_rate, - "impingement_effect": impingement_effect, - "recirculation_effect": recirculation_effect, - "ablation_effect": ablation_effect, - "pintle_length_effect": pintle_length_effect, - "recirculation_intensity": recirculation_intensity, - "turbulence_intensity": turbulence_intensity, - "pintle_ratio": pintle_ratio, - "frequency_shift": frequency_shift, - } - diff --git a/EngineDesign/engine/pipeline/stability/report.py b/EngineDesign/engine/pipeline/stability/report.py index a6f1c64bf..0a008d572 100644 --- a/EngineDesign/engine/pipeline/stability/report.py +++ b/EngineDesign/engine/pipeline/stability/report.py @@ -64,8 +64,15 @@ def _vaporization_profile(inp: Dict[str, Any], Pc: float, n_pts: int = 40) -> Di L_ch = inp["L_ch"] rho_O = float(inp.get("rho_O", 1140.0)) # config-sourced via build_stability_inputs (P2c) eta = inp["eta_inj_O"] - # representative droplet axial speed ~ LOX injection velocity v=sqrt(2*dP/rho) (Cd~0.6) - v_drop = 0.6 * float(np.sqrt(max(2.0 * eta * Pc / rho_O, 1.0))) + # Representative droplet axial speed: the solved oxidizer injection velocity when the closure + # provides it, else Bernoulli with the solved Cd (a fixed Cd of 0.6 used to sit here). + u_O = inp.get("u_O") + if u_O is not None and np.isfinite(float(u_O)) and float(u_O) > 0.0: + v_drop = float(u_O) + else: + Cd = inp.get("Cd_O") + Cd = float(Cd) if (Cd is not None and np.isfinite(float(Cd)) and float(Cd) > 0.0) else 0.6 + v_drop = Cd * float(np.sqrt(max(2.0 * eta * Pc / rho_O, 1.0))) tau_vap = inp["tau_conv_O"] L_vap = v_drop * tau_vap if np.isfinite(tau_vap) else float("nan") x_max = float(max(L_ch, L_vap if np.isfinite(L_vap) else L_ch) * 1.1) @@ -83,11 +90,11 @@ def _vaporization_profile(inp: Dict[str, Any], Pc: float, n_pts: int = 40) -> Di def _sensitivity(inp: Dict[str, Any]) -> Dict[str, Any]: """n / chi sensitivity bands for the acoustic limiting-mode growth rate (cheap sweep).""" - D_ch, L_ch, gas = inp["D_ch"], inp["L_ch"], inp["gas"] + D_ch, L_ch, gas, coeffs = inp["D_ch"], inp["L_ch"], inp["gas"], inp["damping_coeffs"] tv = inp["tau_conv_O"] - a_n = [acoustic.fast_acoustic(D_ch, L_ch, gas, n=nn, tau_sens=inp["tau_sens"])["alpha_max"] + a_n = [acoustic.fast_acoustic(D_ch, L_ch, gas, n=nn, tau_sens=inp["tau_sens"], coeffs=coeffs)["alpha_max"] for nn in (0.3, 0.6)] - a_chi = [acoustic.fast_acoustic(D_ch, L_ch, gas, n=inp["n_interaction"], tau_sens=cc * tv)["alpha_max"] + a_chi = [acoustic.fast_acoustic(D_ch, L_ch, gas, n=inp["n_interaction"], tau_sens=cc * tv, coeffs=coeffs)["alpha_max"] for cc in (0.05, 0.30)] return {"acoustic_alpha_vs_n": [float(min(a_n)), float(max(a_n))], "acoustic_alpha_vs_chi": [float(min(a_chi)), float(max(a_chi))]} @@ -103,7 +110,7 @@ def _chug_pole(chug_rich: Dict[str, Any]) -> Dict[str, float]: def _radar(chug_margin: float, ac: Dict[str, Any], vap: Dict[str, Any], - gate_threshold: float) -> Dict[str, Any]: + gate_threshold: float, alpha_offset: float) -> Dict[str, Any]: """Viz #7: one-glance health radar.""" def mode_alpha(name): for m in ac["modes"]: @@ -112,8 +119,8 @@ def mode_alpha(name): return float("-inf") a1L, a1T = mode_alpha("1L"), mode_alpha("1T") # normalize alphas to a 0..1.3 "margin-like" scale via the same acoustic gate mapping - v1L = analysis._acoustic_gate_margin(a1L) - v1T = analysis._acoustic_gate_margin(a1T) + v1L = analysis._acoustic_gate_margin(a1L, alpha_offset) + v1T = analysis._acoustic_gate_margin(a1T, alpha_offset) vap_complete = float(np.clip(vap["L_ch_m"] / vap["L_vap_m"], 0.0, 1.3)) if ( np.isfinite(vap["L_vap_m"]) and vap["L_vap_m"] > 0) else 1.3 axes = ["chug", "1L", "1T", "vaporization"] @@ -272,9 +279,10 @@ def build_rich_report(config, Pc: float, MR: float, mdot_total: float, cstar: fl # --- acoustic (full mode set with damping budgets) --- ac = acoustic.analyze_acoustic_modes(inp["D_ch"], inp["L_ch"], gas, - n=inp["n_interaction"], tau_sens=inp["tau_sens"]) + n=inp["n_interaction"], tau_sens=inp["tau_sens"], + coeffs=inp["damping_coeffs"]) ac_alpha_max = ac["modes"][0]["alpha"] if ac["modes"] else float("nan") - acoustic_margin = analysis._acoustic_gate_margin(ac_alpha_max) + acoustic_margin = analysis._acoustic_gate_margin(ac_alpha_max, inp["acoustic_gate_alpha_offset"]) acoustic_modes = [{ "name": m["mode"], "freq_hz": m["f_hz"], "alpha": m["alpha"], "driving": m["driving"], @@ -288,7 +296,7 @@ def build_rich_report(config, Pc: float, MR: float, mdot_total: float, cstar: fl vap = _vaporization_profile(inp, Pc) sens = _sensitivity(inp) - radar = _radar(chug_margin, ac, vap, gate_threshold) + radar = _radar(chug_margin, ac, vap, gate_threshold, inp["acoustic_gate_alpha_offset"]) min_margin = float(min(chug_margin, acoustic_margin)) state = ("stable" if (chug_margin >= gate_threshold and acoustic_margin >= gate_threshold @@ -320,6 +328,10 @@ def build_rich_report(config, Pc: float, MR: float, mdot_total: float, cstar: fl "dP_reg_max_psi": float(streams[0].regulator.max_excursion_pa / _PA_PER_PSI), "eta_inj_O": inp["eta_inj_O"], "eta_inj_F": inp["eta_inj_F"], "smd_O_um": float(inp["D32_O"] * 1e6), + "mach_nozzle_entrance": float(inp["mach_nozzle_entrance"]), + "contraction_ratio": float(inp["contraction_ratio"]), + "feed_length_O_m": float(inp["feed_length_O"]), "feed_length_F_m": float(inp["feed_length_F"]), + "acoustic_gate_alpha_offset": float(inp["acoustic_gate_alpha_offset"]), # Every recorded silent-default substitution this process has made (P2c registry). # Empty list = config fully specified the physics. The hardcoded-Cd bug class, surfaced. "fallbacks_used": _fallbacks_used(), diff --git a/EngineDesign/engine/pipeline/time_varying_solver.py b/EngineDesign/engine/pipeline/time_varying_solver.py index ea11d740d..fd367a861 100644 --- a/EngineDesign/engine/pipeline/time_varying_solver.py +++ b/EngineDesign/engine/pipeline/time_varying_solver.py @@ -35,7 +35,6 @@ from engine.pipeline.stability.analysis import ( calculate_chugging_frequency, calculate_acoustic_modes, - analyze_feed_system_stability, # Correct function name ) from engine.pipeline.thermal.regen_cooling import estimate_hot_wall_heat_flux @@ -718,109 +717,34 @@ def solve_time_step( "hotspot_max_intensity": 1.0, } - # Calculate stability with pintle geometry, impingement, and recirculation - # Use enhanced physics-based spatial stability analysis - # NOTE (UNIFICATION P3, FINDING F3): calculate_pintle_stability_enhanced returns a flat - # PLACEHOLDER (margin 0.5, 30 Hz) for non-pintle injectors. The scalars stability_margin and - # chugging_freq are overwritten below by the injector-agnostic comprehensive_stability_analysis, - # but feed_stability["stability_margin"] and the spatial `acoustic` retain placeholder values - # for non-pintle. Clean separation deferred to the legacy_pintle/ refactor (see CONTEXT.md). - try: - from engine.pipeline.stability.enhanced import calculate_pintle_stability_enhanced - from engine.pipeline.localized_ablation import calculate_impingement_zones - - # Create position array for spatial analysis - n_stability_points = 50 - positions_stability = np.linspace(0.0, self.L_chamber, n_stability_points) - - # Calculate local properties (simplified - assume uniform for now) - P_local = np.full(n_stability_points, Pc) - c_local = np.full(n_stability_points, np.sqrt(gamma_chamber * R_chamber * Tc)) - rho_local = np.full(n_stability_points, Pc / (R_chamber * Tc)) - mdot_local = np.full(n_stability_points, mdot_total) - - # Recession profile (spatial variation) - recession_profile = None - if ablative_cfg and ablative_cfg.enabled: - # Create spatial recession profile (more at impingement zones) - impingement_data = calculate_impingement_zones( - config_current, self.L_chamber, D_chamber_new, n_points=n_stability_points - ) - # Recession is enhanced at impingement zones - recession_base = recession_chamber_new - recession_profile = recession_base * impingement_data["impingement_heat_flux_multiplier"] - - # Get injection velocities for recirculation calculation - # These would come from injector solve, but use estimates for now - fuel_velocity = 50.0 # [m/s] - typical fuel injection velocity - lox_velocity = 30.0 # [m/s] - typical LOX injection velocity - - # Calculate enhanced pintle-based stability with recirculation - stability_spatial = calculate_pintle_stability_enhanced( - config_current, - positions_stability, - P_local, - c_local, - rho_local, - mdot_local, - recession_profile=recession_profile, - L_chamber=self.L_chamber, - D_chamber=D_chamber_new, - fuel_velocity=fuel_velocity, - lox_velocity=lox_velocity, - ) - - # Use average values for single-point metrics - chugging_freq = float(np.mean(stability_spatial["chugging_frequency"])) - stability_margin = float(np.mean(stability_spatial["stability_margin"])) - - # Acoustic modes (use base calculation for now, could be enhanced) - acoustic = calculate_acoustic_modes( - self.L_chamber, - D_chamber_new, - gamma_chamber, - R_chamber, - Tc, - ) - - # Feed system stability - feed_stability = { - "pogo_frequency": np.nan, - "surge_frequency": np.nan, - "stability_margin": stability_margin, - } - - except Exception as e: - # Fallback to simple calculation - import warnings - warnings.warn(f"Pintle stability calculation failed, using fallback: {e}") - chugging = calculate_chugging_frequency( - V_chamber_new, - A_throat_new, - cstar_actual, - gamma_chamber, - Pc, - R=R_chamber, - Tc=Tc, - ) - chugging_freq = chugging["frequency"] - # CRITICAL FIX: Remove arbitrary 0.5 default - stability margin should be calculated - # If not available, use neutral (0.0) rather than arbitrary positive value - stability_margin = chugging.get("stability_margin", 0.0) # Neutral if unknown - - acoustic = calculate_acoustic_modes( - self.L_chamber, - D_chamber_new, - Tc, - gamma_chamber, - R_chamber, - ) - - feed_stability = { - "pogo_frequency": np.nan, - "surge_frequency": np.nan, - "stability_margin": 1.0, - } + # Stability. The injector-agnostic comprehensive analysis below is authoritative; these + # are the placeholders it overwrites, kept so the state record is always populated even + # when that analysis raises. (The old pintle-only "enhanced" spatial model that used to + # run here was fed hardcoded 50/30 m/s injection velocities and invented damping + # constants, and every scalar it produced was overwritten anyway -- removed.) + chugging = calculate_chugging_frequency( + V_chamber_new, + A_throat_new, + cstar_actual, + gamma_chamber, + Pc, + R=R_chamber, + Tc=Tc, + ) + chugging_freq = chugging["frequency"] + stability_margin = float("nan") + acoustic = calculate_acoustic_modes( + self.L_chamber, + D_chamber_new, + Tc, + gamma_chamber, + R_chamber, + ) + feed_stability = { + "pogo_frequency": np.nan, + "surge_frequency": np.nan, + "stability_margin": np.nan, + } # Use comprehensive stability analysis if available comprehensive_stability = None @@ -852,9 +776,11 @@ def solve_time_step( diagnostics=stability_diag, ) - # Update stability_margin from comprehensive analysis + # The comprehensive analysis owns every stability scalar in the state record. stability_margin = comprehensive_stability.get("chugging", {}).get("stability_margin", stability_margin) chugging_freq = comprehensive_stability.get("chugging", {}).get("frequency", chugging_freq) + acoustic = comprehensive_stability.get("acoustic", acoustic) + feed_stability = comprehensive_stability.get("feed_system", feed_stability) except Exception as e: import warnings warnings.warn(f"Comprehensive stability analysis failed: {e}") diff --git a/EngineDesign/frontend/src/api/client.ts b/EngineDesign/frontend/src/api/client.ts index afceb4590..b194eaffa 100644 --- a/EngineDesign/frontend/src/api/client.ts +++ b/EngineDesign/frontend/src/api/client.ts @@ -1075,6 +1075,7 @@ export interface Layer1Results { final_change?: number; best_objective?: number; best_objective_breakdown?: Record; + infeasible_reason?: string | null; /** Thrust / O-F / P_exit relative errors and RMS (dimensionless); use for “true” physics convergence. */ primary_relative_residual?: { rel_thrust?: number; diff --git a/EngineDesign/frontend/src/components/ConfigEditor.tsx b/EngineDesign/frontend/src/components/ConfigEditor.tsx index 329b23209..f2bcc3e0d 100644 --- a/EngineDesign/frontend/src/components/ConfigEditor.tsx +++ b/EngineDesign/frontend/src/components/ConfigEditor.tsx @@ -10,28 +10,29 @@ interface ConfigEditorProps { } // Section metadata for better labels and descriptions -const SECTION_META: Record = { - fluids: { label: 'Fluids', icon: '💧', description: 'Oxidizer and fuel properties' }, - injector: { label: 'Injector', icon: '🔧', description: 'Injector geometry (pintle or impinging doublet)' }, - feed_system: { label: 'Feed System', icon: '⚡', description: 'Propellant feed configuration' }, - regen_cooling: { label: 'Regenerative Cooling', icon: '❄️', description: 'Cooling channel parameters' }, - film_cooling: { label: 'Film Cooling', icon: '🌊', description: 'Film cooling settings' }, - ablative_cooling: { label: 'Ablative Cooling', icon: '🔥', description: 'Ablative material properties' }, - graphite_insert: { label: 'Graphite Insert', icon: '⬛', description: 'Throat insert configuration' }, - stainless_steel_case: { label: 'Steel Case', icon: '🔩', description: 'Case material properties' }, - discharge: { label: 'Discharge Coefficients', icon: '📊', description: 'Cd models for oxidizer/fuel' }, - spray: { label: 'Spray Modeling', icon: '💨', description: 'Atomization and spray parameters' }, - combustion: { label: 'Combustion', icon: '🔥', description: 'CEA and efficiency models' }, - chamber_geometry: { label: 'Chamber Geometry (Unified)', icon: '🎯', description: 'Unified chamber and nozzle design parameters' }, - chamber: { label: 'Chamber', icon: '🎯', description: 'Combustion chamber geometry' }, - nozzle: { label: 'Nozzle', icon: '🚀', description: 'Nozzle expansion parameters' }, - solver: { label: 'Solver', icon: '⚙️', description: 'Numerical solver settings' }, - lox_tank: { label: 'LOX Tank', icon: '🛢️', description: 'Oxidizer tank geometry' }, - fuel_tank: { label: 'Fuel Tank', icon: '⛽', description: 'Fuel tank geometry' }, - press_tank: { label: 'Pressurization Tank', icon: '🎈', description: 'Pressurant system' }, - rocket: { label: 'Rocket', icon: '🚀', description: 'Vehicle mass and geometry' }, - environment: { label: 'Environment', icon: '🌍', description: 'Launch site conditions' }, - thrust: { label: 'Thrust Profile', icon: '📈', description: 'Burn duration settings' }, +const SECTION_META: Record = { + fluids: { label: 'Fluids', description: 'Oxidizer and fuel properties' }, + injector: { label: 'Injector', description: 'Injector geometry (pintle or impinging doublet)' }, + feed_system: { label: 'Feed System', description: 'Propellant feed configuration' }, + regen_cooling: { label: 'Regenerative Cooling', description: 'Cooling channel parameters' }, + film_cooling: { label: 'Film Cooling', description: 'Film cooling settings' }, + ablative_cooling: { label: 'Ablative Cooling', description: 'Ablative material properties' }, + graphite_insert: { label: 'Graphite Insert', description: 'Throat insert configuration' }, + stainless_steel_case: { label: 'Steel Case', description: 'Case material properties' }, + discharge: { label: 'Discharge Coefficients', description: 'Cd models for oxidizer/fuel' }, + spray: { label: 'Spray Modeling', description: 'Atomization and spray parameters' }, + combustion: { label: 'Combustion', description: 'CEA and efficiency models' }, + chamber_geometry: { label: 'Chamber Geometry (Unified)', description: 'Unified chamber and nozzle design parameters' }, + chamber: { label: 'Chamber', description: 'Combustion chamber geometry' }, + nozzle: { label: 'Nozzle', description: 'Nozzle expansion parameters' }, + solver: { label: 'Solver', description: 'Numerical solver settings' }, + stability: { label: 'Stability Model', description: 'Combustion-response calibration, regulator dynamics, acoustic damping' }, + lox_tank: { label: 'LOX Tank', description: 'Oxidizer tank geometry' }, + fuel_tank: { label: 'Fuel Tank', description: 'Fuel tank geometry' }, + press_tank: { label: 'Pressurization Tank', description: 'Pressurant system' }, + rocket: { label: 'Rocket', description: 'Vehicle mass and geometry' }, + environment: { label: 'Environment', description: 'Launch site conditions' }, + thrust: { label: 'Thrust Profile', description: 'Burn duration settings' }, }; // Human-readable field labels @@ -60,6 +61,18 @@ const FIELD_LABELS: Record = { spacing: 'Element Spacing (m)', d_inlet: 'Inlet Diameter (m)', line_size: 'Feed Line Size', + // Stability model inputs (StabilityConfig) + n_interaction: 'Interaction Index n', + chi_acoustic: 'Sensitive-Lag Fraction χ', + mach_nozzle_entrance: 'Nozzle-Entrance Mach (blank = from contraction ratio)', + damping_injector_frac: 'Injector Damping Fraction', + damping_twophase_frac: 'Two-Phase Damping Fraction', + droplet_loading: 'Droplet Loading', + acoustic_gate_alpha_offset: 'Acoustic Gate Allowance (1/s)', + regulator_enabled: 'Model Dome Regulator', + regulator_corner_hz: 'Regulator Corner Frequency (Hz)', + regulator_Z_hf: 'Regulator HF Impedance (Pa·s/kg)', + regulator_max_excursion_psi: 'Regulator Excursion Bound (psi)', // Evaporation / spray-length model C_evap: 'Evaporation Calibration Constant', cp_gas: 'Combustion Gas cp (J/kg·K)', @@ -454,7 +467,6 @@ function SectionCard({ sectionKey, data, onEdit }: SectionCardProps) { const [isExpanded, setIsExpanded] = useViewState(`configSection.${sectionKey}`, false); const meta = SECTION_META[sectionKey] || { label: sectionKey.replace(/_/g, ' ').replace(/\b\w/g, c => c.toUpperCase()), - icon: '📄', description: '', }; @@ -462,7 +474,6 @@ function SectionCard({ sectionKey, data, onEdit }: SectionCardProps) { return (
- {meta.icon}

{meta.label}

Not configured

@@ -519,7 +530,6 @@ function SectionCard({ sectionKey, data, onEdit }: SectionCardProps) { className="w-full flex items-center justify-between p-4 hover:bg-[var(--color-bg-tertiary)] transition-colors" >
- {meta.icon}

{meta.label}

diff --git a/EngineDesign/frontend/src/components/ConfigurationSelector.tsx b/EngineDesign/frontend/src/components/ConfigurationSelector.tsx index 03e78dfc8..e67d7a272 100644 --- a/EngineDesign/frontend/src/components/ConfigurationSelector.tsx +++ b/EngineDesign/frontend/src/components/ConfigurationSelector.tsx @@ -14,7 +14,6 @@ import { useReadOnly } from '@stardesign-ui'; const PRETTY: Record = { pintle: 'Pintle', impinging: 'Doublet (unlike-impinging)', - coaxial: 'Coaxial', methalox: 'Methalox (LOX / CH₄)', ethalox: 'Ethalox (LOX / Ethanol)', kerolox: 'Kerolox (LOX / RP-1)', diff --git a/EngineDesign/frontend/src/components/ControllerMode.tsx b/EngineDesign/frontend/src/components/ControllerMode.tsx index b354c437f..c1157213c 100644 --- a/EngineDesign/frontend/src/components/ControllerMode.tsx +++ b/EngineDesign/frontend/src/components/ControllerMode.tsx @@ -542,7 +542,7 @@ export function ControllerMode({ config }: ControllerModeProps) {

Mass Flow Rates

- d.time || 0)).map((t, i) => ({ + d.time || 0)).map((t, i) => ({ time: t, mdot_F: results ? (results.mdot_F?.[i] || 0) : (realtimeData[i]?.mdot_F || 0), mdot_O: results ? (results.mdot_O?.[i] || 0) : (realtimeData[i]?.mdot_O || 0), @@ -847,7 +847,7 @@ export function ControllerMode({ config }: ControllerModeProps) { - + @@ -859,7 +859,7 @@ export function ControllerMode({ config }: ControllerModeProps) {

Altitude & Velocity

- d.time || 0)).map((t, i) => ({ + d.time || 0)).map((t, i) => ({ time: t, altitude: results ? results.altitude[i] : (realtimeData[i]?.altitude || 0), velocity: results ? results.velocity[i] : (realtimeData[i]?.velocity || 0), @@ -869,7 +869,7 @@ export function ControllerMode({ config }: ControllerModeProps) { - + diff --git a/EngineDesign/frontend/src/components/CustomPlotter.tsx b/EngineDesign/frontend/src/components/CustomPlotter.tsx index a3b234763..c5a1a9fb5 100644 --- a/EngineDesign/frontend/src/components/CustomPlotter.tsx +++ b/EngineDesign/frontend/src/components/CustomPlotter.tsx @@ -651,7 +651,7 @@ export function CustomPlotter({ isVisible = true }: CustomPlotterProps) {
- + } /> - + {/* Render series */} {yAxes.map((field, idx) => { diff --git a/EngineDesign/frontend/src/components/DesignRequirements.tsx b/EngineDesign/frontend/src/components/DesignRequirements.tsx index a1c0a477a..268118f4f 100644 --- a/EngineDesign/frontend/src/components/DesignRequirements.tsx +++ b/EngineDesign/frontend/src/components/DesignRequirements.tsx @@ -1,5 +1,6 @@ import { useState, useEffect, useCallback } from 'react'; import { getInjectorSchema } from '../api/client'; +import type { EngineConfig } from '../api/client'; import { useConfigChanged } from '../lib/configBus'; import { useReadOnly } from '@stardesign-ui'; import type { DesignRequirements as DesignRequirementsType, FrozenParameters } from '../api/client'; @@ -34,31 +35,47 @@ interface DesignRequirementsProps { requirements: DesignRequirementsType; onRequirementsChange: (next: DesignRequirementsType) => void; onSave: () => void; + /** Live config, for the propellant's CEA table range under the O/F target. */ + config?: EngineConfig | null; +} + +/** [lo, hi] of the loaded propellant's CEA mixture-ratio table, or null when unknown. */ +function ceaMrRange(config?: EngineConfig | null): [number, number] | null { + const cea = (config?.combustion as { cea?: { MR_range?: unknown } } | undefined)?.cea; + const r = cea?.MR_range; + if (Array.isArray(r) && r.length === 2 && Number.isFinite(Number(r[0])) && Number.isFinite(Number(r[1]))) { + return [Number(r[0]), Number(r[1])]; + } + return null; } // Metadata for every injector-specific frozen field (both families). The frozen-injector UI renders // ONLY the fields the backend says apply to the current injector type (INJECTOR_PARITY_PLAN W5), so // switching to doublet shows doublet freezes — not pintle ones. const FROZEN_INJECTOR_META: Record = { - d_pintle_tip_mm: { label: 'Pintle Tip Ø [mm]', def: 25, min: 10, max: 50, step: 1 }, - h_gap_mm: { label: 'Gap Height [mm]', def: 1.0, min: 0.2, max: 3.0, step: 0.1 }, - n_orifices: { label: '# LOX Orifices', def: 16, min: 4, max: 48, step: 2, int: true }, - d_orifice_mm: { label: 'Orifice Ø [mm]', def: 2.5, min: 0.5, max: 8, step: 0.1 }, - n_doublets: { label: '# Doublets', def: 20, min: 5, max: 40, step: 1, int: true }, - d_jet_O_mm: { label: 'LOX Jet Ø [mm]', def: 2.0, min: 0.5, max: 6, step: 0.1 }, - d_jet_F_mm: { label: 'Fuel Jet Ø [mm]', def: 2.0, min: 0.5, max: 6, step: 0.1 }, - impingement_angle_O_deg: { label: 'LOX Imp. Angle [°]', def: 50, min: 20, max: 90, step: 1 }, - impingement_angle_F_deg: { label: 'Fuel Imp. Angle [°]', def: 60, min: 20, max: 90, step: 1 }, - spacing_O_mm: { label: 'LOX Spacing [mm]', def: 6, min: 1, max: 20, step: 0.5 }, - spacing_F_mm: { label: 'Fuel Spacing [mm]', def: 6, min: 1, max: 20, step: 0.5 }, + d_pintle_tip_mm: { label: 'Pintle Tip Ø [mm]', def: 25, min: 1, max: 500, step: 1 }, + h_gap_mm: { label: 'Gap Height [mm]', def: 1.0, min: 0.05, max: 20, step: 0.1 }, + n_orifices: { label: '# LOX Orifices', def: 16, min: 1, max: 400, step: 1, int: true }, + d_orifice_mm: { label: 'Orifice Ø [mm]', def: 2.5, min: 0.1, max: 30, step: 0.1 }, + n_doublets: { label: '# Doublets', def: 20, min: 1, max: 400, step: 1, int: true }, + d_jet_O_mm: { label: 'LOX Jet Ø [mm]', def: 2.0, min: 0.1, max: 30, step: 0.1 }, + d_jet_F_mm: { label: 'Fuel Jet Ø [mm]', def: 2.0, min: 0.1, max: 30, step: 0.1 }, + impingement_angle_O_deg: { label: 'LOX Imp. Angle [°]', def: 50, min: 1, max: 90, step: 1 }, + impingement_angle_F_deg: { label: 'Fuel Imp. Angle [°]', def: 60, min: 1, max: 90, step: 1 }, + spacing_O_mm: { label: 'LOX Spacing [mm]', def: 6, min: 0.5, max: 200, step: 0.5 }, + spacing_F_mm: { label: 'Fuel Spacing [mm]', def: 6, min: 0.5, max: 200, step: 0.5 }, }; export function DesignRequirements({ requirements, onRequirementsChange, onSave, + config, }: DesignRequirementsProps) { const readOnly = useReadOnly(); + const mrRange = ceaMrRange(config); + const ofOutsideTable = mrRange !== null && Number.isFinite(requirements.optimal_of_ratio) + && (requirements.optimal_of_ratio < mrRange[0] || requirements.optimal_of_ratio > mrRange[1]); // Which injector-specific frozen fields to show — fetched from the backend authority so the UI // matches the live injector type (no hardcoded pintle assumptions). const [injectorFrozenFields, setInjectorFrozenFields] = useState([]); @@ -126,7 +143,7 @@ export function DesignRequirements({ {/* Performance Targets */}
-

🎯 Performance Targets

+

Performance Targets

@@ -206,9 +225,8 @@ export function DesignRequirements({ value={requirements.target_burn_time} onChange={(e) => updateField('target_burn_time', parseFloat(e.target.value))} className="w-full px-3 py-2 bg-[var(--color-bg-primary)] border border-[var(--color-border)] rounded-lg text-[var(--color-text-primary)] focus:outline-none focus:ring-2 focus:ring-blue-500" - min="1" - max="60" - step="1" + min="0.1" + step="0.5" />

Design burn time. Flight sim will truncate if propellant depletes earlier.

@@ -217,7 +235,7 @@ export function DesignRequirements({ {/* Tank Pressures */}
-

🔋 Tank Pressures

+

Tank Pressures

Maximum operating pressure in LOX tank.

@@ -278,8 +295,7 @@ export function DesignRequirements({ ? 'bg-[var(--color-bg-secondary)] border-[var(--color-border)] text-[var(--color-text-secondary)] cursor-not-allowed opacity-50' : 'bg-[var(--color-bg-primary)] border-[var(--color-border)] text-[var(--color-text-primary)]' }`} - min="100" - max="5000" + min="1" step="25" />

@@ -291,7 +307,7 @@ export function DesignRequirements({ {/* Geometry Constraints */}

-

📏 Geometry Constraints

+

Geometry Constraints

{/* Stock-size snapping. Ablative sleeve / chamber tube / case come in fixed sizes, so a continuous optimum like 4.2" is not purchasable. Snapping INSIDE the search means the returned design is already buildable, instead of being rounded afterwards (which moves @@ -388,8 +404,7 @@ export function DesignRequirements({ value={requirements.max_engine_length} onChange={(e) => updateField('max_engine_length', parseFloat(e.target.value))} className="w-full px-3 py-2 bg-[var(--color-bg-primary)] border border-[var(--color-border)] rounded-lg text-[var(--color-text-primary)] focus:outline-none focus:ring-2 focus:ring-blue-500" - min="0.1" - max="3.0" + min="0.01" step="0.05" />

Maximum total engine length (chamber + nozzle). Must fit in vehicle.

@@ -404,8 +419,7 @@ export function DesignRequirements({ value={requirements.max_chamber_outer_diameter} onChange={(e) => updateField('max_chamber_outer_diameter', parseFloat(e.target.value))} className="w-full px-3 py-2 bg-[var(--color-bg-primary)] border border-[var(--color-border)] rounded-lg text-[var(--color-text-primary)] focus:outline-none focus:ring-2 focus:ring-blue-500" - min="0.05" - max="1.0" + min="0.01" step="0.01" />

Maximum chamber outer diameter (including wall thickness and cooling jacket).

@@ -420,8 +434,7 @@ export function DesignRequirements({ value={requirements.max_nozzle_exit_diameter} onChange={(e) => updateField('max_nozzle_exit_diameter', parseFloat(e.target.value))} className="w-full px-3 py-2 bg-[var(--color-bg-primary)] border border-[var(--color-border)] rounded-lg text-[var(--color-text-primary)] focus:outline-none focus:ring-2 focus:ring-blue-500" - min="0.05" - max="1.0" + min="0.01" step="0.01" />

Maximum nozzle exit outer diameter. Constrains expansion ratio.

@@ -431,7 +444,7 @@ export function DesignRequirements({ {/* L* Constraints */}
-

📐 L* (Characteristic Length) Constraints

+

L* (Characteristic Length) Constraints

@@ -458,9 +470,8 @@ export function DesignRequirements({ value={requirements.max_Lstar} onChange={(e) => updateField('max_Lstar', parseFloat(e.target.value))} className="w-full px-3 py-2 bg-[var(--color-bg-primary)] border border-[var(--color-border)] rounded-lg text-[var(--color-text-primary)] focus:outline-none focus:ring-2 focus:ring-blue-500" - min="0.5" - max="3.0" - step="0.1" + min="0.05" + step="0.05" />

Maximum characteristic length. Higher = better combustion but heavier/longer chamber. Typical: 1.0-2.0m for LOX/hydrocarbon.

@@ -469,16 +480,13 @@ export function DesignRequirements({ {/* Stability Requirements */}
-

🛡️ Stability Requirements

+

Stability Requirements

- New Comprehensive Stability Analysis:
- • Uses stability_score (0-1) and stability_state ("stable"/"marginal"/"unstable")
- • Considers chugging, acoustic modes, feed system, and mode coupling
- • Stable: score ≥ 0.75 (recommended for flight)
- • Marginal: 0.4 ≤ score < 0.75 (acceptable with caution)
- • Unstable: score < 0.4 (not acceptable) + Score maps the limiting gate margin (chug gain margin, worst acoustic mode) onto 0–1. + A design is stable when every margin is at least 1.05 and no mode is driven, + marginal above 0.95, unstable below.

@@ -496,7 +504,7 @@ export function DesignRequirements({ max="1.0" step="0.05" /> -

Minimum stability score (0-1). 0.75 = 'stable', 0.4 = 'marginal', <0.4 = 'unstable'

+

Minimum stability score (0–1) the optimizer must reach; the score is 0.44 at the marginal boundary and 1.0 when every margin is 1.3 or better.

@@ -589,7 +597,7 @@ export function DesignRequirements({
- 🔒 Frozen Parameters (Optional) + Frozen Parameters (Optional)
@@ -845,7 +853,7 @@ export function DesignRequirements({ {/* Summary */}
-

📋 Design Summary

+

Design Summary

Target Thrust

diff --git a/EngineDesign/frontend/src/components/FlightSimulation.tsx b/EngineDesign/frontend/src/components/FlightSimulation.tsx index d3367bfe9..d9874ef79 100644 --- a/EngineDesign/frontend/src/components/FlightSimulation.tsx +++ b/EngineDesign/frontend/src/components/FlightSimulation.tsx @@ -37,6 +37,8 @@ interface FlightSimulationProps { config: EngineConfig | null; isVisible?: boolean; onConfigUpdated?: (config: EngineConfig) => void; + /** Design-state slice this mount persists under (the Layer 4 mount uses its own). */ + sliceKey?: string; } // Session storage handled via timeseriesSession utility (shared with TimeSeriesMode) @@ -167,7 +169,7 @@ function MetricCard({ ); } -export function FlightSimulation({ config, isVisible = true, onConfigUpdated }: FlightSimulationProps) { +export function FlightSimulation({ config, isVisible = true, onConfigUpdated, sliceKey = 'flight' }: FlightSimulationProps) { // RocketPy availability const [rocketPyAvailable, setRocketPyAvailable] = useState(null); const [rocketPyMessage, setRocketPyMessage] = useState(''); @@ -220,7 +222,9 @@ export function FlightSimulation({ config, isVisible = true, onConfigUpdated }: // The fields Save Configuration does NOT write. Everything else on this tab // round-trips through config.rocket / config.environment on an explicit save; // these seven have no config field at all, so they were lost on reload. - useDesignSlice('flight', { + // Keyed per mount point: the Flight tab and the optimizer's Layer 4 both render this + // component, and two registrations under one key would leave one of them un-restored. + useDesignSlice(sliceKey, { atmosphereModel: [atmosphereModel, setAtmosphereModel], autoInertia: [autoInertia, setAutoInertia], noseFineness: [noseFineness, setNoseFineness], @@ -1538,7 +1542,7 @@ export function FlightSimulation({ config, isVisible = true, onConfigUpdated }:

Altitude vs Time

- + Velocity vs Time
- +
- +
- + `Time: ${Number(label).toFixed(3)}s`} /> (value === 'lox_pressure' ? 'LOX Tank' : 'Fuel Tank')} />
- + `Time: ${Number(label).toFixed(3)}s`} /> (value === 'lox_fill' ? 'LOX Tank' : 'Fuel Tank')} /> diff --git a/EngineDesign/frontend/src/components/HeatFluxProfileChart.tsx b/EngineDesign/frontend/src/components/HeatFluxProfileChart.tsx index bd48b2a3f..3260850df 100644 --- a/EngineDesign/frontend/src/components/HeatFluxProfileChart.tsx +++ b/EngineDesign/frontend/src/components/HeatFluxProfileChart.tsx @@ -350,7 +350,8 @@ export function HeatFluxProfileChart({ data }: HeatFluxProfileChartProps) { /> { const match = value.match(/t_(\d+)/); if (match && time) { @@ -411,7 +412,8 @@ export function HeatFluxProfileChart({ data }: HeatFluxProfileChartProps) { /> { const match = value.match(/t_(\d+)/); if (match && time) { @@ -495,7 +497,7 @@ export function HeatFluxProfileChart({ data }: HeatFluxProfileChartProps) { /> )} - + -

⚙️ Optimization Settings

+

Optimization Settings

{/* When the throat is solved from the thrust target the engine lands on that target by construction, so a tolerance has nothing to authorise -- showing an @@ -1292,14 +1292,14 @@ export function Layer1Optimization({ : 'bg-blue-600 hover:bg-blue-700 hover:scale-105' }`} > - {isRunning ? '🔄 Running Optimization...' : '🚀 Run Layer 1 Optimization'} + {isRunning ? 'Running Optimization...' : 'Run Layer 1 Optimization'} {isRunning && ( )}
@@ -1307,7 +1307,7 @@ export function Layer1Optimization({ {/* Progress */} {(isRunning || progress > 0) && (
-

📊 Progress

+

Progress

{/* Progress Bar */}
@@ -1336,7 +1336,9 @@ export function Layer1Optimization({ {objectiveHistory.length > 0 ? ( - + {/* Legend on top and a real bottom margin: with both the legend and the axis + title in the strip under the axis they drew over each other mid-run. */} + (Math.abs(v) >= 1e4 || (v !== 0 && Math.abs(v) < 1e-2) ? v.toExponential(0) : v.toLocaleString())} label={{ value: 'Weighted penalty sum (log)', angle: -90, position: 'insideLeft', fill: 'var(--color-text-secondary)' }} /> -

❌ Error: {error}

+

Error: {error}

)} @@ -1454,7 +1460,7 @@ export function Layer1Optimization({ {results && results.performance && (
-

✅ Optimization Results

+

Optimization Results

{results.config_yaml && ( )}
{/* Key Performance Metrics */}
-

🎯 Performance

+

Performance

-

🧠 Objective Diagnostics

+

Objective Diagnostics

+ {results.convergence_info.infeasible_reason && ( +

+ {results.convergence_info.infeasible_reason} +

+ )}
-

🔋 Optimized Tank Pressures

+

Optimized Tank Pressures

-

🫧 Injector Pressure Drops

+

Injector Pressure Drops

{isRunning && ( )}
@@ -630,8 +630,8 @@ export function Layer2Optimization({ />
{message &&

{message}

} - {error &&

❌ {error}

} - {successMessage &&

✅ {successMessage}

} + {error &&

{error}

} + {successMessage &&

{successMessage}

}
)} @@ -640,7 +640,7 @@ export function Layer2Optimization({ {/* Convergence History */}

- 📈 Convergence History + Convergence History

{objectiveHistory.length > 0 ? ( @@ -668,7 +668,7 @@ export function Layer2Optimization({ {/* Pressure Curves */}

- 🌊 Pressure Curves (Current Best) + Pressure Curves (Current Best) {pressureCurves.some(p => p.copv !== undefined) && ( + COPV )} @@ -732,19 +732,19 @@ export function Layer2Optimization({ {results && (
-

✨ Final Optimization Results

+

Final Optimization Results

@@ -765,7 +765,7 @@ export function Layer2Optimization({ results.summary.thrust_curve_time.length > 0 && results.summary.thrust_curve_values.length > 0 && (

- 🚀 Thrust Curve (Time Series, No Ablation/Oxidation) + Thrust Curve (Time Series, No Ablation/Oxidation)

@@ -812,7 +812,7 @@ export function Layer2Optimization({ results.summary.thrust_curve_time.length > 0 && results.summary.of_curve_values.length > 0 && (

- ⚗️ O/F Ratio (Mixture Ratio) Curve (Time Series, No Ablation/Oxidation) + O/F Ratio (Mixture Ratio) Curve (Time Series, No Ablation/Oxidation)

@@ -859,7 +859,7 @@ export function Layer2Optimization({ results.summary.thrust_curve_time.length > 0 && results.summary.delta_p_inj_O_psi.length > 0 && results.summary.delta_p_inj_F_psi.length > 0 && (

- 💧 Injector Pressure Drops (Time Series, No Ablation/Oxidation) + Injector Pressure Drops (Time Series, No Ablation/Oxidation)

@@ -916,7 +916,7 @@ export function Layer2Optimization({ results.summary.copv_time_s.length > 0 && results.summary.copv_pressure_trace_Pa.length > 0 && (

- 🔋 COPV & Tank Pressures + COPV & Tank Pressures

- 🎮 Controller Simulation + Controller Simulation

{controllerLoading && ( @@ -1036,7 +1036,7 @@ export function Layer2Optimization({ onClick={handleStopController} className="px-4 py-2 rounded-lg font-medium transition-colors bg-red-600 hover:bg-red-700 text-white" > - ⏹ Stop + Stop ) : ( <> @@ -1052,7 +1052,7 @@ export function Layer2Optimization({ ▶ Run from Layer 2 Results
)} diff --git a/EngineDesign/frontend/src/components/Layer3Optimization.tsx b/EngineDesign/frontend/src/components/Layer3Optimization.tsx index 38d8c55f3..7d62d30c2 100644 --- a/EngineDesign/frontend/src/components/Layer3Optimization.tsx +++ b/EngineDesign/frontend/src/components/Layer3Optimization.tsx @@ -387,9 +387,9 @@ export function Layer3Optimization({ disabled={isRunning || readOnly} className="px-3 py-2 bg-[var(--color-bg-primary)] border border-[var(--color-border)] rounded-lg text-[var(--color-text-primary)] text-sm focus:border-orange-500 focus:outline-none" > - - - + + + {settings.optimization_method === 'gradient' && 'Exploits monotonic thickness-recession relationship'} @@ -430,7 +430,7 @@ export function Layer3Optimization({ : 'bg-orange-600 hover:bg-orange-700 shadow-lg shadow-orange-500/20' }`} > - {isRunning ? '🔄 Optimizing...' : '🔥 Run Layer 3'} + {isRunning ? 'Optimizing...' : 'Run Layer 3'} {isRunning && ( )}
@@ -464,8 +464,8 @@ export function Layer3Optimization({ />
{message &&

{message}

} - {error &&

❌ {error}

} - {successMessage &&

✅ {successMessage}

} + {error &&

{error}

} + {successMessage &&

{successMessage}

}
)} @@ -473,19 +473,19 @@ export function Layer3Optimization({ {results && (
-

🔥 Final Optimization Results

+

Final Optimization Results

@@ -554,7 +554,7 @@ export function Layer3Optimization({ {/* Convergence History */}

- 📈 Convergence History + Convergence History

{objectiveHistory.length > 0 ? ( @@ -582,7 +582,7 @@ export function Layer3Optimization({ {/* Pressure Curves (Real-time or Final) */}

- 🌊 Pressure Curves {isRunning ? '(Baseline)' : '(Optimized)'} + Pressure Curves {isRunning ? '(Baseline)' : '(Optimized)'}

{pressureCurves.length > 0 ? ( @@ -615,7 +615,7 @@ export function Layer3Optimization({

- 🚀 Thrust Curve + Thrust Curve

@@ -638,7 +638,7 @@ export function Layer3Optimization({

- 🔥 Chamber Pressure + Chamber Pressure

@@ -664,7 +664,7 @@ export function Layer3Optimization({

- ⚗️ Mixture Ratio (O/F) + Mixture Ratio (O/F)

@@ -687,7 +687,7 @@ export function Layer3Optimization({

- 🛡️ Cumulative Recession + Cumulative Recession

@@ -716,7 +716,7 @@ export function Layer3Optimization({

- Recession Rates + Recession Rates

@@ -742,7 +742,7 @@ export function Layer3Optimization({

- 📏 Diameters & L* + Diameters & L*

@@ -774,7 +774,7 @@ export function Layer3Optimization({

- 🧪 Throat Recession Rate Breakdown + Throat Recession Rate Breakdown

Diagnostic split of throat recession rate into oxidation and thermal ablation. @@ -808,7 +808,7 @@ export function Layer3Optimization({

- 🔳 Areas & Contraction Ratio + Areas & Contraction Ratio

diff --git a/EngineDesign/frontend/src/components/Layer4Optimization.tsx b/EngineDesign/frontend/src/components/Layer4Optimization.tsx deleted file mode 100644 index f11028096..000000000 --- a/EngineDesign/frontend/src/components/Layer4Optimization.tsx +++ /dev/null @@ -1,1093 +0,0 @@ -import { useState, useEffect, useMemo, useCallback } from 'react'; -import { - LineChart, - Line, - XAxis, - YAxis, - CartesianGrid, - Tooltip, - Legend, - ResponsiveContainer, - ReferenceLine, -} from 'recharts'; -import { - runFlightSimulation, - checkRocketPy, - type FlightSimRequest, - type FlightSimResponse, - type FlightEnvironmentConfig, - type FlightRocketConfig, - type FlightTankConfig, - type TimeSeriesData, - type DesignRequirements, - type SaveDesignRequirementsResponse, -} from '../api/client'; -import { useDesignSlice } from '../lib/designState'; -import { useReadOnly } from '@stardesign-ui'; -import { useViewState } from '../lib/viewState'; - -interface Layer4OptimizationProps { - requirements: DesignRequirements; - isDirty: boolean; - saveRequirementsToServer: ( - reqs: DesignRequirements - ) => Promise<{ error?: string; data?: SaveDesignRequirementsResponse }>; -} - -// Session storage key -const TIMESERIES_RESULTS_KEY = 'timeseries_results'; - -interface StoredTimeSeriesResults { - data: TimeSeriesData; - timestamp: number; -} - -function loadTimeSeriesFromSession(): TimeSeriesData | null { - try { - const stored = sessionStorage.getItem(TIMESERIES_RESULTS_KEY); - if (!stored) return null; - const parsed: StoredTimeSeriesResults = JSON.parse(stored); - return parsed.data; - } catch { - return null; - } -} - -// Fluid densities (kg/m³) -const LOX_DENSITY = 1141; // Liquid oxygen at boiling point -const RP1_DENSITY = 820; // RP-1 kerosene - -// Fill factor - conservative to account for RocketPy's internal density calculations -const FILL_FACTOR = 0.85; // 85% fill factor - -// Calculate cylindrical tank volume (m³) -function calculateTankVolume(height: number, radius: number): number { - return Math.PI * radius * radius * height; -} - -// Calculate max propellant mass for a tank (kg) -function calculateMaxPropellantMass(tank: FlightTankConfig, density: number): number { - const volume = calculateTankVolume(tank.height, tank.radius); - return volume * density * FILL_FACTOR; -} - -// Generate time-series arrays from manual parameters -function generateThrustCurve( - thrust_N: number, - burn_time_s: number, - mdot_O_kg_s: number, - mdot_F_kg_s: number, - n_points: number = 100 -): { - time: number[]; - thrust: number[]; - mdot_O: number[]; - mdot_F: number[]; -} { - const time: number[] = []; - const thrust: number[] = []; - const mdot_O: number[] = []; - const mdot_F: number[] = []; - - for (let i = 0; i < n_points; i++) { - const t = (i / (n_points - 1)) * burn_time_s; - time.push(t); - // Constant thrust profile (can be modified for more complex profiles) - thrust.push(thrust_N); - mdot_O.push(mdot_O_kg_s); - mdot_F.push(mdot_F_kg_s); - } - - return { time, thrust, mdot_O, mdot_F }; -} - -// Helper component for result cards -function ResultCard({ - label, - value, - unit, - decimals = 2, - color = 'cyan', -}: { - label: string; - value: number | string | undefined; - unit?: string; - decimals?: number; - color?: string; -}) { - const colorClasses: Record = { - cyan: 'bg-cyan-500/10 border-cyan-500/30', - green: 'bg-green-500/10 border-green-500/30', - blue: 'bg-blue-500/10 border-blue-500/30', - purple: 'bg-purple-500/10 border-purple-500/30', - orange: 'bg-orange-500/10 border-orange-500/30', - red: 'bg-red-500/10 border-red-500/30', - yellow: 'bg-yellow-500/10 border-yellow-500/30', - }; - - const textColorClasses: Record = { - cyan: 'text-cyan-400', - green: 'text-green-400', - blue: 'text-blue-400', - purple: 'text-purple-400', - orange: 'text-orange-400', - red: 'text-red-400', - yellow: 'text-yellow-400', - }; - - const displayValue = typeof value === 'number' - ? value.toFixed(decimals) - : value !== undefined && value !== null - ? String(value) - : '-'; - - return ( -
-

{label}

-

- {displayValue} - {unit && {unit}} -

-
- ); -} - -// Collapsible section component -function CollapsibleSection({ - title, - icon, - children, - defaultExpanded = false, -}: { - title: string; - icon: React.ReactNode; - children: React.ReactNode; - defaultExpanded?: boolean; -}) { - const [isExpanded, setIsExpanded] = useViewState(`layer4Section.${title}`, defaultExpanded); - - return ( -
- - {isExpanded &&
{children}
} -
- ); -} - -// Input field component -function InputField({ - label, - value, - onChange, - unit, - help, - min, - max, - step, - disabled, -}: { - label: string; - value: string; - onChange: (value: string) => void; - unit?: string; - help?: string; - min?: number; - max?: number; - step?: number; - disabled?: boolean; -}) { - const readOnly = useReadOnly(); - // Every numeric field on this tab renders through here, and they all edit - // the vehicle definition this tab stores in the design. - return ( -
- - onChange(e.target.value)} - min={min} - max={max} - step={step} - disabled={disabled || readOnly} - className="w-full px-3 py-2 rounded-lg border border-[var(--color-border)] bg-[var(--color-bg-tertiary)] text-[var(--color-text-primary)] focus:ring-2 focus:ring-cyan-500/50 focus:border-cyan-500 disabled:opacity-50" - /> - {help &&

{help}

} -
- ); -} - -// Helper to get tomorrow's date as tuple (for atmospheric data availability) -function getTomorrowDateTuple(): [number, number, number, number] { - const tomorrow = new Date(); - tomorrow.setDate(tomorrow.getDate() + 1); - return [tomorrow.getFullYear(), tomorrow.getMonth() + 1, tomorrow.getDate(), 12]; -} - -type DataSource = 'manual' | 'timeseries'; - -export function Layer4Optimization({ - requirements, - isDirty, - saveRequirementsToServer, -}: Layer4OptimizationProps) { - // RocketPy availability - const [rocketPyAvailable, setRocketPyAvailable] = useState(null); - const [rocketPyMessage, setRocketPyMessage] = useState(''); - - // Data source selection - const readOnly = useReadOnly(); - const [dataSource, setDataSource] = useState('manual'); - - // Time series data from session (optional) - const [timeSeriesData, setTimeSeriesData] = useState(null); - - // Manual engine parameters - const [thrust, setThrust] = useState(7000); // N - const [burnTime, setBurnTime] = useState(10); // s - const [mdotO, setMdotO] = useState(2.0); // kg/s - const [mdotF, setMdotF] = useState(0.87); // kg/s (for O/F ~2.3) - - // Environment config - const [envConfig, setEnvConfig] = useState({ - latitude: 32.99, - longitude: -106.97, - elevation: 1401, - date: getTomorrowDateTuple(), - }); - - // Rocket config - const [rocketConfig, setRocketConfig] = useState({ - airframe_mass: 50, - engine_mass: 10, - lox_tank_structure_mass: 3, - fuel_tank_structure_mass: 2, - radius: 0.1, - rocket_length: 3.0, - motor_position: 0.5, - inertia: [10, 10, 0.5], - }); - - // Propellant masses - const [loxMass, setLoxMass] = useState(20); - const [fuelMass, setFuelMass] = useState(10); - - // Tank config - const [loxTankConfig, setLoxTankConfig] = useState({ - mass: 20, - height: 0.5, - radius: 0.08, - position: 1.5, - }); - - const [fuelTankConfig, setFuelTankConfig] = useState({ - mass: 10, - height: 0.3, - radius: 0.08, - position: 0.8, - }); - - // Part of the design. This tab duplicates config.rocket / config.environment - // and syncs with neither, so its vehicle, launch site and tank geometry only - // ever lived in React state. Registering it is what makes them persist; the - // duplication itself is left alone deliberately -- reconciling it with the - // config is a behaviour change, not a persistence fix. - useDesignSlice('layer4', { - dataSource: [dataSource, setDataSource], - thrust: [thrust, setThrust], - burnTime: [burnTime, setBurnTime], - mdotO: [mdotO, setMdotO], - mdotF: [mdotF, setMdotF], - envConfig: [envConfig, setEnvConfig], - rocketConfig: [rocketConfig, setRocketConfig], - loxMass: [loxMass, setLoxMass], - fuelMass: [fuelMass, setFuelMass], - loxTankConfig: [loxTankConfig, setLoxTankConfig], - fuelTankConfig: [fuelTankConfig, setFuelTankConfig], - }); - - - // Results - const [results, setResults] = useState(null); - const [isRunning, setIsRunning] = useState(false); - const [error, setError] = useState(null); - - // Check RocketPy availability on mount - useEffect(() => { - checkRocketPy().then((response) => { - setRocketPyAvailable(response.data?.available ?? false); - setRocketPyMessage(response.data?.message ?? ''); - }); - }, []); - - // Load time series data from session - useEffect(() => { - const data = loadTimeSeriesFromSession(); - setTimeSeriesData(data); - }, []); - - // Update from requirements when available - useEffect(() => { - if (requirements) { - // Update thrust and burn time from requirements - if (requirements.target_thrust) { - setThrust(requirements.target_thrust); - } - if (requirements.target_burn_time) { - setBurnTime(requirements.target_burn_time); - } - // Calculate mdot from thrust and O/F ratio - if (requirements.target_thrust && requirements.optimal_of_ratio) { - // Approximate: F = mdot_total * Ve, assume Ve ~ 2500 m/s for LOX/RP-1 - const Ve_approx = 2500; - const mdot_total = requirements.target_thrust / Ve_approx; - const OF = requirements.optimal_of_ratio; - setMdotO(mdot_total * OF / (1 + OF)); - setMdotF(mdot_total / (1 + OF)); - } - // Update tank masses from requirements if available - if (requirements.lox_tank_capacity_kg) { - setLoxMass(requirements.lox_tank_capacity_kg); - setLoxTankConfig(prev => ({ ...prev, mass: requirements.lox_tank_capacity_kg! })); - } - if (requirements.fuel_tank_capacity_kg) { - setFuelMass(requirements.fuel_tank_capacity_kg); - setFuelTankConfig(prev => ({ ...prev, mass: requirements.fuel_tank_capacity_kg! })); - } - } - }, [requirements]); - - // Check if we have time-series data - const hasTimeSeriesData = timeSeriesData !== null && - timeSeriesData.time && timeSeriesData.time.length > 0 && - timeSeriesData.thrust_kN && timeSeriesData.thrust_kN.length > 0; - - // Calculate tank capacities based on geometry - const loxTankCapacity = useMemo(() => ({ - volume: calculateTankVolume(loxTankConfig.height, loxTankConfig.radius), - maxMass: calculateMaxPropellantMass(loxTankConfig, LOX_DENSITY), - }), [loxTankConfig.height, loxTankConfig.radius]); - - const fuelTankCapacity = useMemo(() => ({ - volume: calculateTankVolume(fuelTankConfig.height, fuelTankConfig.radius), - maxMass: calculateMaxPropellantMass(fuelTankConfig, RP1_DENSITY), - }), [fuelTankConfig.height, fuelTankConfig.radius]); - - // Check if propellant masses exceed tank capacity - const loxOverfilled = loxMass > loxTankCapacity.maxMass; - const fuelOverfilled = fuelMass > fuelTankCapacity.maxMass; - - // Track mass adjustments made during simulation - const [massAdjustments, setMassAdjustments] = useState<{ - loxOriginal?: number; - loxAdjusted?: number; - fuelOriginal?: number; - fuelAdjusted?: number; - } | null>(null); - - // Run flight simulation - const runSimulation = useCallback(async () => { - if (isDirty) { - const saveResp = await saveRequirementsToServer(requirements); - if (saveResp.error) { - setError(saveResp.error); - return; - } - } - - setIsRunning(true); - setError(null); - setResults(null); - setMassAdjustments(null); - - try { - let time_array: number[]; - let thrust_array: number[]; - let mdot_O_array: number[]; - let mdot_F_array: number[]; - - if (dataSource === 'timeseries' && hasTimeSeriesData && timeSeriesData) { - // Use time-series data - time_array = timeSeriesData.time; - thrust_array = timeSeriesData.thrust_kN.map(t => t * 1000); // kN to N - mdot_O_array = timeSeriesData.mdot_O_kg_s || timeSeriesData.time.map(() => mdotO); - mdot_F_array = timeSeriesData.mdot_F_kg_s || timeSeriesData.time.map(() => mdotF); - } else { - // Generate from manual parameters - const generated = generateThrustCurve(thrust, burnTime, mdotO, mdotF, 100); - time_array = generated.time; - thrust_array = generated.thrust; - mdot_O_array = generated.mdot_O; - mdot_F_array = generated.mdot_F; - } - - // Auto-cap propellant masses to tank capacity to prevent overfill errors - let effectiveLoxMass = loxMass; - let effectiveFuelMass = fuelMass; - const adjustments: typeof massAdjustments = {}; - - if (loxMass > loxTankCapacity.maxMass) { - adjustments.loxOriginal = loxMass; - adjustments.loxAdjusted = loxTankCapacity.maxMass; - effectiveLoxMass = loxTankCapacity.maxMass; - } - - if (fuelMass > fuelTankCapacity.maxMass) { - adjustments.fuelOriginal = fuelMass; - adjustments.fuelAdjusted = fuelTankCapacity.maxMass; - effectiveFuelMass = fuelTankCapacity.maxMass; - } - - if (Object.keys(adjustments).length > 0) { - setMassAdjustments(adjustments); - } - - // Build request with capped masses - const request: FlightSimRequest = { - time_array, - thrust_array, - mdot_O_array, - mdot_F_array, - lox_mass_kg: effectiveLoxMass, - fuel_mass_kg: effectiveFuelMass, - lox_tank: loxTankConfig, - fuel_tank: fuelTankConfig, - environment: envConfig, - rocket: rocketConfig, - }; - - const response = await runFlightSimulation(request); - - if (response.error) { - setError(response.error); - } else if (response.data) { - setResults(response.data); - } - } catch (err) { - setError(err instanceof Error ? err.message : 'Unknown error'); - } finally { - setIsRunning(false); - } - }, [ - dataSource, - hasTimeSeriesData, - timeSeriesData, - thrust, - burnTime, - mdotO, - mdotF, - loxMass, - fuelMass, - loxTankConfig, - fuelTankConfig, - envConfig, - rocketConfig, - loxTankCapacity.maxMass, - fuelTankCapacity.maxMass, - isDirty, - requirements, - saveRequirementsToServer, - ]); - - // Prepare chart data - const altitudeChartData = useMemo(() => { - if (!results?.trajectory?.time || !results?.trajectory?.altitude) return []; - return results.trajectory.time.map((t, i) => ({ - time: t, - altitude: results.trajectory!.altitude[i], - })); - }, [results]); - - const velocityChartData = useMemo(() => { - if (!results?.trajectory?.time || !results?.trajectory?.velocity) return []; - return results.trajectory.time.map((t, i) => ({ - time: t, - velocity: results.trajectory!.velocity[i], - })); - }, [results]); - - const targetApogee = requirements?.target_apogee ?? 3048; - - // Helper to format date tuple for display - const formatDateForInput = (dateTuple: [number, number, number, number]): string => { - const [year, month, day] = dateTuple; - return `${year}-${String(month).padStart(2, '0')}-${String(day).padStart(2, '0')}`; - }; - - // Helper to parse date input to tuple - const parseDateInput = (dateStr: string): [number, number, number, number] => { - const [year, month, day] = dateStr.split('-').map(Number); - return [year, month, day, 12]; - }; - - // Calculate derived values for display - const totalMdot = mdotO + mdotF; - const ofRatio = mdotF > 0 ? mdotO / mdotF : 0; - const estimatedIsp = thrust / (totalMdot * 9.81); - - return ( -
- {/* Header */} -
-

✈️ Layer 4: Flight Simulation

-

- Run trajectory simulation using RocketPy to validate apogee and flight performance. - Configure engine parameters directly or use data from Time-Series Analysis. -

-
- - {/* RocketPy Status */} - {rocketPyAvailable === false && ( -
-

❌ RocketPy is not available

-

- {rocketPyMessage || 'Install RocketPy to enable flight simulation:'} pip install rocketpy -

-
- )} - - {rocketPyAvailable === true && ( -
-

✅ RocketPy is available

-
- )} - - {/* Data Source Selection */} -
-

📊 Data Source

- -
- - -
-
- - {/* Configuration Sections */} -
- {/* Engine Parameters (only shown for manual mode) */} - {dataSource === 'manual' && ( - 🔥} defaultExpanded={true}> -
- setThrust(parseFloat(v) || 0)} - unit="N" - help="Engine thrust (constant profile)" - /> - setBurnTime(parseFloat(v) || 0)} - unit="s" - help="Total burn duration" - /> - setMdotO(parseFloat(v) || 0)} - unit="kg/s" - help="Oxidizer mass flow rate" - /> - setMdotF(parseFloat(v) || 0)} - unit="kg/s" - help="Fuel mass flow rate" - /> -
- - {/* Derived values */} -
-
-

O/F Ratio

-

{ofRatio.toFixed(2)}

-
-
-

Total mdot

-

{totalMdot.toFixed(2)} kg/s

-
-
-

Est. Isp

-

{estimatedIsp.toFixed(0)} s

-
-
-
- )} - - {/* Propellant Masses */} - ⛽} defaultExpanded={true}> -
-
- setLoxMass(parseFloat(v) || 0)} - unit="kg" - help="Liquid oxygen propellant mass" - /> - {/* Tank capacity info */} -
-

- Tank: {(loxTankCapacity.volume * 1000).toFixed(1)}L ({loxTankCapacity.volume.toFixed(4)} m³) -

-

- Max capacity: {loxTankCapacity.maxMass.toFixed(1)} kg (85% fill) -

- {loxOverfilled && ( -

- ⚠️ Overfilled by {(loxMass - loxTankCapacity.maxMass).toFixed(1)} kg - will be auto-capped -

- )} -
-
-
- setFuelMass(parseFloat(v) || 0)} - unit="kg" - help="Fuel propellant mass" - /> - {/* Tank capacity info */} -
-

- Tank: {(fuelTankCapacity.volume * 1000).toFixed(1)}L ({fuelTankCapacity.volume.toFixed(4)} m³) -

-

- Max capacity: {fuelTankCapacity.maxMass.toFixed(1)} kg (85% fill) -

- {fuelOverfilled && ( -

- ⚠️ Overfilled by {(fuelMass - fuelTankCapacity.maxMass).toFixed(1)} kg - will be auto-capped -

- )} -
-
-
- - {/* Propellant consumption estimate */} - {dataSource === 'manual' && ( -
-

- At current flow rates, burn will consume{' '} - {(mdotO * burnTime).toFixed(1)} kg LOX and{' '} - {(mdotF * burnTime).toFixed(1)} kg fuel. - {mdotO * burnTime > loxMass && ⚠️ LOX will run out!} - {mdotF * burnTime > fuelMass && ⚠️ Fuel will run out!} -

-
- )} -
- - {/* Environment */} - 🌍} defaultExpanded={false}> -
- setEnvConfig(prev => ({ ...prev, latitude: parseFloat(v) || 0 }))} - unit="°" - /> - setEnvConfig(prev => ({ ...prev, longitude: parseFloat(v) || 0 }))} - unit="°" - /> - setEnvConfig(prev => ({ ...prev, elevation: parseFloat(v) || 0 }))} - unit="m" - /> -
- - setEnvConfig(prev => ({ ...prev, date: parseDateInput(e.target.value) }))} - className="w-full px-3 py-2 rounded-lg border border-[var(--color-border)] bg-[var(--color-bg-tertiary)] text-[var(--color-text-primary)]" - /> -

Must be within forecast range (tomorrow or later)

-
-
-
- - {/* Rocket */} - 🚀} defaultExpanded={false}> -
- setRocketConfig(prev => ({ ...prev, airframe_mass: parseFloat(v) || 0 }))} - unit="kg" - /> - setRocketConfig(prev => ({ ...prev, engine_mass: parseFloat(v) || 0 }))} - unit="kg" - /> - setRocketConfig(prev => ({ ...prev, lox_tank_structure_mass: parseFloat(v) || 0 }))} - unit="kg" - /> - setRocketConfig(prev => ({ ...prev, fuel_tank_structure_mass: parseFloat(v) || 0 }))} - unit="kg" - /> - setRocketConfig(prev => ({ ...prev, radius: parseFloat(v) || 0 }))} - unit="m" - /> - setRocketConfig(prev => ({ ...prev, rocket_length: parseFloat(v) || 0 }))} - unit="m" - /> - setRocketConfig(prev => ({ ...prev, motor_position: parseFloat(v) || 0 }))} - unit="m" - /> -
-
- - {/* Tanks */} - 🛢️} defaultExpanded={false}> -
-
-

LOX Tank

-
- setLoxTankConfig(prev => ({ ...prev, height: parseFloat(v) || 0 }))} - unit="m" - /> - setLoxTankConfig(prev => ({ ...prev, radius: parseFloat(v) || 0 }))} - unit="m" - /> - setLoxTankConfig(prev => ({ ...prev, position: parseFloat(v) || 0 }))} - unit="m" - /> -
-
-
-

Fuel Tank

-
- setFuelTankConfig(prev => ({ ...prev, height: parseFloat(v) || 0 }))} - unit="m" - /> - setFuelTankConfig(prev => ({ ...prev, radius: parseFloat(v) || 0 }))} - unit="m" - /> - setFuelTankConfig(prev => ({ ...prev, position: parseFloat(v) || 0 }))} - unit="m" - /> -
-
-
-
-
- - {/* Run Button */} -
- {isDirty && ( - - Unsaved DR - auto-save on Run - - )} - -
- - {/* Error */} - {error && ( -
-

❌ Error

-

{error}

-
- )} - - {/* Mass Adjustment Notification */} - {massAdjustments && (massAdjustments.loxAdjusted || massAdjustments.fuelAdjusted) && ( -
-

⚡ Propellant Masses Auto-Adjusted

-

- Propellant masses were capped to prevent tank overfill: -

-
    - {massAdjustments.loxAdjusted && ( -
  • - • LOX: {massAdjustments.loxOriginal?.toFixed(1)} kg → {massAdjustments.loxAdjusted.toFixed(1)} kg - (reduced by {((massAdjustments.loxOriginal ?? 0) - massAdjustments.loxAdjusted).toFixed(1)} kg) -
  • - )} - {massAdjustments.fuelAdjusted && ( -
  • - • Fuel: {massAdjustments.fuelOriginal?.toFixed(1)} kg → {massAdjustments.fuelAdjusted.toFixed(1)} kg - (reduced by {((massAdjustments.fuelOriginal ?? 0) - massAdjustments.fuelAdjusted).toFixed(1)} kg) -
  • - )} -
-

- 💡 Increase tank height/radius in Tank Geometry to fit more propellant. -

-
- )} - - {/* Results */} - {results && ( -
- {/* Key Metrics */} -
-

📊 Flight Results

- -
- - - - -
- -
- - -
- - {results.truncation?.truncated && ( -
-

- ⚠️ Burn truncated due to {results.truncation.reason} at {results.truncation.cutoff_time?.toFixed(2)}s -

-
- )} - - {results.error && ( -
-

⚠️ {results.error}

-
- )} -
- - {/* Charts */} - {altitudeChartData.length > 0 && ( -
-

📈 Altitude vs Time

-
- - - - - - - - - - - -
-
- )} - - {velocityChartData.length > 0 && ( -
-

📈 Velocity vs Time

-
- - - - - - - - - - -
-
- )} - - {/* Rocket Diagram */} - {results.rocket_diagram && ( -
-

🚀 Rocket Diagram

-
- Rocket diagram -
-
- )} -
- )} -
- ); -} diff --git a/EngineDesign/frontend/src/components/Optimizer.tsx b/EngineDesign/frontend/src/components/Optimizer.tsx index e8cfca739..176e61b64 100644 --- a/EngineDesign/frontend/src/components/Optimizer.tsx +++ b/EngineDesign/frontend/src/components/Optimizer.tsx @@ -3,7 +3,8 @@ import { DesignRequirements, DEFAULT_DESIGN_REQUIREMENTS } from './DesignRequire import { Layer1Optimization } from './Layer1Optimization'; import { Layer2Optimization } from './Layer2Optimization'; import { Layer3Optimization } from './Layer3Optimization'; -import { Layer4Optimization } from './Layer4Optimization'; +import { FlightSimulation } from './FlightSimulation'; +import { emitConfigChanged } from '../lib/configBus'; import { useReadOnly } from '@stardesign-ui'; import { useViewState } from '../lib/viewState'; import { @@ -98,7 +99,7 @@ export function Optimizer({ config }: OptimizerProps) {
{/* Main Header */}
-

🚀 Engine Design Optimization

+

Engine Design Optimization

Goal: Size optimal injector and chamber geometry to meet your:

@@ -116,7 +117,7 @@ export function Optimizer({ config }: OptimizerProps) { : 'bg-red-500/10 border-red-500/30 text-red-400' }`}>

- {saveStatus.type === 'success' ? '✅' : '❌'} {saveStatus.message} + {saveStatus.message}

)} @@ -130,7 +131,7 @@ export function Optimizer({ config }: OptimizerProps) { : 'border-transparent text-[var(--color-text-secondary)] hover:text-[var(--color-text-primary)] hover:border-[var(--color-border)]' }`} > - 📋 Design Requirements + Design Requirements @@ -177,6 +178,7 @@ export function Optimizer({ config }: OptimizerProps) { requirements={requirements} onRequirementsChange={setRequirements} onSave={handleSave} + config={config} />
@@ -201,10 +203,14 @@ export function Optimizer({ config }: OptimizerProps) { />
- emitConfigChanged(c)} + sliceKey="layer4" />
diff --git a/EngineDesign/frontend/src/components/PressureCurveChart.tsx b/EngineDesign/frontend/src/components/PressureCurveChart.tsx index ea3a3dfdf..8ab54c921 100644 --- a/EngineDesign/frontend/src/components/PressureCurveChart.tsx +++ b/EngineDesign/frontend/src/components/PressureCurveChart.tsx @@ -314,7 +314,7 @@ export function PressureCurveChart({ data, summary }: PressureCurveChartProps) { Thrust vs Time

- + - + - + - + - + - + - + - + - + - + - + - + {data.lox_mass_remaining_kg && ( - + - + - + - + - + - + {data.recession_rate_ablative_um_s && ( - + - + {data.recession_cumulative_ablative_mm && ( - + - + 🟢 - : 🔴 + } >
diff --git a/EngineDesign/frontend/src/components/TimeSeriesMode.tsx b/EngineDesign/frontend/src/components/TimeSeriesMode.tsx index ebf170cae..bac874521 100644 --- a/EngineDesign/frontend/src/components/TimeSeriesMode.tsx +++ b/EngineDesign/frontend/src/components/TimeSeriesMode.tsx @@ -1,4 +1,4 @@ -import { useState, useCallback, useEffect } from 'react'; +import { useState, useCallback, useEffect, useRef } from 'react'; import { PressureProfileForm } from './PressureProfileForm'; import { SegmentCurveBuilder } from './SegmentCurveBuilder'; import { PressureCurveChart } from './PressureCurveChart'; @@ -141,6 +141,7 @@ export function TimeSeriesMode({ config, onConfigLoaded }: TimeSeriesModeProps) const [duration, setDuration] = useState(() => getConfigBurnTime(config)); const [nSteps, setNSteps] = useState(101); const [loxProfile, setLoxProfile] = useState(defaultLoxProfile); + const lastPressureKey = useRef(null); const [fuelProfile, setFuelProfile] = useState(defaultFuelProfile); // Segment builder state @@ -232,6 +233,26 @@ export function TimeSeriesMode({ config, onConfigLoaded }: TimeSeriesModeProps) if (typeof fuelTank?.initial_pressure_psi === 'number') { setFuelInitialPressure(fuelTank.initial_pressure_psi); } + + // The Simple Profile starts from the design's own tank pressures. It used to start from a + // 750/600 psi constant regardless of the loaded config, so the first time-series a user ran + // was for some other rocket's tanks. Re-seeded only when the config's pressures change, so a + // profile the user is editing is not clobbered by an unrelated autosave. + const pKey = `${loxTank?.initial_pressure_psi ?? ''}|${fuelTank?.initial_pressure_psi ?? ''}`; + if (lastPressureKey.current !== pKey) { + lastPressureKey.current = pKey; + const seed = (p0: unknown, set: typeof setLoxProfile) => { + if (typeof p0 !== 'number' || !(p0 > 0)) return; + const start = Math.round(p0 * 10) / 10; + set((prev) => ({ + ...prev, + start_pressure_psi: start, + end_pressure_psi: prev.end_pressure_psi < start ? prev.end_pressure_psi : Math.round(0.7 * start), + })); + }; + seed(loxTank?.initial_pressure_psi, setLoxProfile); + seed(fuelTank?.initial_pressure_psi, setFuelProfile); + } }, [config]); // Sync local input states diff --git a/EngineDesign/frontend/src/lib/gating.test.ts b/EngineDesign/frontend/src/lib/gating.test.ts index 74474c97b..a8440e82d 100644 --- a/EngineDesign/frontend/src/lib/gating.test.ts +++ b/EngineDesign/frontend/src/lib/gating.test.ts @@ -90,8 +90,6 @@ const VIEW_ONLY: Record = { 'Optimizer.tsx:setActiveSubTab': 'which optimizer sub-tab is shown', 'DemoLayerCard.tsx:setIsExpanded': 'expand/collapse a layer card', 'FlightSimulation.tsx:setIsExpanded': 'expand/collapse a section', - 'Layer4Optimization.tsx:setIsExpanded': 'expand/collapse a section', - 'Layer4Optimization.tsx:onClick={runSimulation}': 'runs the flight sim; writes nothing', 'ForwardMode.tsx:handleEvaluate()': 'runs evaluate(); reads the config, never writes it', 'FlightSimulation.tsx:handleOptimize : handleSimulate': 'runs the flight sim; writes nothing', 'CustomPlotter.tsx:setShowDataPreview': 'shows the raw data table under the chart', diff --git a/EngineDesign/tests/test_flight_propellant_iteration.py b/EngineDesign/tests/test_flight_propellant_iteration.py index 3ba43c796..5871bfceb 100644 --- a/EngineDesign/tests/test_flight_propellant_iteration.py +++ b/EngineDesign/tests/test_flight_propellant_iteration.py @@ -1,13 +1,19 @@ -"""Flight sim propellant iteration — regression tests. +"""Flight-sim propellant iteration -- regression tests against the shipped default config. -Verifies that: -- Apogee rises when propellant increases below truncation threshold -- Apogee falls (or plateaus) when excess propellant is loaded above full-burn requirement -- Tank mass caps are surfaced (fuel tank smaller than config mass) +- Tank caps: a fuel tank too small for the requested load is capped by the flight router, a tank + that fits the load is not, and an explicit design_requirements capacity wins over geometry. +- Apogee rises with propellant while the burn is truncated. +- Excess propellant above the full-burn requirement does not raise apogee. +- A longer burn with enough propellant raises impulse and apogee. + +The previous version of this module read a config out of one developer's ~/Downloads folder and +asserted a property of that file's tank, so it skipped on every other machine and failed on his +once the file changed. Everything here is built from configs/default.yaml. """ from __future__ import annotations +import copy from pathlib import Path import numpy as np @@ -15,36 +21,30 @@ pytest.importorskip("rocketpy") -USER_CONFIG = Path("/Users/carlton/Downloads/Liquid Engine Designer Config.yaml") +CONFIG = Path(__file__).resolve().parents[1] / "configs" / "default.yaml" -def _load_user_config(): - if not USER_CONFIG.is_file(): - pytest.skip(f"User config not found: {USER_CONFIG}") +@pytest.fixture(scope="module") +def config(): from engine.pipeline.io import load_config - return load_config(str(USER_CONFIG)) + cfg = load_config(str(CONFIG)) + # Tanks sized so that no volume cap can bite in the propellant-response tests below; the cap + # logic has its own test. (default.yaml's fuel tank holds ~4.7 kg of methane at 100%.) + for tank, h in ((cfg.lox_tank, "lox_h"), (cfg.fuel_tank, "rp1_h")): + tank.tank_volume_m3 = None + setattr(tank, h, 2.0) + return cfg -def _simple_pressure_profiles(duration_s: float, n_points: int = 101): +def _pressure_profiles(config, duration_s: float, n_points: int = 101): + """Exponential blowdown from each tank's configured initial pressure to 70% of it.""" from engine.pipeline.time_series import generate_pressure_profile - times, lox_psi = generate_pressure_profile( - "exponential", - 560.8013772234059, - 400.0, - duration_s, - n_points, - decay_constant=3.0, - ) - _, fuel_psi = generate_pressure_profile( - "exponential", - 568.9552250017153, - 350.0, - duration_s, - n_points, - decay_constant=3.0, - ) + P_O0 = float(config.lox_tank.initial_pressure_psi or config.design_requirements.max_lox_tank_pressure_psi) + P_F0 = float(config.fuel_tank.initial_pressure_psi or config.design_requirements.max_fuel_tank_pressure_psi) + times, lox_psi = generate_pressure_profile("exponential", P_O0, 0.7 * P_O0, duration_s, n_points, decay_constant=3.0) + _, fuel_psi = generate_pressure_profile("exponential", P_F0, 0.7 * P_F0, duration_s, n_points, decay_constant=3.0) return times, lox_psi, fuel_psi @@ -53,20 +53,11 @@ def _run_timeseries(config, duration_s: float): from backend.routers.timeseries import compute_timeseries_results runner = PintleEngineRunner(config) - times, lox_psi, fuel_psi = _simple_pressure_profiles(duration_s) - data, summary = compute_timeseries_results( - runner, - times, - lox_psi, - fuel_psi, - run_copv=False, - ) - return data, summary + times, lox_psi, fuel_psi = _pressure_profiles(config, duration_s) + return compute_timeseries_results(runner, times, lox_psi, fuel_psi, run_copv=False) def _run_flight(config, data, lox_kg: float, fuel_kg: float): - import copy - from engine.optimizer.copv_flight_helpers import run_flight_simulation cfg = copy.deepcopy(config) @@ -74,43 +65,60 @@ def _run_flight(config, data, lox_kg: float, fuel_kg: float): cfg.fuel_tank.mass = fuel_kg times = np.asarray(data["time"], dtype=float) times = times - times[0] - burn_time = float(times[-1]) pressure_curves = { "time": times, "thrust": np.asarray(data["thrust_kN"], dtype=float) * 1000.0, "mdot_O": np.asarray(data["mdot_O_kg_s"], dtype=float), "mdot_F": np.asarray(data["mdot_F_kg_s"], dtype=float), } - return run_flight_simulation(cfg, pressure_curves, burn_time) + return run_flight_simulation(cfg, pressure_curves, float(times[-1])) -def test_fuel_tank_cap_below_config_mass(): - config = _load_user_config() +def test_fuel_mass_is_capped_to_tank_and_left_alone_when_it_fits(config): + """The flight router caps a load the tank cannot hold, matches the shared resolver, does not + touch a load that fits, and honours an explicit design_requirements capacity over geometry.""" + from backend.routers.flight import _apply_propellant_mass_caps + from engine.pipeline.config_schemas import PintleEngineConfig from engine.pipeline.tank_capacity import resolve_fuel_tank_limits - fuel_density = float(config.fluids["fuel"].density) - max_fill, _, fill_factor, _ = resolve_fuel_tank_limits(config, fuel_density) - config_mass = float(config.fuel_tank.mass) - assert config_mass > max_fill, ( - f"Expected config fuel mass ({config_mass} kg) to exceed tank cap ({max_fill:.2f} kg at {fill_factor*100:.0f}% fill)" - ) - - -def test_apogee_increases_with_propellant_when_truncated(): - config = _load_user_config() + rho_F = float(config.fluids["fuel"].density) + base = config.model_dump() + requested = 7.0 + + small = copy.deepcopy(base) + small["fuel_tank"].update(mass=requested, tank_volume_m3=None, rp1_h=0.3, rp1_radius=0.05) + adj, _, fuel_max, fill_factor = _apply_propellant_mass_caps(small, config) + expected_max, _, expected_ff, explicit = resolve_fuel_tank_limits(PintleEngineConfig(**small), rho_F) + assert not explicit + assert fuel_max == pytest.approx(expected_max) and fill_factor == pytest.approx(expected_ff) + assert fuel_max < requested, "test premise: this tank must be too small for the load" + assert adj["fuel"]["was_capped"] is True + assert small["fuel_tank"]["mass"] == pytest.approx(fuel_max), "router must write the capped mass back" + assert adj["fuel"]["original"] == pytest.approx(requested) and adj["fuel"]["capped"] == pytest.approx(fuel_max) + + big = copy.deepcopy(base) + big["fuel_tank"].update(mass=requested, tank_volume_m3=None, rp1_h=2.0, rp1_radius=0.15) + adj, _, fuel_max, _ = _apply_propellant_mass_caps(big, config) + assert fuel_max > requested + assert adj["fuel"]["was_capped"] is False + assert big["fuel_tank"]["mass"] == pytest.approx(requested), "a load that fits must not be touched" + + capped = copy.deepcopy(big) + capped["design_requirements"]["fuel_tank_capacity_kg"] = 3.0 + adj, _, fuel_max, _ = _apply_propellant_mass_caps(capped, config) + assert fuel_max == pytest.approx(3.0) and adj["fuel"]["explicit_capacity_kg"] == pytest.approx(3.0) + assert adj["fuel"]["was_capped"] is True and capped["fuel_tank"]["mass"] == pytest.approx(3.0) + + +def test_apogee_increases_with_propellant_when_truncated(config): data, summary = _run_timeseries(config, duration_s=6.8) lox_required = float(summary.get("lox_propellant_kg") or 0) fuel_required = float(summary.get("fuel_propellant_kg") or 0) assert lox_required > 0 and fuel_required > 0 - low_lox = max(0.5, lox_required * 0.45) - low_fuel = max(0.3, fuel_required * 0.45) - high_lox = lox_required * 1.05 - high_fuel = fuel_required * 1.05 - - low = _run_flight(config, data, low_lox, low_fuel) - high = _run_flight(config, data, high_lox, high_fuel) + low = _run_flight(config, data, max(0.5, lox_required * 0.45), max(0.3, fuel_required * 0.45)) + high = _run_flight(config, data, lox_required * 1.05, fuel_required * 1.05) assert low.get("success"), low.get("error") assert high.get("success"), high.get("error") @@ -120,9 +128,7 @@ def test_apogee_increases_with_propellant_when_truncated(): ) -def test_excess_propellant_does_not_increase_apogee(): - config = _load_user_config() - # Use a shorter burn so full-burn is achievable within tank caps +def test_excess_propellant_does_not_increase_apogee(config): data, summary = _run_timeseries(config, duration_s=3.5) lox_required = float(summary["lox_propellant_kg"]) @@ -134,14 +140,12 @@ def test_excess_propellant_does_not_increase_apogee(): assert optimal.get("success"), optimal.get("error") assert heavy.get("success"), heavy.get("error") assert optimal["truncation_info"].get("truncated") is False - assert heavy["apogee"] <= optimal["apogee"] + 15.0, ( f"Excess propellant should not increase apogee: optimal={optimal['apogee']:.1f}m heavy={heavy['apogee']:.1f}m" ) -def test_longer_burn_time_changes_impulse_and_apogee_with_enough_propellant(): - config = _load_user_config() +def test_longer_burn_time_changes_impulse_and_apogee_with_enough_propellant(config): short_data, short_summary = _run_timeseries(config, duration_s=4.0) long_data, long_summary = _run_timeseries(config, duration_s=8.0) @@ -150,16 +154,12 @@ def test_longer_burn_time_changes_impulse_and_apogee_with_enough_propellant(): assert long_imp > short_imp * 1.15, "Longer burn should deliver materially more impulse" short_flight = _run_flight( - config, - short_data, - float(short_summary["lox_propellant_kg"]) * 1.1, - float(short_summary["fuel_propellant_kg"]) * 1.1, + config, short_data, + float(short_summary["lox_propellant_kg"]) * 1.1, float(short_summary["fuel_propellant_kg"]) * 1.1, ) long_flight = _run_flight( - config, - long_data, - float(long_summary["lox_propellant_kg"]) * 1.1, - float(long_summary["fuel_propellant_kg"]) * 1.1, + config, long_data, + float(long_summary["lox_propellant_kg"]) * 1.1, float(long_summary["fuel_propellant_kg"]) * 1.1, ) assert short_flight.get("success"), short_flight.get("error") diff --git a/EngineDesign/tests/test_layer1_of_target_range.py b/EngineDesign/tests/test_layer1_of_target_range.py new file mode 100644 index 000000000..bbc475d52 --- /dev/null +++ b/EngineDesign/tests/test_layer1_of_target_range.py @@ -0,0 +1,39 @@ +"""Layer 1 refuses a target O/F outside the propellant's CEA table instead of optimizing garbage.""" + +import pytest + +from engine.optimizer.layers.layer1_static_optimization import _layer1_check_of_target_in_cea_range +from engine.pipeline.config_switch import load_canonical_config, apply_propellant +from engine.pipeline.config_schemas import PintleEngineConfig + + +def _ethalox_pintle(): + return PintleEngineConfig(**load_canonical_config("pintle")) # ethalox, MR_range [1.0, 2.5] + + +def test_in_range_target_passes(): + cfg = _ethalox_pintle() + lo, hi = cfg.combustion.cea.MR_range + _layer1_check_of_target_in_cea_range(cfg, 0.5 * (lo + hi)) + _layer1_check_of_target_in_cea_range(cfg, lo) + _layer1_check_of_target_in_cea_range(cfg, hi) + + +def test_out_of_range_target_is_refused_with_the_fix_in_the_message(): + cfg = _ethalox_pintle() + with pytest.raises(ValueError, match=r"outside the CEA table for ethalox .*\[1\.00, 2\.50\]"): + _layer1_check_of_target_in_cea_range(cfg, 3.5) + + +def test_stale_target_after_propellant_switch_is_caught(): + """The ethalox canonical carries optimal_of_ratio 1.4; overlay methalox (table 2.4-4.2) and the + untouched target must be refused, not silently pinned at the table edge.""" + switched = PintleEngineConfig(**apply_propellant(load_canonical_config("pintle"), "methalox")) + with pytest.raises(ValueError, match="outside the CEA table for methalox"): + _layer1_check_of_target_in_cea_range(switched, switched.design_requirements.optimal_of_ratio) + + +def test_missing_table_range_does_not_block(): + cfg = _ethalox_pintle() + cfg.combustion.cea.MR_range = None + _layer1_check_of_target_in_cea_range(cfg, 99.0) diff --git a/lib/stardesign-ui/src/useCheckout.ts b/lib/stardesign-ui/src/useCheckout.ts index eac8b452c..74d8b7444 100644 --- a/lib/stardesign-ui/src/useCheckout.ts +++ b/lib/stardesign-ui/src/useCheckout.ts @@ -10,6 +10,7 @@ * its inputs on. Greyed fields, not merely a refused save -- so there is no * state where someone believes they have it and does not. * - **It lapses** after inactivity server-side, and is released on tab close. + * While held, a heartbeat re-takes it well inside the server's lock_ttl. * * Two things in here exist to prevent data loss rather than to be tidy: * @@ -60,6 +61,8 @@ export interface UseCheckoutOptions { reload?: () => Promise | void; /** How often to re-check while somebody else holds it. */ pollMs?: number; + /** How often to re-stamp the checkout while we hold it (must beat the server's lock_ttl). */ + heartbeatMs?: number; } export function useCheckout({ @@ -67,6 +70,7 @@ export function useCheckout({ ref, reload, pollMs = 10_000, + heartbeatMs = 60_000, }: UseCheckoutOptions): Checkout { const [state, setState] = useState(FREE); const [busy, setBusy] = useState(false); @@ -87,9 +91,41 @@ export function useCheckout({ setError(null); }, [key]); + // Keep it alive while we hold it. The server treats a checkout with no + // heartbeat inside lock_ttl (5 min) as free, and autosave only writes when the + // design CHANGED -- so a long read, or a long optimizer run, let the token + // lapse under the holder and the next write came back 423 ("your checkout + // has lapsed"). Re-taking is the heartbeat: for the holder it is idempotent + // and re-stamps lockHeartbeat. A 423 here means somebody took it in a gap; + // drop to read-only at once instead of finding out on the next save. + useEffect(() => { + if (!ref || !state.lockedByMe) return; + let cancelled = false; + const beat = () => { + const r = refRef.current; + if (!r) return; + api + .takeCheckout(r) + .then((s) => !cancelled && setState(s)) + .catch((e) => { + if (cancelled) return; + if (e instanceof ApiError && e.status === 423) { + setState((prev) => ({ ...prev, lockedByMe: false })); + api.getCheckout(r).then((s) => !cancelled && setState(s)).catch(() => {}); + } + /* anything else is transient (network blip); the next beat retries */ + }); + }; + const id = setInterval(beat, heartbeatMs); + return () => { + cancelled = true; + clearInterval(id); + }; + }, [api, key, state.lockedByMe, heartbeatMs]); // eslint-disable-line react-hooks/exhaustive-deps + // Poll only while we do NOT hold it. A chip reading "taken" after the holder // has released is worse than no chip; this is what makes Take light up on its - // own. Once we hold it there is nothing to learn -- our own saves keep it. + // own. While we hold it the heartbeat above is what keeps the state honest. useEffect(() => { if (!ref || state.lockedByMe) return; let cancelled = false; @@ -110,20 +146,19 @@ export function useCheckout({ }, [api, key, state.lockedByMe, pollMs]); // eslint-disable-line react-hooks/exhaustive-deps // Give it back when the tab goes away, so a colleague is not left waiting out - // the inactivity timeout for a design nobody has open. + // the inactivity timeout for a design nobody has open. Only on pagehide: this + // used to fire on visibilitychange too, which released the checkout the moment + // the user looked at another tab -- mid-optimization, that turned the run's + // result write into "Take Design 1 before saving". A hidden tab still holds + // and heartbeats; a closed one lapses after lock_ttl. useEffect(() => { const drop = () => { const r = refRef.current; if (r && heldRef.current) api.releaseCheckoutOnUnload(r); }; - const onHide = () => { - if (document.visibilityState === 'hidden') drop(); - }; window.addEventListener('pagehide', drop); - document.addEventListener('visibilitychange', onHide); return () => { window.removeEventListener('pagehide', drop); - document.removeEventListener('visibilitychange', onHide); }; }, [api]);