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/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/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 2fc34bd4b..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) @@ -1739,6 +1741,115 @@ 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_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``. @@ -2575,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: @@ -2732,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: @@ -2877,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. @@ -2987,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 @@ -3039,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 @@ -3150,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): @@ -3327,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( @@ -3865,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, @@ -3975,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, @@ -4075,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]: @@ -7697,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: @@ -8402,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) @@ -8563,7 +8791,11 @@ 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) + _layer1_warn_stale_pressure_curves(optimized_config, 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/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/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/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/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/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 4aa7f619c..a8440e82d 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', @@ -89,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/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() 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_derived_dofs.py b/EngineDesign/tests/test_layer1_derived_dofs.py index ac6cbdb74..e0a68d0d9 100644 --- a/EngineDesign/tests/test_layer1_derived_dofs.py +++ b/EngineDesign/tests/test_layer1_derived_dofs.py @@ -804,3 +804,97 @@ 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) + + +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] 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]);