4405 cea_cache: Any =
None,
4407 eval_error: Optional[str] =
None,
4408 return_terms: bool =
False,
4410 """The Layer-1 objective; with ``return_terms`` every weighted term and the values behind it.
4412 The ONE scoring function. The worker pool (``_eval_candidate``) and the inline
4413 ``objective()`` closure both call it on the SOLVED vector, so CMA, the block stage,
4414 L-BFGS-B and the reported breakdown rank the same surface. It used to exist twice and the
4415 two drifted by ~300 points on the same design (chamber length, ΔP centre term, L* keys).
4417 Structure: infeasibility (BASE_INFEAS + W_INFEAS x residual) dominates; below it the
4418 requirement terms (thrust, Pc, O/F with deadbands), the preference bands, and the
4419 figure of merit -- delivered Isp (``isp_penalty``) against chamber mass
4420 (``chamber_mass_penalty``).
4422 Returns the scalar, or with ``return_terms`` a dict: ``objective``, ``terms`` (weighted;
4423 sums to ``objective`` on a feasible design) and the context values.
4425 inj_type = constants.get(
"injector_type",
"pintle")
4426 x = np.asarray(x, dtype=float)
4428 idx_P_O = int(constants.get(
"idx_P_O", 8))
4429 idx_P_F = int(constants.get(
"idx_P_F", 9))
4431 target_thrust = float(constants.get(
"target_thrust", 1000))
4432 optimal_of = float(constants.get(
"optimal_of", 2.5))
4433 target_P_exit = float(constants.get(
"P_ambient", 101325.0))
4434 max_lox_P_psi = float(constants.get(
"max_lox_P_psi", 500))
4435 max_fuel_P_psi = float(constants.get(
"max_fuel_P_psi", 500))
4436 TOTAL_WALL_THICKNESS_M = float(constants.get(
"TOTAL_WALL_THICKNESS_M", 0.0254))
4437 max_nozzle_exit = float(constants.get(
"max_nozzle_exit", 1.0))
4438 min_stability = float(requirements.get(
"min_stability_margin", _LAYER1_DEFAULT_MIN_STABILITY_MARGIN))
4440 A_throat = float(x[0])
4442 expansion_ratio = float(x[2])
4444 float(x[3]), constants.get(
"layer1_chamber_od_increment_in", 0.0),
4445 wall_m=TOTAL_WALL_THICKNESS_M,
4446 snap_target=constants.get(
"layer1_chamber_od_snap_target",
"outer"),
4451 d_jet_O = d_jet_F = ang_O = ang_F = sp_O = sp_F = float(
"nan")
4452 if inj_type ==
"impinging":
4456 d_jet_O, ang_O, sp_O = float(x[5]), float(x[6]), float(x[7])
4457 d_jet_F, ang_F, sp_F = float(x[8]), float(x[9]), float(x[10])
4458 A_lox_injector = float(n_el_O * np.pi * (d_jet_O / 2.0) ** 2)
4459 A_fuel_injector = float(n_el_O * np.pi * (d_jet_F / 2.0) ** 2)
4461 d_pintle_tip = float(x[4])
4463 n_orifices = int(x[6])
4464 d_orifice = float(x[7])
4465 A_lox_injector = float(n_orifices * np.pi * (d_orifice / 2.0) ** 2)
4466 R_inner = float(d_pintle_tip / 2.0)
4467 A_fuel_injector = float(np.pi * ((R_inner + h_gap) ** 2 - R_inner ** 2))
4469 D_chamber_inner = D_chamber_outer - TOTAL_WALL_THICKNESS_M
4470 if D_chamber_inner <= 0:
4471 D_chamber_inner =
max(D_chamber_outer * 0.3, 0.01)
4472 A_chamber_check = np.pi * (D_chamber_inner / 2.0) ** 2
4473 A_throat_check = A_throat
4475 A_exit = A_throat * expansion_ratio
4477 A_exit = A_throat * 10.0
4478 D_exit_check = np.sqrt(
max(0.0, 4.0 * A_exit / np.pi))
4479 if D_exit_check > max_nozzle_exit:
4480 D_exit_check = max_nozzle_exit
4481 D_throat_check = np.sqrt(4.0 * A_throat_check / np.pi)
if A_throat_check > 0
else 0.0
4483 out: Dict[str, Any] = {
"eval_success":
False,
"terms": {}}
4486 A_lox_flow, A_fuel_flow = A_lox_injector, A_fuel_injector
4487 if inj_type
in (
"impinging",
"pintle")
and isinstance(result, dict):
4488 A_lox_flow, A_fuel_flow, eff_w = effective_flow_areas_from_cd(
4489 result.get(
"diagnostics"), A_lox_injector, A_fuel_injector)
4490 if eff_w
and isinstance(result.get(
"diagnostics"), dict):
4491 merge_effective_area_warnings(result[
"diagnostics"], eff_w)
4492 lox_ratio = A_lox_flow / A_throat_check
if A_throat_check > 0
else np.nan
4493 fuel_ratio = A_fuel_flow / A_throat_check
if A_throat_check > 0
else np.nan
4495 P_O_psi = float(x[idx_P_O])
4496 P_F_psi = float(x[idx_P_F])
4497 P_O_ratio = P_O_psi / max_lox_P_psi
if max_lox_P_psi > 0
else 0.0
4498 P_F_ratio = P_F_psi / max_fuel_P_psi
if max_fuel_P_psi > 0
else 0.0
4500 infeasibility_score = 0.0
4502 _tr = {}
if _INFEAS_TRACE_ON
else None
4504 _INFEAS_TRACE.append(_tr)
4506 if A_chamber_check > 0
and A_throat_check > 0:
4507 cr_check = A_chamber_check / A_throat_check
4508 infeasibility_score +=
max(0.0, (A_throat_check * 1.1) / A_chamber_check - 1.0) ** 2
4510 infeasibility_score += (1.5 - cr_check) ** 2
4511 elif cr_check > 15.0:
4512 infeasibility_score += (cr_check - 15.0) ** 2
4513 if inj_type ==
"pintle" and D_chamber_inner > 0:
4514 infeasibility_score +=
max(0.0, (d_pintle_tip * 1.1) / D_chamber_inner - 1.0) ** 2
4515 if D_throat_check > 0
and D_chamber_inner > 0:
4516 infeasibility_score +=
max(0.0, D_throat_check - D_chamber_inner * 0.95) ** 2 * 10.0
4517 if D_exit_check > 0
and D_throat_check > 0:
4518 infeasibility_score +=
max(0.0, D_throat_check / D_exit_check - 1.0) ** 2
4520 if inj_type ==
"pintle" and A_throat_check > 0:
4521 infeasibility_score +=
max(0.0, lox_ratio - 1.0) ** 2
4522 infeasibility_score +=
max(0.0, fuel_ratio - 1.0) ** 2
4524 area_ratio = A_lox_flow / A_fuel_flow
4525 Cd_ratio = float(constants.get(
"discharge_cd_inf_ratio", 1.0))
4526 rho_ratio = np.sqrt(
4527 max(float(constants.get(
"rho_oxidizer", 1140.0)), 1e-9)
4528 /
max(float(constants.get(
"rho_fuel", 780.0)), 1e-9))
4529 area_ratio_factor = Cd_ratio * rho_ratio * np.sqrt(1.2)
4530 required = optimal_of / area_ratio_factor
if area_ratio_factor > 0
else np.inf
4531 if required > 0
and np.isfinite(required):
4532 infeasibility_score +=
max(0.0, abs(area_ratio - required) / required - 0.5) ** 2
4533 elif inj_type ==
"impinging" and A_throat_check > 0:
4536 n_elements=n_el_O, spacing_O_m=sp_O, spacing_F_m=sp_F,
4537 d_jet_O_m=d_jet_O, d_jet_F_m=d_jet_F, D_chamber_inner_m=D_chamber_inner,
4538 angle_O_deg=ang_O, angle_F_deg=ang_F,
4539 center_clear_dia_m=float(constants.get(
"layer1_injector_center_clear_dia_m", 0.0)
or 0.0),
4540 min_web_m=float(constants.get(
"layer1_injector_min_web_m", 0.0)
or 0.0),
4541 wall_clearance_m=float(constants.get(
"layer1_injector_wall_clearance_m", 0.0)
or 0.0),
4542 spray_radius_frac=float(constants.get(
"layer1_injector_spray_radius_frac", 0.0)
or 0.0),
4543 spray_radius_tol=float(constants.get(
"layer1_injector_spray_radius_tol", 0.08)
or 0.08),
4544 min_face_incidence_deg=float(
4545 constants.get(
"layer1_injector_min_face_incidence_deg", 0.0)
or 0.0),
4546 face_contoured=bool(constants.get(
"layer1_injector_face_contoured",
False)),
4547 exit_land=float(constants.get(
"layer1_injector_exit_land", EXIT_LAND_DEFAULT)),
4550 n_elements=n_el_O, spacing_O_m=sp_O, spacing_F_m=sp_F,
4551 d_jet_O_m=d_jet_O, d_jet_F_m=d_jet_F, D_chamber_inner_m=D_chamber_inner,
4552 angle_O_deg=ang_O, angle_F_deg=ang_F,
4553 plate_thickness_m=float(constants.get(
"layer1_injector_plate_thickness_m", 0.0)
or 0.0),
4554 counterbore_dia_m=float(constants.get(
"layer1_injector_counterbore_dia_m", 0.0)
or 0.0),
4555 land_ld_O=float(constants.get(
"layer1_injector_land_ld_O", 4.0)
or 4.0),
4556 land_ld_F=float(constants.get(
"layer1_injector_land_ld_F", 4.0)
or 4.0),
4557 min_back_web_m=float(constants.get(
"layer1_injector_min_back_web_m", 0.0)
or 0.0),
4558 back_keepout_dia_m=float(constants.get(
"layer1_injector_back_keepout_dia_m", 0.0)
or 0.0),
4559 back_keepout_wall_m=float(constants.get(
"layer1_injector_back_keepout_wall_m", 0.0)
or 0.0),
4560 face_contoured=bool(constants.get(
"layer1_injector_face_contoured",
False)),
4561 exit_land=float(constants.get(
"layer1_injector_exit_land", EXIT_LAND_DEFAULT)),
4562 back_channels=bool(constants.get(
"layer1_injector_back_channels",
False)),
4564 passage_ld_O=float(constants.get(
"layer1_injector_land_ld_O", 4.0)
or 4.0),
4565 passage_ld_F=float(constants.get(
"layer1_injector_land_ld_F", 4.0)
or 4.0),
4566 channel_width_m=float(constants.get(
"layer1_injector_channel_width_m", 0.0)
or 0.0),
4567 channel_floor=str(constants.get(
"layer1_injector_channel_floor",
"flat")),
4568 r_plate_m=float(constants.get(
"layer1_injector_r_plate_m", 0.0)
or 0.0),
4569 r_port_m=float(constants.get(
"layer1_injector_r_port_m", 0.0)
or 0.0),
4570 liner_thickness_m=float(constants.get(
"layer1_injector_liner_thickness_m", 0.0)
or 0.0),
4574 n_elements=n_el_O, spacing_O_m=sp_O, spacing_F_m=sp_F,
4575 d_jet_O_m=d_jet_O, d_jet_F_m=d_jet_F, angle_O_deg=ang_O, angle_F_deg=ang_F)
4580 if bool(constants.get(
"layer1_injector_back_channels",
False))
and isinstance(result, dict):
4581 _dg = result.get(
"diagnostics")
if isinstance(result.get(
"diagnostics"), dict)
else {}
4583 n_elements=n_el_O, spacing_O_m=sp_O, spacing_F_m=sp_F,
4584 d_jet_O_m=d_jet_O, d_jet_F_m=d_jet_F, angle_O_deg=ang_O, angle_F_deg=ang_F,
4585 plate_thickness_m=float(constants.get(
"layer1_injector_plate_thickness_m", 0.0)
or 0.0),
4586 face_contoured=bool(constants.get(
"layer1_injector_face_contoured",
False)),
4587 exit_land=float(constants.get(
"layer1_injector_exit_land", EXIT_LAND_DEFAULT)),
4588 passage_ld_O=float(constants.get(
"layer1_injector_land_ld_O", 4.0)
or 4.0),
4589 passage_ld_F=float(constants.get(
"layer1_injector_land_ld_F", 4.0)
or 4.0),
4590 channel_width_m=float(constants.get(
"layer1_injector_channel_width_m", 0.0)
or 0.0),
4591 channel_floor=str(constants.get(
"layer1_injector_channel_floor",
"flat")),
4592 channel_inlets=int(constants.get(
"layer1_injector_channel_inlets", 1)
or 1),
4593 mdot_O=result.get(
"mdot_O"), mdot_F=result.get(
"mdot_F"),
4594 dp_O=_dg.get(
"delta_p_injector_O"), dp_F=_dg.get(
"delta_p_injector_F"),
4595 rho_O=constants.get(
"rho_O"), rho_F=constants.get(
"rho_F"))
4599 d_jet_O=d_jet_O, d_jet_F=d_jet_F, sp_O=sp_O, sp_F=sp_F,
4600 n_el_O=n_el_O, n_el_F=n_el_O, D_chamber_inner=D_chamber_inner)
4602 A_lox_flow=A_lox_flow, A_fuel_flow=A_fuel_flow, A_throat_check=A_throat_check)
4607 eval_success = bool(result.get(
"success",
False))
if isinstance(result, dict)
else False
4608 if isinstance(result, dict)
and not eval_success:
4609 eval_success = bool(np.isfinite(float(result.get(
"F", np.nan)))
4610 and np.isfinite(float(result.get(
"Pc", np.nan))))
4612 "lox_ratio": lox_ratio,
"fuel_ratio": fuel_ratio,
4613 "D_chamber_inner": D_chamber_inner,
"D_throat": D_throat_check,
4616 if not eval_success:
4618 infeasibility_score += 1.0
4619 err = (eval_error
or "").lower()
4620 if (
"supply > demand" in err)
or (
"invalid bracket" in err)
or (
"no solution" in err):
4621 if np.isfinite(lox_ratio):
4622 infeasibility_score +=
max(0.0, lox_ratio - 0.90) ** 2
4623 if np.isfinite(fuel_ratio):
4624 infeasibility_score +=
max(0.0, fuel_ratio - 0.90) ** 2
4626 infeasibility_score +=
max(0.0, 0.90 - P_O_ratio) ** 2 +
max(0.0, 0.90 - P_F_ratio) ** 2
4627 if np.isfinite(geom_hard)
and geom_hard > 0.0:
4628 infeasibility_score += float(geom_hard)
4630 "objective": float(_LAYER1_BASE_INFEAS + _LAYER1_W_INFEAS * float(infeasibility_score)),
4631 "infeasibility_score": float(infeasibility_score),
4632 "inf_residual": float(infeasibility_score),
4633 "length_violation":
False,
4634 "thrust_error": 1.0,
"of_error": 1.0,
4635 "stability_state":
"unstable",
"stability_score": 0.0,
4636 "chugging_margin": 0.0,
"acoustic_margin": 0.0,
"feed_margin": 0.0,
4638 return out
if return_terms
else float(out[
"objective"])
4640 F_actual = float(result.get(
"F", np.nan))
4641 MR_actual = float(result.get(
"MR", np.nan))
4642 Pc_actual = float(result.get(
"Pc", np.nan))
4643 Isp_actual = float(result.get(
"Isp", np.nan))
4644 Cf_actual = float(result.get(
"Cf_actual", result.get(
"Cf", np.nan)))
4645 P_exit_actual = float(result.get(
"P_exit", np.nan))
4646 stability = result.get(
"stability_results", {})
or {}
4647 diagnostics = result.get(
"diagnostics")
if isinstance(result.get(
"diagnostics"), dict)
else {}
4649 thrust_error = (abs(F_actual - target_thrust) / target_thrust
4650 if (target_thrust > 0
and np.isfinite(F_actual))
else 1.0)
4652 _of_tol = float(constants.get(
"layer1_of_deadband_rel", 0.0)
or 0.0)
4653 of_error_raw = (abs(MR_actual - optimal_of) / optimal_of
4654 if (optimal_of > 0
and np.isfinite(MR_actual))
else 1.0)
4655 of_error =
max(0.0, of_error_raw - _of_tol)
4657 thrust_penalty_sq_term = 0.0
4658 if target_thrust > 0
and np.isfinite(F_actual):
4661 thrust_penalty_sq_term = excess ** 2
4663 thrust_penalty_sq_term = 1.0
4668 pc_penalty_term = 0.0
4669 pc_target_pa = constants.get(
"layer1_target_Pc_pa")
4670 if pc_target_pa
and float(pc_target_pa) > 0:
4671 if np.isfinite(Pc_actual):
4672 pc_penalty_term =
max(0.0, abs(Pc_actual - float(pc_target_pa)) / float(pc_target_pa)
4673 - _LAYER1_PC_AIM_REL)
4675 pc_penalty_term = 1.0
4677 _ins_q_raw = constants.get(
"layer1_exit_pressure_inside_quad_scale")
4678 _ins_q = (float(_ins_q_raw)
if _ins_q_raw
is not None and np.isfinite(float(_ins_q_raw))
4679 else _LAYER1_DEFAULT_EXIT_PRESSURE_INSIDE_QUAD)
4681 P_exit_actual, target_P_exit,
4682 deadband_rel=float(constants.get(
"layer1_exit_pressure_deadband_rel", 0.05)),
4683 inside_rel_quad_scale=_ins_q))
4686 stability_state = stability.get(
"stability_state",
"unstable")
4687 stability_score = float(stability.get(
"stability_score", 0.0))
4688 chugging_margin =
max(0.0, float(stability.get(
"chugging", {}).get(
"stability_margin", 0.0)))
4689 acoustic_margin =
max(0.0, float(stability.get(
"acoustic", {}).get(
"stability_margin", 0.0)))
4690 feed_margin =
max(0.0, float(stability.get(
"feed_system", {}).get(
"stability_margin", 0.0)))
4691 handicap = float(requirements.get(
"stability_margin_handicap", 0.0))
4692 effective_min_score = float(requirements.get(
"min_stability_score", 0.75)) *
max(0.0, 1.0 - handicap)
4693 effective_margin = float(min_stability) *
max(0.0, 1.0 - handicap)
4694 require_stable_state = bool(requirements.get(
"require_stable_state",
True))
4695 state_ok = (stability_state ==
"stable")
if require_stable_state
else (stability_state !=
"unstable")
4697 infeasibility_score += 1.0
4698 if effective_min_score > 0:
4699 infeasibility_score +=
max(0.0, (effective_min_score - stability_score) / effective_min_score) ** 2
4700 if effective_margin > 0:
4701 for _m
in (chugging_margin, acoustic_margin, feed_margin):
4702 infeasibility_score +=
max(0.0, (effective_margin - _m) / effective_margin) ** 2
4705 Cf_lo, Cf_hi = 1.3, 1.8
4706 _cf = float(Cf_actual)
if np.isfinite(Cf_actual)
else 0.0
4707 cf_hinge = (
max(0.0, Cf_lo - _cf) ** 2 +
max(0.0, _cf - Cf_hi) ** 2) / (Cf_hi - Cf_lo) ** 2
4712 L_chamber_curr = _lens[
"length"]
4713 L_cylindrical = _lens[
"length_cylindrical"]
4714 eps_eff = (np.pi * D_exit_check ** 2 / 4.0) / A_throat
if A_throat > 0
else expansion_ratio
4717 max_engine_length = float(requirements.get(
"max_engine_length", 0.50)
or 0.50)
4719 length_violation =
False
4720 if np.isfinite(L_engine_curr)
and max_engine_length > 0:
4721 if L_engine_curr > max_engine_length:
4722 length_violation =
True
4723 length_term = ((L_engine_curr - max_engine_length) / max_engine_length) ** 2
4725 length_term =
max(0.0, (L_engine_curr - max_engine_length * 0.9)
4726 / (max_engine_length * 0.1)) ** 2
4728 smd_eff_um =
_effective_smd_um(diagnostics, MR_actual
if np.isfinite(MR_actual)
else None)
4733 W_LSTAR = float(constants.get(
"layer1_W_LSTAR", 0.0)
or 0.0)
4735 lstar_target_curr = float(
"nan")
4738 if k
in (
"min_Lstar",
"max_Lstar")
and constants.get(k)
is not None:
4739 return constants.get(k)
4740 return requirements.get(k)
4742 smd_eff_um, requirements, _lstar_get,
4743 float(constants.get(
"layer1_Lstar_target_m", 1.0)
or 1.0))
4748 W_ISP = float(constants.get(
"layer1_W_ISP", 0.0)
or 0.0)
4750 isp_ref = float(constants.get(
"layer1_Isp_ref_s", np.nan)
or np.nan)
4753 W_MASS = float(constants.get(
"layer1_W_MASS", 0.0)
or 0.0)
4755 A_chamber_check, L_chamber_curr, TOTAL_WALL_THICKNESS_M,
4756 float(constants.get(
"layer1_chamber_wall_density_kg_m3", 2000.0)),
4757 Pc_pa=float(Pc_actual)
if np.isfinite(Pc_actual)
else 0.0,
4758 plate_support=str(constants.get(
"layer1_injector_plate_support",
"clamped")),
4760 r_plate_m=(0.5 * float(D_chamber_inner)
4761 + float(constants.get(
"layer1_injector_liner_thickness_m", 0.0)
or 0.0)
4762 if constants.get(
"layer1_injector_plate_support")
is not None else 0.0))
4764 chamber_mass_kg, float(constants.get(
"layer1_chamber_mass_ref_kg", 5.0)))
4765 if W_MASS > 0.0
else 0.0)
4769 constants, L_cylindrical=L_cylindrical, D_chamber_inner=D_chamber_inner,
4770 A_chamber=A_chamber_check, n_elements=n_el_O
if inj_type ==
"impinging" else 0)
4772 W_THRUST = float(constants.get(
"layer1_W_THRUST", 1e4))
4773 W_PC = float(constants.get(
"layer1_W_PC", 0.0))
4774 W_OF = float(constants.get(
"layer1_W_OF", 1e4))
4775 w_of_low_mr =
max(1.0, float(constants.get(
"layer1_W_OF_low_MR_scale", 1.0)))
4776 w_of_high_mr =
max(1.0, float(constants.get(
"layer1_W_OF_high_MR_scale", 1.0)))
4778 W_EXIT = float(constants.get(
"layer1_W_EXIT", 2.0e2))
4780 W_MOM = float(constants.get(
"W_MOM", 75.0))
4781 W_DP = float(constants.get(
"W_DP", 160.0))
4782 W_SMD = float(constants.get(
"W_SMD", 0.0))
4783 W_ANGLE = float(constants.get(
"W_IMPINGING_ANGLE", 0.0))
4784 W_JET_ASYM = float(constants.get(
"W_IMPINGING_JET_ASYM", 0.0))
4785 _jet_asym_mx = constants.get(
"layer1_impinging_jet_angle_max_asym_deg")
4786 jet_asym_max = float(_jet_asym_mx)
if _jet_asym_mx
is not None else 28.0
4787 W_CHAMBER_SHAPE = float(constants.get(
"W_CHAMBER_SHAPE", 2500.0))
4788 W_TANK_EQUAL = float(constants.get(
"W_TANK_EQUAL", 0.0))
4789 W_IMP_GEOM = float(constants.get(
"W_IMP_GEOM", 0.0))
4790 W_geom_ao_af = float(constants.get(
"W_geom_ao_af_momentum", 0.0))
4792 momentum_term = angle_term = jet_asym_term = tilt_term = 0.0
4793 W_TILT = float(constants.get(
"layer1_W_TILT", 0.0)
or 0.0)
4794 R_val = float(
"nan")
4795 tilt_deg = float(
"nan")
4796 tilt_breakeven = float(
"nan")
4797 imp_angle_deg = float(
"nan")
4798 if inj_type ==
"impinging":
4800 imp_angle_deg = float(diagnostics.get(
"impingement_angle_deg", np.nan))
4801 except (TypeError, ValueError):
4802 imp_angle_deg = float(
"nan")
4804 n_elements=n_el_O, spacing_O_m=sp_O, spacing_F_m=sp_F,
4805 angle_O_deg=ang_O, angle_F_deg=ang_F,
4806 D_chamber_inner_m=D_chamber_inner, L_chamber_m=L_chamber_curr)
4807 _Rv = diagnostics.get(
"momentum_ratio_R")
4808 if _Rv
is not None and np.isfinite(float(_Rv))
and float(_Rv) > 0:
4812 result.get(
"mdot_O"), result.get(
"mdot_F"),
4813 constants.get(
"rho_O", 1141.0), constants.get(
"rho_F", 789.0),
4814 n_el_O, d_jet_O, d_jet_F, ang_O, ang_F,
4815 lox_inboard=bool(constants.get(
"layer1_ring_order_fuel_outboard",
True)),
4816 spacing_O_m=sp_O, spacing_F_m=sp_F)
4820 from_reach=bool(constants.get(
"layer1_resultant_tilt_from_reach",
False)),
4821 constant_deg=float(constants.get(
"layer1_resultant_tilt_max_deg", 0.0)),
4822 breakeven_deg=tilt_breakeven,
4823 margin=float(constants.get(
"layer1_resultant_tilt_reach_margin", 1.5))),
4824 scale_deg=float(constants.get(
"layer1_resultant_tilt_scale_deg", 2.0)))
4830 scale_deg=float(constants.get(
"layer1_resultant_tilt_scale_deg", 2.0)),
4831 outward_multiplier=float(
4832 constants.get(
"layer1_resultant_tilt_outward_multiplier", 25.0)
or 25.0))
4834 _mlo = constants.get(
"impinging_momentum_R_min")
4835 _mhi = constants.get(
"impinging_momentum_R_max")
4838 wall_side_multiplier=float(constants.get(
"layer1_momentum_wall_side_multiplier", 1.0)),
4839 low_side_multiplier=float(constants.get(
"layer1_momentum_low_side_multiplier", 10.0)),
4840 band_width=float(constants.get(
"layer1_momentum_band_width", 0.05)),
4841 safe_side_is_below_one=bool(constants.get(
"layer1_ring_order_fuel_outboard",
True)),
4842 scale=float(constants.get(
"layer1_momentum_scale", 0.10)),
4843 band_lo=float(_mlo)
if _mlo
is not None else None,
4844 band_hi=float(_mhi)
if _mhi
is not None else None,
4847 _alo = constants.get(
"layer1_impinging_angle_deg_min")
4848 _ahi = constants.get(
"layer1_impinging_angle_deg_max")
4851 angle_band_lo_deg=float(_alo)
if _alo
is not None else None,
4852 angle_band_hi_deg=float(_ahi)
if _ahi
is not None else None)
4853 if W_JET_ASYM > 0.0:
4855 ang_O, ang_F, max_delta_deg=jet_asym_max)
4857 chamber_shape_term = 0.0
4858 dt_ratio = D_chamber_inner / D_throat_check
if (D_chamber_inner > 0
and D_throat_check > 0)
else float(
"nan")
4859 ld_ratio = L_chamber_curr / D_chamber_inner
if (D_chamber_inner > 0
and L_chamber_curr > 0)
else float(
"nan")
4860 if np.isfinite(dt_ratio):
4862 dt_ratio, float(constants.get(
"layer1_chamber_dt_ratio_min", 2.2)),
4863 float(constants.get(
"layer1_chamber_dt_ratio_max", 3.2)))
4864 if np.isfinite(ld_ratio):
4866 ld_ratio, float(constants.get(
"layer1_chamber_ld_ratio_min", 1.0)),
4867 float(constants.get(
"layer1_chamber_ld_ratio_max", 3.2)))
4869 geom_ao_af_term = 0.0
4870 geom_ao_af_ratio = expected_ao_af = float(
"nan")
4871 if inj_type ==
"impinging" and W_geom_ao_af > 0.0
and A_fuel_injector > 0:
4873 float(A_lox_injector), float(A_fuel_injector), float(optimal_of),
4874 float(constants.get(
"rho_oxidizer", 1140.0)), float(constants.get(
"rho_fuel", 422.6)))
4876 _of_tol_for_geom =
max(0.05, float(constants.get(
"layer1_of_validation_tol", 0.15)))
4877 geom_ao_af_scale =
min(1.0, float(of_error / _of_tol_for_geom) ** 2)
4878 geom_ao_af_weighted = W_geom_ao_af * geom_ao_af_term * geom_ao_af_scale
4880 injector_dp_weighted = 0.0
4881 ratio_o = ratio_f =
None
4882 if np.isfinite(Pc_actual)
and Pc_actual > 0:
4883 ratio_o, ratio_f = injector_dp_ratios_from_eval_result(Pc_actual, result)
4884 injector_dp_weighted = injector_dp_ratio_penalty_weighted(
4885 ratio_o, ratio_f, W_DP, float(constants.get(
"W_DP_HIGH", 480.0)),
4886 o_band=(float(constants.get(
"injector_dp_ratio_O_min", _LAYER1_DEFAULT_DP_O_BAND[0])),
4887 float(constants.get(
"injector_dp_ratio_O_max", _LAYER1_DEFAULT_DP_O_BAND[1]))),
4888 f_band=(float(constants.get(
"injector_dp_ratio_F_min", _LAYER1_DEFAULT_DP_F_BAND[0])),
4889 float(constants.get(
"injector_dp_ratio_F_max", _LAYER1_DEFAULT_DP_F_BAND[1]))),
4890 w_dp_o=float(constants.get(
"W_DP_O", W_DP)),
4891 w_dp_f=float(constants.get(
"W_DP_F", W_DP)),
4892 o_soft_floor=(float(constants[
"injector_dp_ratio_O_soft_floor"])
4893 if constants.get(
"injector_dp_ratio_O_soft_floor")
is not None else None),
4894 w_dp_o_floor=float(constants.get(
"W_DP_O_FLOOR", 0.0)
or 0.0),
4895 w_dp_center=float(constants.get(
"W_DP_CENTER", 500.0)),
4899 if inj_type ==
"impinging" and W_SMD > 0.0:
4902 target_smd_microns=float(constants.get(
"target_smd_microns", 50.0)),
4903 smd_rel_tol=float(constants.get(
"layer1_smd_rel_tol", 0.20)),
4904 mr_mass=float(MR_actual)
if np.isfinite(MR_actual)
and MR_actual > 0
else None)
4906 tank_equal_term = 0.0
4907 if W_TANK_EQUAL > 0.0:
4909 P_O_psi, P_F_psi, scale_psi=float(constants.get(
"layer1_tank_equal_scale_psi", 100.0)),
4910 tol_psi=float(constants.get(
"layer1_tank_equal_tol_psi", 0.0)
or 0.0),
4911 inband_frac=float(constants.get(
"layer1_tank_equal_inband_frac", 0.0)
or 0.0))
4914 if inj_type ==
"impinging" and W_IMP_GEOM > 0.0:
4916 diagnostics, D_chamber_inner_m=D_chamber_inner, L_chamber_m=L_chamber_curr)
4919 if constants.get(
"layer1_enforce_ring_geometry",
True):
4921 n_elements=n_el_O, spacing_O_m=sp_O, spacing_F_m=sp_F,
4922 d_jet_O_m=d_jet_O, d_jet_F_m=d_jet_F, D_chamber_inner_m=D_chamber_inner,
4923 angle_O_deg=ang_O, angle_F_deg=ang_F,
4924 Ld_min=float(constants.get(
"layer1_impingement_Ld_min", _LD_BAND_DEFAULT.lo)),
4925 Ld_max=float(constants.get(
"layer1_impingement_Ld_max", _LD_BAND_DEFAULT.hi)),
4926 ring_order_fuel_outboard=bool(constants.get(
"layer1_ring_order_fuel_outboard",
True)),
4927 center_clear_dia_m=float(constants.get(
"layer1_injector_center_clear_dia_m", 0.0)
or 0.0),
4928 min_web_m=float(constants.get(
"layer1_injector_min_web_m", 0.0)
or 0.0),
4929 wall_clearance_m=float(constants.get(
"layer1_injector_wall_clearance_m", 0.0)
or 0.0),
4930 spray_radius_frac=float(constants.get(
"layer1_injector_spray_radius_frac", 0.0)
or 0.0),
4931 spray_radius_tol=float(constants.get(
"layer1_injector_spray_radius_tol", 0.08)
or 0.08),
4932 face_contoured=bool(constants.get(
"layer1_injector_face_contoured",
False)),
4933 exit_land=float(constants.get(
"layer1_injector_exit_land", EXIT_LAND_DEFAULT)),
4936 if not np.isfinite(infeasibility_score)
or infeasibility_score < 0:
4937 infeasibility_score = 1.0
4940 if np.isfinite(geom_hard)
and geom_hard > 0.0:
4941 inf_residual = float(inf_residual) + float(geom_hard)
4943 of_sq = float(of_error) ** 2
4944 if np.isfinite(MR_actual)
and optimal_of > 0:
4945 if MR_actual < optimal_of:
4946 of_sq *= w_of_low_mr
4947 elif MR_actual > optimal_of:
4948 of_sq *= w_of_high_mr
4951 "thrust_penalty": float(W_THRUST * thrust_penalty_sq_term),
4952 "pc_target_penalty": float(W_PC * pc_penalty_term),
4953 "of_penalty": float(W_OF * of_sq),
4954 "cf_penalty": float(W_CF * cf_hinge),
4955 "exit_pressure_penalty": float(W_EXIT * exit_pressure_sq_term),
4956 "injector_dp_penalty": float(injector_dp_weighted),
4957 "length_penalty": float(W_LEN * length_term),
4958 "chamber_shape_penalty": float(W_CHAMBER_SHAPE * chamber_shape_term),
4959 "lstar_penalty": float(W_LSTAR * lstar_term),
4960 "isp_penalty": float(W_ISP * isp_term),
4961 "chamber_mass_penalty": float(W_MASS * mass_term),
4962 "momentum_balance_penalty": float(W_MOM * momentum_term),
4963 "resultant_tilt_penalty": float(W_TILT * tilt_term),
4964 "impingement_angle_penalty": float(W_ANGLE * angle_term),
4965 "jet_angle_asymmetry_penalty": float(W_JET_ASYM * jet_asym_term),
4966 "smd_penalty": float(W_SMD * smd_term),
4967 "tank_pressure_equal_penalty": float(W_TANK_EQUAL * tank_equal_term),
4968 "impinging_geometry_fit_penalty": float(W_IMP_GEOM * geom_fit_term),
4969 "geom_ao_af_momentum_penalty": float(geom_ao_af_weighted),
4971 obj_quality = float(sum(terms.values()))
4975 if length_violation:
4976 obj = _LAYER1_BASE_INFEAS + _LAYER1_W_INFEAS * length_term + obj_quality
4977 elif inf_residual > 0.0:
4978 obj = _LAYER1_BASE_INFEAS + _LAYER1_W_INFEAS * float(inf_residual) + obj_quality
4981 if not np.isfinite(obj):
4982 obj = _LAYER1_BASE_INFEAS
4984 t_b = float(constants.get(
"target_burn_time", 0.0)
or 0.0)
4986 "objective": float(obj),
4987 "objective_quality": obj_quality,
4989 "eval_success":
True,
4990 "infeasibility_score": float(infeasibility_score),
4991 "inf_residual": float(inf_residual),
4992 "length_violation": bool(length_violation),
4993 "F": F_actual,
"MR": MR_actual,
"Pc": Pc_actual,
"Isp": Isp_actual,
4994 "Cf": Cf_actual,
"P_exit": P_exit_actual,
4995 "thrust_error": float(thrust_error),
"of_error": float(of_error),
4996 "thrust_penalty_sq": float(thrust_penalty_sq_term),
4997 "cf_hinge": float(cf_hinge),
4998 "stability_results": stability,
"stability_state": stability_state,
4999 "stability_score": stability_score,
"state_ok": bool(state_ok),
5000 "chugging_margin": chugging_margin,
"acoustic_margin": acoustic_margin,
5001 "feed_margin": feed_margin,
"effective_min_score": effective_min_score,
5002 "effective_margin": effective_margin,
5003 "injector_dp_ratio_O": ratio_o,
"injector_dp_ratio_F": ratio_f,
5004 "momentum_ratio_R": R_val,
"resultant_tilt_deg": tilt_deg,
5005 "resultant_tilt_breakeven_deg": tilt_breakeven,
5006 "effective_smd_microns": smd_eff_um,
"impingement_angle_deg": imp_angle_deg,
5007 "chamber_D_over_Dt": dt_ratio,
"chamber_L_over_D": ld_ratio,
5008 "L_chamber": L_chamber_curr,
"L_cylindrical": L_cylindrical,
5009 "L_face_to_throat": _lens[
"face_to_throat"],
"L_engine": float(L_engine_curr),
5010 "Lstar_m": Lstar,
"Lstar_target_m": lstar_target_curr,
5011 "chamber_mass_kg": float(chamber_mass_kg),
5012 "Isp_merit_s": isp_merit,
"Isp_ref_s": isp_ref,
5013 "propellant_mass_kg": (float(target_thrust * t_b / (_G0 * isp_merit))
5014 if (t_b > 0
and np.isfinite(isp_merit)
and isp_merit > 0)
else float(
"nan")),
5015 "geom_ao_af_momentum_scale": float(geom_ao_af_scale),
5016 "geom_ao_af": geom_ao_af_ratio,
"expected_ao_af_for_R1": expected_ao_af,
5018 return out
if return_terms
else float(obj)
5588 config_obj: PintleEngineConfig,
5589 runner: PintleEngineRunner,
5590 requirements: Dict[str, Any],
5591 target_burn_time: float,
5592 tolerances: Dict[str, float],
5593 pressure_config: Dict[str, Any],
5594 update_progress: Optional[Callable[[str, float, str],
None]] =
None,
5595 log_status: Optional[Callable[[str, str],
None]] =
None,
5596 objective_callback: Optional[Callable[[int, float, float],
None]] =
None,
5597 stop_event: Optional[Any] =
None,
5598 layer1_max_iterations: Optional[int] =
None,
5599 layer1_cma_restarts: Optional[int] =
None,
5600 layer1_smoke: bool =
False,
5601) -> Tuple[PintleEngineConfig, Dict[str, Any]]:
5603 Run complete Layer 1 optimization: geometry + initial tank pressures.
5605 This function contains ALL Layer 1 optimization logic:
5606 - Setup (bounds, initial guess)
5607 - Objective function definition
5608 - Optimization loop (CMA-ES with random restarts + L-BFGS-B)
5613 Default ``max_iterations`` is 150 for robust convergence (override via ``layer1_max_iterations``
5614 for smoke tests). Design vector is **10** DOFs for pintle and **13** for impinging injectors
5615 (paired unlike doublets: shared ``n_doublets`` for LOX/fuel counts, independent jets per side).
5616 Random CMA restarts default to 3 (override via ``layer1_cma_restarts``).
5617 Set ``layer1_smoke=True`` to force pure CMA-ES (ignore ``hybrid_cma_blocks``) and,
5618 unless ``layer1_cma_restarts`` is set explicitly, use a single CMA restart for short runs.
5621 config_obj: Base engine configuration
5622 runner: Engine runner (for validation)
5623 requirements: Design requirements dict
5624 target_burn_time: Target burn time [s]
5625 tolerances: Tolerance dict (thrust, apogee)
5626 pressure_config: Pressure configuration dict
5627 update_progress: Optional progress callback (stage, progress, message)
5628 log_status: Optional status logging callback (stage, message)
5629 objective_callback: Optional callback for objective history (iteration, objective, best_objective)
5630 layer1_smoke: When True, force pure CMA-ES and (if ``layer1_cma_restarts`` is unset) a single restart.
5633 optimized_config: Optimized engine configuration
5634 results: Results dict with performance, validation, history, etc.
5638 timestamp = datetime.now().strftime(
"%Y%m%d_%H%M%S")
5640 output_logs_dir = Path(__file__).resolve().parents[3] /
"output" /
"logs"
5641 output_logs_dir.mkdir(parents=
True, exist_ok=
True)
5643 log_file_path = output_logs_dir / f
"layer1_static_{timestamp}_{os.getpid()}.log"
5646 layer1_logger = logging.getLogger(
'layer1_static')
5647 layer1_logger.setLevel(logging.INFO)
5650 layer1_logger.handlers.clear()
5653 file_handler = logging.FileHandler(log_file_path, mode=
'w', encoding=
'utf-8')
5654 file_handler.setLevel(logging.INFO)
5655 file_formatter = logging.Formatter(
'%(asctime)s - %(levelname)s - %(message)s')
5656 file_handler.setFormatter(file_formatter)
5657 layer1_logger.addHandler(file_handler)
5660 layer1_logger.propagate =
False
5662 layer1_logger.info(
"="*70)
5663 layer1_logger.info(
"Layer 1: Static Optimization")
5664 layer1_logger.info(
"="*70)
5665 layer1_logger.info(f
"Log file: {log_file_path}")
5668 if update_progress
is None:
5669 def update_progress(stage: str, progress: float, message: str):
5671 if log_status
is None:
5672 def log_status(stage: str, message: str):
5676 class OptimizationStopped(Exception):
5677 """Exception raised when optimization is stopped by user."""
5681 """Check if optimization should stop and raise exception if so."""
5682 if stop_event
is not None and stop_event.is_set():
5683 raise OptimizationStopped(
"Optimization stopped by user")
5686 target_thrust = requirements.get(
"target_thrust", 7000.0)
5687 optimal_of = requirements.get(
"optimal_of_ratio", 2.3)
5688 layer1_warnings: List[str] = []
5689 min_stability = float(requirements.get(
"min_stability_margin", _LAYER1_DEFAULT_MIN_STABILITY_MARGIN))
5691 def _resolve_Lstar_bounds_from_req_and_config() -> Tuple[float, float]:
5692 """Layer‑1 DOF ``x[1]`` is L* [m]. Prefer ``requirements``; fill missing keys from ``config_obj.design_requirements``."""
5693 mn_raw = requirements.get(
"min_Lstar",
None)
5694 mx_raw = requirements.get(
"max_Lstar",
None)
5696 dr = getattr(config_obj,
"design_requirements",
None)
5697 if mn_raw
is None and dr
is not None:
5698 mn_raw = getattr(dr,
"min_Lstar",
None)
5699 if mx_raw
is None and dr
is not None:
5700 mx_raw = getattr(dr,
"max_Lstar",
None)
5701 except (TypeError, ValueError, AttributeError):
5704 mn = float(mn_raw)
if mn_raw
is not None else 0.95
5705 except (TypeError, ValueError):
5708 mx = float(mx_raw)
if mx_raw
is not None else 1.27
5709 except (TypeError, ValueError):
5711 if not np.isfinite(mn)
or mn <= 0:
5713 if not np.isfinite(mx)
or mx <= 0:
5719 min_Lstar, max_Lstar = _resolve_Lstar_bounds_from_req_and_config()
5721 "Layer 1 optimizing L* [m] within [%.5f, %.5f] (x[1]); freeze ``design_requirements.frozen_parameters.Lstar_mm`` for a single target.",
5726 max_chamber_od = requirements.get(
"max_chamber_outer_diameter", 0.15)
5732 "Layer 1 wall (OD - inner): %.2f mm total (%.2f mm/side = metal case + ablative liner).",
5733 wall_total_m * 1e3, wall_total_m * 0.5e3,
5735 max_nozzle_exit = requirements.get(
"max_nozzle_exit_diameter", 0.101)
5736 thrust_tol = tolerances.get(
"thrust", 0.10)
5742 layer1_logger.info(f
"Site ambient (target exit pressure): {target_P_exit:.1f} Pa ({target_P_exit/101325.0:.3f} atm)")
5746 _caps = [v
for v
in (requirements.get(
"max_lox_tank_pressure_psi"),
5747 requirements.get(
"max_fuel_tank_pressure_psi"))
if v]
5748 _of_pc_pa = (float(requirements.get(
"target_chamber_pressure_psi")
or 0.0)
5749 or (
min(float(v)
for v
in _caps) / 1.3
if _caps
else 435.0)) * 6894.757
5750 _of_cea = getattr(runner,
"cea_cache",
None)
5752 config_obj, optimal_of, _of_cea
if _layer1_cea_ok(_of_cea)
else None, _of_pc_pa,
5753 float(target_P_exit))
5755 requirements = dict(requirements)
5756 requirements[
"optimal_of_ratio"] = optimal_of
5757 layer1_warnings.append(_of_msg)
5758 layer1_logger.warning(_of_msg)
5759 log_status(
"warning", _of_msg)
5761 layer1_logger.info(f
"Target thrust: {target_thrust:.1f} N")
5762 layer1_logger.info(f
"Target O/F ratio: {optimal_of:.2f}")
5763 layer1_logger.info(f
"Min stability margin: {min_stability:.2f}")
5764 layer1_logger.info(f
"Target burn time: {target_burn_time:.2f} s")
5775 max_iterations = 1042
5776 if layer1_max_iterations
is not None:
5777 max_iterations =
max(1, int(layer1_max_iterations))
5779 layer1_logger.info(f
"Max iterations: {max_iterations}")
5780 layer1_logger.info(
"")
5783 max_lox_P_psi = float(pressure_config.get(
"max_lox_pressure_psi", 700))
5784 max_fuel_P_psi = float(pressure_config.get(
"max_fuel_pressure_psi", 850))
5788 _dr_cap = getattr(config_obj,
"design_requirements",
None)
5789 if _dr_cap
is not None:
5790 _lox_cap = getattr(_dr_cap,
"max_lox_tank_pressure_psi",
None)
5791 if _lox_cap
is not None and float(_lox_cap) > 0.0:
5792 max_lox_P_psi =
min(max_lox_P_psi, float(_lox_cap))
5793 _f_cap = getattr(_dr_cap,
"max_fuel_tank_pressure_psi",
None)
5794 if _f_cap
is not None and float(_f_cap) > 0.0:
5795 max_fuel_P_psi =
min(max_fuel_P_psi, float(_f_cap))
5796 except (TypeError, ValueError, AttributeError):
5799 "Layer 1 tank stagnation caps (psi): max LOX=%.1f, max fuel=%.1f",
5808 _lbfgs_2 = requirements.get(
"layer1_lbfgs_second_pass")
5809 layer1_lbfgs_second_pass =
True if _lbfgs_2
is None else bool(_lbfgs_2)
5811 update_progress(
"Layer 1: Setup", 0.01,
"Initializing Layer 1 optimization...")
5814 config_base = copy.deepcopy(config_obj)
5815 if hasattr(config_base,
'injector')
and config_base.injector.type ==
"pintle":
5816 if hasattr(config_base.injector.geometry,
'lox'):
5817 config_base.injector.geometry.lox.theta_orifice = 90.0
5819 l1_injector_type = getattr(getattr(config_base,
"injector",
None),
"type",
"pintle")
5820 if l1_injector_type
not in (
"pintle",
"impinging"):
5822 "Layer 1 static optimization supports injector types 'pintle' and 'impinging' only; "
5823 f
"got {l1_injector_type!r}."
5825 idx_P_O = 11
if l1_injector_type ==
"impinging" else 8
5826 idx_P_F = 12
if l1_injector_type ==
"impinging" else 9
5827 layer1_integer_dims = [4]
if l1_injector_type ==
"impinging" else [6]
5839 requirements,
"layer1_integer_jet_angles",
True
5841 layer1_integer_dims = layer1_integer_dims + [6, 9]
5843 "Layer 1: impingement angles snapped to whole degrees (machining limit); "
5844 "set layer1_integer_jet_angles=false to search fractional angles."
5848 if hasattr(config_base,
'combustion')
and hasattr(config_base.combustion,
'efficiency'):
5849 config_base.combustion.efficiency.use_turbulence_coupling =
True
5855 _pc_tgt_seed = requirements.get(
"target_chamber_pressure_psi")
5856 if _pc_tgt_seed
is not None and float(_pc_tgt_seed) > 0:
5857 Pc_est_psi = float(_pc_tgt_seed)
5860 Pc_est = Pc_est_psi * psi_to_Pa
5862 A_throat_init = target_thrust / (Cf_est * Pc_est)
if Pc_est > 0
else 0.001
5863 A_throat_init =
max(float(A_throat_init), 5e-5)
5866 min_A_throat_safe = 5e-5
5881 _fmin = requirements.get(
"layer1_stagnation_pressure_frac_min")
5882 _fmax = requirements.get(
"layer1_stagnation_pressure_frac_max")
5883 min_P_ratio = float(_fmin)
if _fmin
is not None else 0.35
5884 max_P_ratio = float(_fmax)
if _fmax
is not None else 1.00
5885 if max_P_ratio <= min_P_ratio:
5886 max_P_ratio = min_P_ratio + 1e-3
5889 if exp_ratio_hi <= exp_ratio_lo:
5890 exp_ratio_hi = exp_ratio_lo + 1e-3
5898 _cache_eps = getattr(getattr(config_obj.combustion,
"cea",
None),
"eps_range",
None)
5899 if _cache_eps
and len(_cache_eps) == 2:
5900 _ce_lo, _ce_hi = float(_cache_eps[0]), float(_cache_eps[1])
5901 _clamped_lo, _clamped_hi =
max(exp_ratio_lo, _ce_lo),
min(exp_ratio_hi, _ce_hi)
5902 if (_clamped_lo, _clamped_hi) != (exp_ratio_lo, exp_ratio_hi):
5903 layer1_logger.warning(
5904 "Expansion-ratio search [%g, %g] exceeds the CEA cache grid [%g, %g]; "
5905 "clamping to [%g, %g] to avoid extrapolating the table.",
5906 exp_ratio_lo, exp_ratio_hi, _ce_lo, _ce_hi, _clamped_lo, _clamped_hi,
5908 exp_ratio_lo, exp_ratio_hi = _clamped_lo, _clamped_hi
5909 if exp_ratio_hi <= exp_ratio_lo:
5910 exp_ratio_hi = exp_ratio_lo + 1e-3
5911 except (AttributeError, TypeError, ValueError):
5916 _cf_hi_for_floor =
_requirement_float(requirements,
"layer1_cf_upper_bound_for_throat_floor", 1.8)
5917 _pc_floor_factor =
_requirement_float(requirements,
"layer1_pc_fraction_for_throat_floor", 0.75)
5918 _pc_cap_psi =
min(max_lox_P_psi * max_P_ratio, max_fuel_P_psi * max_P_ratio)
5919 _pc_cap_pa = _pc_cap_psi * psi_to_Pa *
max(0.1, _pc_floor_factor)
5920 if target_thrust > 0
and _pc_cap_pa > 0
and _cf_hi_for_floor > 0:
5921 _a_req = target_thrust / (_cf_hi_for_floor * _pc_cap_pa)
5922 if np.isfinite(_a_req)
and _a_req > 0:
5924 min_A_throat_safe =
max(min_A_throat_safe, 0.7 * float(_a_req))
5926 "A_throat lower bound raised by thrust-feasibility floor: %.2f mm^2 "
5927 "(target %.0f N, Pc_cap≈%.1f psi, Cf_hi=%.2f)",
5928 min_A_throat_safe * 1e6,
5937 max_chamber_od, wall_total_m, min_A_throat_safe,
5939 if not (min_outer_diameter < max_chamber_od):
5940 layer1_logger.warning(
5941 "Chamber OD cap %.1f mm cannot hold the %.1f mm wall plus a bore at the D/Dt floor "
5942 "for the smallest reachable throat; the chamber-shape term will report it.",
5943 max_chamber_od * 1e3, wall_total_m * 1e3)
5944 min_outer_diameter = max_chamber_od * 0.5
5947 A_throat_hi =
max(4.0e-3, float(target_thrust) / (
5948 1.3 *
max(min_P_ratio, 0.05) *
min(max_lox_P_psi, max_fuel_P_psi) * psi_to_Pa))
5950 if l1_injector_type ==
"impinging":
5951 D_inner_bounds = impinging_chamber_inner_diameter_for_bounds(
5953 max_chamber_outer_diameter_m=max_chamber_od,
5954 wall_thickness_m=wall_total_m,
5956 n_hi_int = impinging_n_elements_hi_int(D_inner_bounds)
5957 _nd_cap = requirements.get(
"layer1_impinging_n_doublets_max")
5958 if _nd_cap
is not None:
5962 n_hi_int =
min(n_hi_int, _c)
5963 except (TypeError, ValueError):
5966 n_hi_upper = float(n_hi_int) + 0.499
5967 d_jet_hi = impinging_d_jet_upper_bound_m(D_inner_bounds)
5968 spacing_hi = impinging_spacing_upper_bound_m(D_inner_bounds)
5982 _ang_lo_req = requirements.get(
"layer1_impinging_angle_deg_min")
5983 _ang_hi_req = requirements.get(
"layer1_impinging_angle_deg_max")
5984 _incl_lo = float(_ang_lo_req)
if _ang_lo_req
is not None else 30.0
5985 _incl_hi = float(_ang_hi_req)
if _ang_hi_req
is not None else 170.0
5986 _asym =
_requirement_float(requirements,
"layer1_impinging_jet_angle_max_asym_deg", 28.0)
5987 imp_angle_lo_bound = float(np.clip(0.5 * (_incl_lo - _asym), 5.0, 89.0))
5988 imp_angle_hi_bound = float(np.clip(0.5 * (_incl_hi + _asym), 6.0, 89.0))
5994 imp_angle_lo_bound =
min(imp_angle_lo_bound, 15.0)
5995 imp_angle_hi_bound =
max(imp_angle_hi_bound, 85.0)
6002 _jet_lo = requirements.get(
"layer1_impinging_jet_angle_min_deg")
6003 if _jet_lo
is not None:
6005 _jet_lo = float(_jet_lo)
6006 if np.isfinite(_jet_lo)
and 0.0 < _jet_lo < imp_angle_hi_bound:
6007 imp_angle_lo_bound =
max(imp_angle_lo_bound, _jet_lo)
6008 except (TypeError, ValueError):
6017 _inc_req = requirements.get(
"layer1_injector_min_face_incidence_deg")
6019 if _inc_req
is not None and _pl.get(
"face") ==
"contoured" and _pl.get(
"back") ==
"channels":
6024 "Contoured injector face: layer1_injector_min_face_incidence_deg (%s) not applied",
6027 if _inc_req
is not None:
6029 _inc_req = float(_inc_req)
6030 if np.isfinite(_inc_req)
and _inc_req > 0.0:
6031 _th_ceiling = 90.0 - _inc_req
6032 if _th_ceiling <= imp_angle_lo_bound:
6034 f
"layer1_injector_min_face_incidence_deg = {_inc_req:g} caps the jet "
6035 f
"angle at {_th_ceiling:g} deg, at or below the "
6036 f
"layer1_impinging_jet_angle_min_deg floor of {imp_angle_lo_bound:g}. "
6037 f
"No jet angle satisfies both."
6039 imp_angle_hi_bound =
min(imp_angle_hi_bound, _th_ceiling)
6041 "Face-incidence floor %g deg -> jet angle capped at %.1f deg "
6042 "(drill meets the face at >= %g deg)",
6043 _inc_req, imp_angle_hi_bound, _inc_req)
6044 except (TypeError, ValueError)
as _e:
6045 if isinstance(_e, ValueError)
and "min_face_incidence" in str(_e):
6047 if imp_angle_hi_bound <= imp_angle_lo_bound:
6048 imp_angle_hi_bound =
min(89.0, imp_angle_lo_bound + 1.0)
6050 "Impingement-angle bounds: included [%g, %g] deg (asym <= %g) -> per-jet [%.1f, %.1f]",
6051 _incl_lo, _incl_hi, _asym, imp_angle_lo_bound, imp_angle_hi_bound,
6061 if np.isfinite(_cc_req)
and _cc_req > 0.0
and n_hi_int >= 1:
6062 _spacing_lo =
max(_spacing_lo, float(np.pi) * float(_cc_req) / float(n_hi_int))
6063 if _spacing_lo >= spacing_hi:
6065 f
"layer1_injector_center_clear_dia_m = {_cc_req*1000:.1f} mm needs a hole "
6066 f
"pitch of at least {_spacing_lo*1000:.2f} mm at n = {n_hi_int}, but the "
6067 f
"pitch ceiling for a {D_inner_bounds*1000:.1f} mm bore is "
6068 f
"{spacing_hi*1000:.2f} mm. Reduce the reserved centre circle, or raise "
6069 f
"layer1_impinging_n_doublets_max."
6072 "Centre clearance %.1f mm reserved -> hole-pitch floor %.2f mm (was 3.00) at n=%d",
6073 _cc_req * 1000.0, _spacing_lo * 1000.0, n_hi_int,
6076 f
"Impinging injector bounds from chamber bore D_inner≈{D_inner_bounds*1000:.2f} mm: "
6077 f
"n_doublets ≤ {n_hi_int} (hard int cap via search box), "
6078 f
"d_jet ≤ {d_jet_hi*1000:.2f} mm, spacing ≤ {spacing_hi*1000:.2f} mm"
6081 (min_A_throat_safe, A_throat_hi),
6082 (min_Lstar, max_Lstar),
6083 (exp_ratio_lo, exp_ratio_hi),
6084 (min_outer_diameter, max_chamber_od),
6087 (imp_angle_lo_bound, imp_angle_hi_bound),
6088 (_spacing_lo, spacing_hi),
6090 (imp_angle_lo_bound, imp_angle_hi_bound),
6091 (_spacing_lo, spacing_hi),
6092 (max_lox_P_psi * min_P_ratio, max_lox_P_psi * max_P_ratio),
6093 (max_fuel_P_psi * min_P_ratio, max_fuel_P_psi * max_P_ratio),
6097 max_d_orifice = 0.003
6098 max_LOX_area = max_n_orifices * np.pi * (max_d_orifice / 2) ** 2
6099 max_d_pintle = 0.040
6101 R_inner_max = max_d_pintle / 2
6102 R_outer_max = R_inner_max + max_h_gap
6103 max_fuel_area = np.pi * (R_outer_max ** 2 - R_inner_max ** 2)
6104 max_injector_area =
max(max_LOX_area, max_fuel_area)
6107 (min_A_throat_safe, A_throat_hi),
6108 (min_Lstar, max_Lstar),
6109 (exp_ratio_lo, exp_ratio_hi),
6110 (min_outer_diameter, max_chamber_od),
6115 (max_lox_P_psi * min_P_ratio, max_lox_P_psi * max_P_ratio),
6116 (max_fuel_P_psi * min_P_ratio, max_fuel_P_psi * max_P_ratio),
6121 _p_o_min = requirements.get(
"layer1_P_O_start_psi_min")
6122 _p_o_max = requirements.get(
"layer1_P_O_start_psi_max")
6123 if _p_o_min
is not None and _p_o_max
is not None:
6124 lo_p, hi_p = float(_p_o_min),
min(float(_p_o_max), max_lox_P_psi)
6125 lo_p =
min(lo_p, hi_p)
6127 bounds[idx_P_O] = (lo_p, hi_p)
6128 _p_f_min = requirements.get(
"layer1_P_F_start_psi_min")
6129 _p_f_max = requirements.get(
"layer1_P_F_start_psi_max")
6130 if _p_f_min
is not None and _p_f_max
is not None:
6131 lo_pf, hi_pf = float(_p_f_min),
min(float(_p_f_max), max_fuel_P_psi)
6132 lo_pf =
min(lo_pf, hi_pf)
6134 bounds[idx_P_F] = (lo_pf, hi_pf)
6137 _o_lo, _o_hi = bounds[idx_P_O]
6138 bounds[idx_P_O] = (float(
min(_o_lo, max_lox_P_psi)), float(
min(_o_hi, max_lox_P_psi)))
6139 if bounds[idx_P_O][0] > bounds[idx_P_O][1]:
6140 bounds[idx_P_O] = (bounds[idx_P_O][1], bounds[idx_P_O][0])
6141 _f_lo, _f_hi = bounds[idx_P_F]
6142 bounds[idx_P_F] = (float(
min(_f_lo, max_fuel_P_psi)), float(
min(_f_hi, max_fuel_P_psi)))
6143 if bounds[idx_P_F][0] > bounds[idx_P_F][1]:
6144 bounds[idx_P_F] = (bounds[idx_P_F][1], bounds[idx_P_F][0])
6149 has_valid_start_geom =
False
6154 def _get_val(obj, attr, default=None):
6155 val = getattr(obj, attr, default)
6156 return float(val)
if val
is not None and val > 0
else None
6160 cg_in = getattr(config_obj,
'chamber_geometry',
None)
6161 c_in = getattr(config_obj,
'chamber',
None)
6162 n_in = getattr(config_obj,
'nozzle',
None)
6164 A_throat_in = _get_val(cg_in,
'A_throat')
or _get_val(c_in,
'A_throat')
or _get_val(n_in,
'A_throat')
6165 Lstar_in = _get_val(cg_in,
'Lstar')
or _get_val(c_in,
'Lstar')
6166 eps_in = _get_val(cg_in,
'expansion_ratio')
or _get_val(n_in,
'expansion_ratio')
6168 D_inner_in = _get_val(cg_in,
'chamber_diameter')
or _get_val(c_in,
'chamber_inner_diameter')
6171 inj_in = getattr(config_obj,
'injector',
None)
6172 inj_geom = getattr(inj_in,
'geometry',
None)
6173 is_pintle = (getattr(inj_in,
'type',
'') ==
'pintle')
6174 is_impinging = (getattr(inj_in,
'type',
'') ==
'impinging')
6176 d_pintle_in = h_gap_in = n_orifices_in = d_orifice_in =
None
6177 n_O_in = d_jet_O_in = ang_O_in = sp_O_in =
None
6178 n_F_in = d_jet_F_in = ang_F_in = sp_F_in =
None
6180 if is_pintle
and inj_geom:
6181 fuel_in = getattr(inj_geom,
'fuel',
None)
6182 lox_in = getattr(inj_geom,
'lox',
None)
6184 d_pintle_in = _get_val(fuel_in,
'd_pintle_tip')
6185 h_gap_in = _get_val(fuel_in,
'h_gap')
6186 n_orifices_in = _get_val(lox_in,
'n_orifices')
6187 d_orifice_in = _get_val(lox_in,
'd_orifice')
6188 elif is_impinging
and inj_geom:
6189 ox_in = getattr(inj_geom,
'oxidizer',
None)
6190 fu_in = getattr(inj_geom,
'fuel',
None)
6192 n_O_in = getattr(ox_in,
'n_elements',
None)
6193 n_F_in = getattr(fu_in,
'n_elements',
None)
6194 n_doublets_in =
None
6195 if n_O_in
is not None and float(n_O_in) > 0:
6197 if n_F_in
is not None and float(n_F_in) > 0:
6200 layer1_logger.warning(
6201 "Impinging YAML n_elements_O (%s) != n_elements_F (%s); paired doublets "
6202 "use one n_doublets — Layer-1 x0 uses min(O,F)=%s.",
6207 n_doublets_in =
min(no_i, nf_i)
6209 n_doublets_in = no_i
6210 elif n_F_in
is not None and float(n_F_in) > 0:
6211 n_doublets_in = int(n_F_in)
6212 d_jet_O_in = _get_val(ox_in,
'd_jet')
6213 ang_O_in = _get_val(ox_in,
'impingement_angle')
6214 sp_O_in = _get_val(ox_in,
'spacing')
6215 d_jet_F_in = _get_val(fu_in,
'd_jet')
6216 ang_F_in = _get_val(fu_in,
'impingement_angle')
6217 sp_F_in = _get_val(fu_in,
'spacing')
6256 A_throat_init = A_throat_in
6257 Lstar_init = Lstar_in
6259 outer_diameter_init = D_inner_in + wall_total_m
6261 default_d_pintle = d_pintle_in
6262 default_h_gap = h_gap_in
6263 default_n_orifices = int(n_orifices_in)
6264 default_d_orifice = d_orifice_in
6266 has_valid_start_geom =
True
6267 layer1_logger.info(
"Using values from input configuration as initial guess (x0).")
6270 A_throat_init = A_throat_in
6271 Lstar_init = Lstar_in
6273 outer_diameter_init = D_inner_in + wall_total_m
6275 default_n_doublets = int(n_doublets_in)
6276 default_d_jet_O = float(d_jet_O_in)
6277 default_d_jet_F = float(d_jet_F_in)
6278 default_ang_O = float(ang_O_in)
6279 default_ang_F = float(ang_F_in)
6280 default_sp_O = float(sp_O_in)
6281 default_sp_F = float(sp_F_in)
6283 has_valid_start_geom =
True
6284 layer1_logger.info(
"Using impinging injector geometry from input configuration as initial guess (x0).")
6287 layer1_logger.info(
"Input configuration incomplete/invalid for x0. Using heuristic defaults.")
6289 except Exception
as e:
6290 layer1_logger.warning(f
"Failed to extract x0 from config: {e}. Using defaults.")
6291 has_valid_start_geom =
False
6293 if not has_valid_start_geom:
6294 if l1_injector_type ==
"impinging":
6295 default_n_doublets = 12
6296 default_d_jet_O = 0.002
6297 default_d_jet_F = 0.002
6298 default_ang_O = 45.0
6299 default_ang_F = 45.0
6300 default_sp_O = 0.012
6301 default_sp_F = 0.012
6302 A_lox_est = default_n_doublets * np.pi * (default_d_jet_O / 2.0) ** 2
6303 A_fuel_est = default_n_doublets * np.pi * (default_d_jet_F / 2.0) ** 2
6304 max_injector_area_est =
max(A_lox_est, A_fuel_est)
6305 A_throat_min_safe = max_injector_area_est * 1.15
6306 A_throat_init =
max(A_throat_init, A_throat_min_safe)
6307 Lstar_init = (min_Lstar + max_Lstar) / 2
6309 outer_diameter_init = np.clip(max_chamber_od * 0.55, min_outer_diameter, max_chamber_od)
6311 default_n_orifices = 14
6312 default_d_orifice = 0.002
6313 default_d_pintle = 0.016
6314 default_h_gap = 0.0006
6316 A_lox_est = default_n_orifices * np.pi * (default_d_orifice / 2) ** 2
6317 R_inner_est = default_d_pintle / 2
6318 R_outer_est = R_inner_est + default_h_gap
6319 A_fuel_est = np.pi * (R_outer_est ** 2 - R_inner_est ** 2)
6320 max_injector_area_est =
max(A_lox_est, A_fuel_est)
6322 A_throat_min_safe = max_injector_area_est * 1.1
6323 A_throat_init =
max(A_throat_init, A_throat_min_safe)
6325 Lstar_init = (min_Lstar + max_Lstar) / 2
6327 outer_diameter_init = np.clip(max_chamber_od * 0.55, min_outer_diameter, max_chamber_od)
6330 A_throat_init = float(np.clip(A_throat_init, bounds[0][0], bounds[0][1]))
6332 P_O_start_init = max_lox_P_psi * 0.80
6333 P_F_start_init = max_fuel_P_psi * 0.80
6334 P_O_start_init = float(np.clip(P_O_start_init, bounds[idx_P_O][0], bounds[idx_P_O][1]))
6335 P_F_start_init = float(np.clip(P_F_start_init, bounds[idx_P_F][0], bounds[idx_P_F][1]))
6337 if l1_injector_type ==
"impinging":
6340 d_jet_O=default_d_jet_O,
6341 d_jet_F=default_d_jet_F,
6342 n_doublets=int(default_n_doublets),
6343 optimal_of=float(optimal_of),
6345 logger=layer1_logger,
6346 also_config=config_base,
6349 if l1_injector_type ==
"impinging":
6352 Lstar_init
if 'Lstar_init' in locals()
else (min_Lstar + max_Lstar) / 2,
6353 eps_init
if 'eps_init' in locals()
else 8.0,
6354 outer_diameter_init,
6355 float(default_n_doublets),
6368 Lstar_init
if 'Lstar_init' in locals()
else (min_Lstar + max_Lstar) / 2,
6369 eps_init
if 'eps_init' in locals()
else 8.0,
6370 outer_diameter_init,
6390 def _frozen_param_nonempty(fp: Any) -> bool:
6391 return isinstance(fp, dict)
and any(v
is not None for v
in fp.values())
6393 _fp_req = requirements.get(
"frozen_parameters")
6394 _fp_cfg: Dict[str, Any] = {}
6395 dr0 = getattr(config_base,
"design_requirements",
None)
6396 if dr0
is not None and dr0.frozen_parameters
is not None:
6397 _fp_cfg = dr0.frozen_parameters.model_dump(exclude_none=
True)
6399 if not _frozen_param_nonempty(_fp_req):
6400 requirements = dict(requirements)
6401 requirements[
"frozen_parameters"] = dict(_fp_cfg)
6402 elif isinstance(_fp_req, dict):
6403 merged_fp = dict(_fp_cfg)
6404 for _k, _v
in _fp_req.items():
6407 requirements = dict(requirements)
6408 requirements[
"frozen_parameters"] = merged_fp
6410 frozen = requirements.get(
"frozen_parameters", {})
or {}
6411 frozen_param_names = []
6433 if l1_injector_type ==
"impinging":
6435 (
"A_throat_mm2", 1e-6),
6437 (
"expansion_ratio", 1.0),
6438 (
"D_chamber_outer_mm", 1e-3),
6439 (
"n_doublets", 1.0),
6440 (
"d_jet_O_mm", 1e-3),
6441 (
"impingement_angle_O_deg", 1.0),
6442 (
"spacing_O_mm", 1e-3),
6443 (
"d_jet_F_mm", 1e-3),
6444 (
"impingement_angle_F_deg", 1.0),
6445 (
"spacing_F_mm", 1e-3),
6446 (
"P_O_start_psi", 1.0),
6447 (
"P_F_start_psi", 1.0),
6451 (
"A_throat_mm2", 1e-6),
6453 (
"expansion_ratio", 1.0),
6454 (
"D_chamber_outer_mm", 1e-3),
6455 (
"d_pintle_tip_mm", 1e-3),
6457 (
"n_orifices", 1.0),
6458 (
"d_orifice_mm", 1e-3),
6459 (
"P_O_start_psi", 1.0),
6460 (
"P_F_start_psi", 1.0),
6463 int_keys = {
"n_orifices",
"n_doublets"}
6465 for idx, (key, conversion)
in enumerate(frozen_mapping):
6466 frozen_val = frozen.get(key)
6467 if frozen_val
is not None:
6468 val_si = frozen_val * conversion
6470 val_si = int(round(frozen_val))
6471 eps = 1e-12
if key
not in int_keys
else 0.01
6472 bounds[idx] = (val_si - eps, val_si + eps)
6474 frozen_param_names.append(f
"{key}={frozen_val}")
6476 if frozen_param_names:
6477 layer1_logger.info(f
"Frozen parameters: {', '.join(frozen_param_names)}")
6478 layer1_logger.info(
"These values will be fixed during optimization.")
6482 _atm_lo = requirements.get(
"layer1_A_throat_mm2_min")
6483 _atm_hi = requirements.get(
"layer1_A_throat_mm2_max")
6484 if _atm_lo
is not None and _atm_hi
is not None:
6485 lo_m2 = float(_atm_lo) * 1e-6
6486 hi_m2 = float(_atm_hi) * 1e-6
6489 "layer1_A_throat_mm2_max must exceed layer1_A_throat_mm2_min "
6490 f
"(got min={_atm_lo}, max={_atm_hi})."
6492 b0_lo =
max(float(bounds[0][0]), lo_m2)
6493 b0_hi =
min(float(bounds[0][1]), hi_m2)
6496 "layer1 A_throat band does not intersect default bounds "
6497 f
"(band mm² [{_atm_lo}, {_atm_hi}], intersect [{b0_lo}, {b0_hi}] m²)."
6499 bounds[0] = (b0_lo, b0_hi)
6501 "Layer 1 A_throat interval from requirements "
6502 "layer1_A_throat_mm2_min/max: [%.6f, %.6f] mm²",
6515 _pc_tgt_psi = requirements.get(
"target_chamber_pressure_psi")
6516 if _pc_tgt_psi
is not None and float(_pc_tgt_psi) > 0:
6517 _pc_t = float(_pc_tgt_psi)
6518 for _pidx, _cap
in ((idx_P_O, max_lox_P_psi), (idx_P_F, max_fuel_P_psi)):
6519 _lo =
max(_pc_t * 1.12, 50.0)
6520 _hi =
min(_pc_t * 1.75, float(_cap))
6524 _lo =
max(0.5 * _hi,
min(_pc_t, _hi - 1.0))
6525 bounds[_pidx] = (_lo, _hi)
6526 x0[_pidx] = 0.5 * (_lo + _hi)
6528 "Layer 1 tank-pressure window re-anchored to Pc target %.0f psi: tanks ∈ [%.0f, %.0f] psi "
6529 "(was cap-relative). This is what lets Pc actually reach the target.",
6530 _pc_t, bounds[idx_P_O][0], bounds[idx_P_O][1],
6547 requirements,
"layer1_derive_expansion_ratio",
True)
6549 requirements,
"layer1_derive_throat_from_thrust",
True)
6556 requirements,
"layer1_derive_max_iters", float(_LAYER1_DERIVE_MAX_ITERS_DEFAULT))))
6558 requirements,
"layer1_derive_thrust_tol_rel", 1.0e-3))
6560 if derive_throat_from_thrust
and not (target_thrust > 0):
6561 layer1_logger.warning(
6562 "Layer 1: A_throat derivation requested but no positive thrust target; "
6563 "searching A_throat instead."
6565 derive_throat_from_thrust =
False
6567 derive_At_range = (float(bounds[0][0]), float(bounds[0][1]))
6568 derive_eps_range = (float(bounds[2][0]), float(bounds[2][1]))
6569 derive_dF_range = ((float(bounds[8][0]), float(bounds[8][1]))
if len(bounds) > 8
else (1e-5, 1.0))
6570 _tanks_locked_flag =
_requirement_bool(requirements,
"layer1_lock_tank_pressures",
6578 requirements,
"layer1_derive_tank_from_dp_ratio",
False)
6581 dp_ratio_target =
_requirement_float(requirements,
"layer1_dp_ratio_target", 0.5 * (_dp_lo + _dp_hi))
6582 lock_P_T_range = (float(bounds[idx_P_O][0]), float(bounds[idx_P_O][1]))
6588 derive_fuel_jet_from_of = (
6589 l1_injector_type ==
"impinging"
6595 derive_spacing_F_range = (
6596 (float(bounds[10][0]), float(bounds[10][1]))
if len(bounds) > 10
else (0.0, float(
"inf"))
6598 lock_P_F_range = (float(bounds[idx_P_F][0]), float(bounds[idx_P_F][1]))
6609 _l1_fixed_vars: Dict[int, float] = {}
6612 _pf_seed = float(np.clip(x0[idx_P_O], *lock_P_F_range))
6613 x0[idx_P_F] = _pf_seed
6614 _l1_fixed_vars[idx_P_F] = _pf_seed
6615 _pinned_msgs.append(
"fuel tank pressure (locked to the LOX tank)")
6616 if derive_throat_from_thrust:
6617 _at_seed = float(np.clip(x0[0], *derive_At_range))
6619 _l1_fixed_vars[0] = _at_seed
6620 _pinned_msgs.append(
6621 "A_throat (solved from thrust target %.0f N, search range %.1f-%.1f mm²)"
6622 % (target_thrust, derive_At_range[0] * 1e6, derive_At_range[1] * 1e6)
6624 if derive_expansion_ratio:
6625 _eps_seed = float(np.clip(x0[2], *derive_eps_range))
6627 _l1_fixed_vars[2] = _eps_seed
6628 _pinned_msgs.append(
6629 "expansion_ratio (solved for Pe = %.2f psi, clamp %.2f-%.2f)"
6630 % (target_P_exit / 6894.76, derive_eps_range[0], derive_eps_range[1])
6632 if (l1_injector_type ==
"impinging" and len(bounds) > 10
6633 and _requirement_bool(requirements,
"layer1_derive_impingement_spacing",
True)):
6634 _spf_seed = float(np.clip(x0[10], *derive_spacing_F_range))
6636 _l1_fixed_vars[10] = _spf_seed
6637 if derive_fuel_jet_from_of
and len(bounds) > 8:
6638 _l1_fixed_vars[8] = float(np.clip(x0[8], *derive_dF_range))
6639 _pinned_msgs.append(
"fuel jet diameter (solved from the O/F target)")
6640 if derive_tank_from_dp:
6641 _l1_fixed_vars[idx_P_O] = float(np.clip(x0[idx_P_O], *lock_P_T_range))
6642 _pinned_msgs.append(
"tank pressure (solved from dP/Pc = %.2f)" % dp_ratio_target)
6643 _pinned_msgs.append(
6644 "fuel ring pitch (solved so the jets meet at %.2f jet diameters)"
6645 %
_requirement_float(requirements,
"layer1_impingement_Ld_target", IMPINGEMENT_LD_TARGET_DEFAULT)
6649 "Layer 1 dimensionality: %d of %d DOFs are solved, not searched -- %s",
6650 len(_pinned_msgs), len(bounds),
"; ".join(_pinned_msgs),
6654 requirements,
"layer1_chamber_od_increment_in", 0.0)
6655 _od_snap_target = str(
6656 requirements.get(
"layer1_chamber_od_snap_target")
or "outer").strip().lower()
6657 if _od_snap_target
not in (
"outer",
"bore"):
6658 layer1_logger.warning(
6659 "Layer 1: unknown layer1_chamber_od_snap_target=%r; using 'outer'.",
6661 _od_snap_target =
"outer"
6664 _derive_constants = {
6665 "layer1_derive_expansion_ratio": derive_expansion_ratio,
6666 "layer1_derive_throat_from_thrust": derive_throat_from_thrust,
6667 "layer1_derive_max_iters": derive_max_iters,
6668 "layer1_derive_thrust_tol_rel": derive_thrust_tol,
6669 "derive_eps_min": derive_eps_range[0],
6670 "derive_eps_max": derive_eps_range[1],
6671 "derive_At_min": derive_At_range[0],
6672 "derive_At_max": derive_At_range[1],
6673 "derive_dF_min": derive_dF_range[0],
6674 "derive_dF_max": derive_dF_range[1],
6675 "layer1_derive_fuel_jet_from_of": derive_fuel_jet_from_of,
6676 "layer1_derive_tank_from_dp_ratio": derive_tank_from_dp,
6677 "layer1_dp_ratio_target": dp_ratio_target,
6678 "lock_P_T_min": lock_P_T_range[0],
6679 "lock_P_T_max": lock_P_T_range[1],
6680 "optimal_of": optimal_of,
6681 "P_ambient": target_P_exit,
6682 "target_thrust": target_thrust,
6683 "max_nozzle_exit": max_nozzle_exit,
6684 "TOTAL_WALL_THICKNESS_M": wall_total_m,
6685 "layer1_chamber_od_increment_in": _od_increment_in,
6686 "layer1_chamber_od_snap_target": _od_snap_target,
6687 "injector_type": l1_injector_type,
6691 requirements,
"layer1_lock_tank_pressures",
6693 "lock_P_F_min": lock_P_F_range[0],
6694 "lock_P_F_max": lock_P_F_range[1],
6696 requirements,
"layer1_derive_impingement_spacing",
True),
6698 requirements,
"layer1_impingement_Ld_target", IMPINGEMENT_LD_TARGET_DEFAULT),
6700 requirements,
"layer1_ring_order_fuel_outboard",
True),
6701 "derive_spacing_F_min": derive_spacing_F_range[0],
6702 "derive_spacing_F_max": derive_spacing_F_range[1],
6706 for i, (lo, hi)
in enumerate(bounds):
6707 x0[i] = np.clip(x0[i], lo, hi)
6709 update_progress(
"Layer 1: Setup", 0.05,
"Creating apply_x_to_config function...")
6711 bounds, max_chamber_od, max_nozzle_exit, l1_injector_type,
6712 chamber_od_increment_in=_od_increment_in,
6713 chamber_od_snap_target=_od_snap_target,
6714 derive_constants=_derive_constants,
6715 derive_spacing_F_range=derive_spacing_F_range,
6719 lower_bounds = np.array([b[0]
for b
in bounds], dtype=float)
6720 upper_bounds = np.array([b[1]
for b
in bounds], dtype=float)
6729 lower_bounds_eff = lower_bounds.copy()
6730 upper_bounds_eff = upper_bounds.copy()
6732 (0, derive_At_range
if derive_throat_from_thrust
else None),
6733 (2, derive_eps_range
if derive_expansion_ratio
else None),
6734 (10, derive_spacing_F_range
6735 if (l1_injector_type ==
"impinging" and len(bounds) > 10
6738 (idx_P_F, lock_P_F_range
6743 if _rng
is not None and _di < len(lower_bounds_eff):
6744 lower_bounds_eff[_di] = float(_rng[0])
6745 upper_bounds_eff[_di] = float(_rng[1])
6746 span = np.maximum(upper_bounds - lower_bounds, 1e-9)
6749 _rs = requirements.get(
"layer1_random_seed")
6750 layer1_seed_base = int(_rs)
if _rs
is not None else int(np.random.SeedSequence().entropy % (2**31))
6751 rng = np.random.default_rng(layer1_seed_base)
6753 "Layer 1 RNG/CMA reproducibility: layer1_random_seed base=%s "
6754 "(CMA seed per restart = base + restart_idx × 1_000_003)",
6761 _ren = requirements.get(
"report_every_n")
6762 report_every_n = int(_ren)
if _ren
is not None else 1
6763 report_every_n =
max(1, report_every_n)
6768 "function_evaluations": 0,
6769 "best_objective": float(
'inf'),
6771 "best_config":
None,
6772 "best_config_x":
None,
6773 "best_lox_end_ratio":
None,
6774 "best_fuel_end_ratio":
None,
6775 "best_pressures":
None,
6776 "last_eval_config":
None,
6777 "last_eval_config_x":
None,
6778 "best_validation_tank_pa":
None,
6779 "last_eval_validation_tank_pa":
None,
6780 "last_good_eval_bundle":
None,
6781 "best_results_for_validation":
None,
6783 "objective_satisfied":
False,
6784 "satisfied_obj": float(
'inf'),
6785 "satisfied_eval_count": 0,
6786 "consecutive_failures": 0,
6787 "last_valid_obj": float(
'inf'),
6789 "objective_buffer": [],
6790 "stop_optimization":
False,
6791 "force_maxfun_1":
False,
6792 "last_best_eval": 0,
6793 "valley_escape_tier": 0,
6794 "cooldown_until": 0,
6800 "layer1_momentum_log_deadband_rel",
6801 _LAYER1_DEFAULT_MOMENTUM_LOG_DEADBAND_REL,
6803 _irm_lo = requirements.get(
"impinging_momentum_R_min")
6804 _irm_hi = requirements.get(
"impinging_momentum_R_max")
6805 layer1_impinging_R_mom_lo = float(_irm_lo)
if _irm_lo
is not None else None
6806 layer1_impinging_R_mom_hi = float(_irm_hi)
if _irm_hi
is not None else None
6813 layer1_R_mom_aim_lo = layer1_R_mom_aim_hi =
None
6814 if layer1_impinging_R_mom_lo
is not None and layer1_impinging_R_mom_hi
is not None:
6815 _r_half = 0.5 * (layer1_impinging_R_mom_hi - layer1_impinging_R_mom_lo)
6816 _r_mid = 0.5 * (layer1_impinging_R_mom_hi + layer1_impinging_R_mom_lo)
6817 layer1_R_mom_aim_lo = _r_mid - 0.8 * _r_half
6818 layer1_R_mom_aim_hi = _r_mid + 0.8 * _r_half
6820 layer1_chamber_dt_ratio_min =
_requirement_float(requirements,
"layer1_chamber_dt_ratio_min", 2.2)
6821 layer1_chamber_dt_ratio_max =
_requirement_float(requirements,
"layer1_chamber_dt_ratio_max", 3.2)
6822 layer1_chamber_ld_ratio_min =
_requirement_float(requirements,
"layer1_chamber_ld_ratio_min", 1.0)
6823 layer1_chamber_ld_ratio_max =
_requirement_float(requirements,
"layer1_chamber_ld_ratio_max", 3.2)
6824 _iang_lo = requirements.get(
"layer1_impinging_angle_deg_min")
6825 _iang_hi = requirements.get(
"layer1_impinging_angle_deg_max")
6826 layer1_impinging_angle_deg_lo = float(_iang_lo)
if _iang_lo
is not None else None
6827 layer1_impinging_angle_deg_hi = float(_iang_hi)
if _iang_hi
is not None else None
6828 layer1_W_impinging_angle =
_requirement_float(requirements,
"W_IMPINGING_ANGLE", 0.0)
6829 layer1_W_impinging_jet_asym =
_requirement_float(requirements,
"W_IMPINGING_JET_ASYM", 160.0)
6831 requirements,
"layer1_impinging_jet_angle_max_asym_deg", 28.0
6834 requirements,
"layer1_exit_pressure_inside_quad_scale", 0.35
6842 requirements,
"layer1_exit_pressure_deadband_rel", 0.05
6847 layer1_Lstar_target_m =
_requirement_float(requirements,
"layer1_Lstar_target_m", 1.0)
6849 requirements,
"layer1_chamber_wall_density_kg_m3", 2000.0
6852 requirements,
"layer1_chamber_mass_ref_kg", 5.0
6862 float(requirements[
"W_DP_O"])
6863 if requirements.get(
"W_DP_O")
is not None
6867 float(requirements[
"W_DP_F"])
6868 if requirements.get(
"W_DP_F")
is not None
6872 layer1_dp_o_band, layer1_dp_f_band = injector_dp_bands_from_requirements(
6874 default_o=_LAYER1_DEFAULT_DP_O_BAND,
6875 default_f=_LAYER1_DEFAULT_DP_F_BAND,
6883 def _interior_band(band, frac: float = 0.8):
6884 lo, hi = float(band[0]), float(band[1])
6885 half = 0.5 * (hi - lo)
6886 mid = 0.5 * (hi + lo)
6887 return (mid - frac * half, mid + frac * half)
6889 layer1_dp_o_aim = _interior_band(layer1_dp_o_band)
6890 layer1_dp_f_aim = _interior_band(layer1_dp_f_band)
6896 _l1_cea_cache = getattr(runner,
"cea_cache",
None)
6898 _l1_cea_cache =
None
6899 _pc_ref_psi = (float(_pc_tgt_seed)
if (_pc_tgt_seed
is not None and float(_pc_tgt_seed) > 0)
6900 else min(max_lox_P_psi, max_fuel_P_psi) * max_P_ratio
6901 / (1.0 + 0.25 * (sum(layer1_dp_o_band) + sum(layer1_dp_f_band))))
6902 _t_b = float(target_burn_time
or 0.0)
6904 requirements, _l1_cea_cache, float(optimal_of), _pc_ref_psi * psi_to_Pa,
6905 float(target_P_exit), (float(exp_ratio_lo), float(exp_ratio_hi)),
6906 float(target_thrust), _t_b, layer1_W_MASS, layer1_chamber_mass_ref_kg)
6907 layer1_W_ISP, layer1_W_MASS = _mw[
"W_ISP"], _mw[
"W_MASS"]
6908 layer1_Isp_ref_s, _eps_ref, layer1_m_prop_ref_kg = _mw[
"Isp_ref_s"], _mw[
"eps_ref"], _mw[
"m_prop_ref_kg"]
6909 if not np.isfinite(layer1_Isp_ref_s)
and _requirement_float(requirements,
"layer1_W_ISP", 1.0) > 0:
6910 layer1_logger.warning(
"Layer 1: no CEA table for the Isp reference; performance term off.")
6912 "Layer 1 figure of merit: W_ISP=%g on Isp_ref %.2f s (ideal, O/F %.3f, Pc %.0f psi, eps %.3f); "
6913 "m_prop_ref %.3f kg for %.0f N x %.2f s; W_MASS=%g (%s)",
6914 layer1_W_ISP, layer1_Isp_ref_s, float(optimal_of), _pc_ref_psi, _eps_ref,
6915 layer1_m_prop_ref_kg, float(target_thrust), _t_b, layer1_W_MASS,
6916 "config" if requirements.get(
"layer1_W_MASS")
is not None else "priced like propellant")
6917 _lio_sf_raw = requirements.get(
"injector_dp_ratio_O_soft_floor")
6918 layer1_dp_o_soft_floor: Optional[float] = (
6920 if _lio_sf_raw
is not None and np.isfinite(float(_lio_sf_raw))
6925 layer1_target_smd_microns =
_requirement_float(requirements,
"target_smd_microns", 50.0)
6928 layer1_tank_equal_scale_psi =
_requirement_float(requirements,
"layer1_tank_equal_scale_psi", 100.0)
6929 layer1_tank_equal_tol_psi =
_requirement_float(requirements,
"layer1_tank_equal_tol_psi", 0.0)
6930 layer1_of_deadband_rel =
_requirement_float(requirements,
"layer1_of_deadband_rel", 0.0)
6936 requirements,
"layer1_enforce_ring_geometry",
True)
6939 requirements,
"layer1_injector_center_clear_dia_m", 0.0)
6941 requirements,
"layer1_injector_min_web_m", 0.0)
6943 requirements,
"layer1_injector_wall_clearance_m", 0.0)
6945 requirements,
"layer1_injector_spray_radius_frac", 0.0)
6947 requirements,
"layer1_injector_spray_radius_tol", 0.08)
6949 requirements,
"layer1_injector_min_face_incidence_deg", 0.0)
6951 requirements,
"layer1_injector_plate_thickness_m", 0.0)
6953 requirements,
"layer1_injector_counterbore_dia_m", 0.0)
6955 requirements,
"layer1_injector_min_back_web_m", 0.0)
6959 layer1_injector_land_ld_O = _pc[
"layer1_injector_land_ld_O"]
6960 layer1_injector_land_ld_F = _pc[
"layer1_injector_land_ld_F"]
6961 layer1_injector_face_contoured = _pc[
"layer1_injector_face_contoured"]
6962 layer1_injector_back_channels = _pc[
"layer1_injector_back_channels"]
6963 layer1_injector_exit_land = _pc[
"layer1_injector_exit_land"]
6964 layer1_injector_channel_width_m = _pc[
"layer1_injector_channel_width_m"]
6965 layer1_injector_channel_floor = _pc[
"layer1_injector_channel_floor"]
6966 layer1_injector_channel_inlets = _pc[
"layer1_injector_channel_inlets"]
6967 layer1_injector_plate_support = _pc[
"layer1_injector_plate_support"]
6968 if layer1_injector_face_contoured
and layer1_injector_back_channels:
6970 layer1_injector_min_face_incidence_deg = 0.0
6972 _env = _assembly_envelope(config_obj)
6973 layer1_injector_r_plate_m = float(_env[
"r_sleeve_id"])
6974 layer1_injector_liner_thickness_m = float(_env[
"liner_thickness"])
6980 layer1_injector_back_keepout_dia_m = float(_ign[
"back_keepout_dia"])
if _ign
else 0.0
6981 layer1_injector_back_keepout_wall_m = float(_ign[
"face_wall"])
if _ign
else 0.0
6982 layer1_injector_r_port_m = 0.5 * float(_ign[
"thread_od"])
if _ign
else 0.0
6984 requirements,
"layer1_resultant_tilt_from_reach",
False)
6986 requirements,
"layer1_resultant_tilt_reach_margin", 1.5)
6988 requirements,
"layer1_chamber_od_increment_in", 0.0
6990 layer1_chamber_od_snap_target = _od_snap_target
6992 requirements,
"layer1_tank_equal_inband_frac", 0.0)
6997 requirements,
"layer1_lock_tank_pressures", layer1_W_TANK_EQUAL > 0.0)
6999 requirements,
"layer1_ring_order_fuel_outboard",
True)
7007 requirements,
"layer1_momentum_wall_side_multiplier", 1.0)
7009 requirements,
"layer1_resultant_tilt_max_deg", 0.0)
7011 requirements,
"layer1_resultant_tilt_scale_deg", 2.0)
7014 requirements,
"layer1_resultant_tilt_outward_multiplier", 25.0)
7016 requirements,
"layer1_resultant_tilt_gate_tol_deg", 1.0)
7018 requirements,
"layer1_momentum_low_side_multiplier", 10.0)
7020 requirements,
"layer1_momentum_band_width", 0.05)
7021 layer1_momentum_scale =
_requirement_float(requirements,
"layer1_momentum_scale", 0.10)
7024 requirements,
"layer1_momentum_gate_safe_slack", 1.0))
7029 requirements,
"layer1_derive_impingement_spacing",
True)
7035 _ld_target, layer1_impingement_Ld_min, layer1_impingement_Ld_max = impingement_ld_band(
7038 layer1_W_geom_ao_af = float(requirements.get(
"W_geom_ao_af_momentum", 0.0))
7046 requirements,
"layer1_lock_tank_pressures",
7050 _w_ol = requirements.get(
"layer1_W_OF_low_MR_scale")
7051 layer1_W_OF_low_MR_scale =
max(1.0, float(_w_ol))
if _w_ol
is not None else 1.0
7052 _w_oh = requirements.get(
"layer1_W_OF_high_MR_scale")
7053 layer1_W_OF_high_MR_scale =
max(1.0, float(_w_oh))
if _w_oh
is not None else 1.0
7056 if l1_injector_type ==
"impinging" and layer1_W_geom_ao_af <= 0.0:
7057 layer1_W_geom_ao_af = 0.0
if _tanks_locked_for_geom
else max(400.0, 0.05 * layer1_W_OF_obj)
7058 if l1_injector_type ==
"impinging" and _tanks_locked_for_geom
and layer1_W_geom_ao_af > 0.0:
7063 "Layer 1: W_geom_ao_af_momentum=%g ignored -- tank pressures are locked, so the "
7064 "momentum ratio is fixed by construction and O/F alone owns the injector area ratio.",
7065 layer1_W_geom_ao_af)
7066 layer1_W_geom_ao_af = 0.0
7067 layer1_of_validation_tol =
_requirement_float(requirements,
"layer1_of_validation_tol", 0.15)
7069 _thr_gate_yaml = requirements.get(
"layer1_thrust_validation_rel_tol")
7070 layer1_gate_thrust_tol = (
7071 float(_thr_gate_yaml)
7072 if _thr_gate_yaml
is not None and np.isfinite(float(_thr_gate_yaml))
and float(_thr_gate_yaml) > 0
7073 else float(thrust_tol)
7076 "Injector ΔP_inj/Pc preferred bands (soft hinge): "
7077 f
"oxidizer [{layer1_dp_o_band[0]:.3f}, {layer1_dp_o_band[1]:.3f}], "
7078 f
"fuel [{layer1_dp_f_band[0]:.3f}, {layer1_dp_f_band[1]:.3f}]"
7081 "Injector ΔP weights: W_DP=%g (fallback); W_DP_O=%g; W_DP_F=%g",
7086 if layer1_dp_o_soft_floor
is not None and layer1_dp_o_soft_floor > 0
and layer1_W_DP_O_FLOOR > 0:
7088 "LOX dP_inj/Pc soft floor: penalize ratios below %.4f using W_DP_O_FLOOR=%.4g*(floor-r_O)^2",
7089 layer1_dp_o_soft_floor,
7090 layer1_W_DP_O_FLOOR,
7092 if layer1_W_geom_ao_af > 0.0:
7094 f
"Geometry AO/AF momentum hint weight W_geom_ao_af_momentum={layer1_W_geom_ao_af:g} "
7095 f
"(targets A_O/A_F ≈ MR/√(ρ_O/ρ_F) for R≈1)"
7097 if layer1_impinging_R_mom_lo
is not None and layer1_impinging_R_mom_hi
is not None:
7099 "Impinging momentum_ratio_R: validation band "
7100 f
"[{layer1_impinging_R_mom_lo:g}, {layer1_impinging_R_mom_hi:g}], optimizer AIM band "
7101 f
"[{layer1_R_mom_aim_lo:g}, {layer1_R_mom_aim_hi:g}] "
7102 f
"(penalty pushes inside the aim band so validation passes with margin)"
7104 if layer1_W_impinging_angle > 0.0
and layer1_impinging_angle_deg_lo
is not None and layer1_impinging_angle_deg_hi
is not None:
7106 "Impinging effective-angle preferred band (deg): "
7107 f
"[{layer1_impinging_angle_deg_lo:g}, {layer1_impinging_angle_deg_hi:g}] "
7108 f
"with W_IMPINGING_ANGLE={layer1_W_impinging_angle:g}"
7110 if l1_injector_type ==
"impinging" and layer1_W_impinging_jet_asym > 0.0:
7112 "Impinging paired jet |Δθ_O/F| soft cap: max Δ=%.2f° with W_IMPINGING_JET_ASYM=%g",
7113 float(layer1_impinging_jet_max_asym_deg),
7114 float(layer1_W_impinging_jet_asym),
7117 "Ambient-matched exit pressure: inside-deadband quadratic scale "
7118 "layer1_exit_pressure_inside_quad_scale=%g (W_EXIT×scale×rel²; 0=flat inside band)",
7119 float(layer1_exit_pressure_inside_quad),
7122 "Layer 1 primary weights: layer1_W_THRUST=%g layer1_W_OF=%g low_MR_scale=%g high_MR_scale=%g; "
7123 "of_validation_tol=%.4f",
7124 layer1_W_THRUST_obj,
7126 layer1_W_OF_low_MR_scale,
7127 layer1_W_OF_high_MR_scale,
7128 layer1_of_validation_tol,
7132 "marginal_candidate_logged":
False,
7135 update_progress(
"Layer 1: Objective", 0.10,
"Defining objective function...")
7149 _ocr = requirements.get(
"objective_cache_rel")
7150 cache_rel = float(_ocr)
if _ocr
is not None else 1e-5
7151 cache_rel = float(np.clip(cache_rel, 1e-6, 1e-3))
7152 cache_steps = np.maximum(span * cache_rel, 1e-12)
7153 eval_cache: Dict[Tuple[int, ...], Dict[str, Any]] = {}
7155 def _make_eval_cache_key(x_raw: np.ndarray) -> Tuple[int, ...]:
7156 """Quantize x to a stable hashable key for caching evaluate() results."""
7157 x_arr = np.clip(np.asarray(x_raw, dtype=float), lower_bounds, upper_bounds)
7158 x_arr = x_arr.copy()
7159 for di
in layer1_integer_dims:
7160 if 0 <= di < len(x_arr):
7161 x_arr[di] = int(round(x_arr[di]))
7163 for i, v
in enumerate(x_arr):
7164 if i
in layer1_integer_dims:
7165 key_parts.append(int(v))
7167 step = float(cache_steps[i])
7169 key_parts.append(int(round((float(v) - float(lower_bounds[i])) / step)))
7170 return tuple(key_parts)
7172 def _check_valley_escape_tier() -> int:
7174 opt_state[
"function_evaluations"], opt_state[
"best_objective"],
7175 opt_state[
"function_evaluations"] - opt_state[
"last_best_eval"])
7178 def objective(x: np.ndarray) -> float:
7179 """Layer 1 objective function: optimize geometry + initial pressures.
7181 Staged/lexicographic structure (approximate but smooth):
7182 1) Feasibility (hard constraints + stability gates)
7185 4) Regularization (Cf band, chamber length, exit pressure)
7187 Note: we still return a single scalar for SciPy/CMA-ES, but the scaling
7188 preserves the priority ordering and reduces "weight fights".
7192 for key
in [
"consecutive_failures",
"last_valid_obj",
"iteration",
"function_evaluations",
"best_objective",
"best_x"]:
7193 if key
not in opt_state:
7194 opt_state[key] = 0
if key
in [
"iteration",
"function_evaluations",
"consecutive_failures"]
else (float(
'inf')
if key
in [
"last_valid_obj",
"best_objective"]
else None)
7196 opt_state[
"iteration"] += 1
7197 iteration = opt_state[
"iteration"]
7198 opt_state[
"function_evaluations"] += 1
7206 eval_budget = int(opt_state.get(
"eval_budget", 0))
7208 progress = 0.10 + 0.40 *
min(iteration / eval_budget, 1.0)
7209 eval_str = f
"{iteration}/{eval_budget}"
7212 eval_str = f
"{iteration}"
7213 if iteration <= 3
or iteration % 25 == 0:
7215 _bo = float(opt_state[
"best_objective"])
7216 except (TypeError, ValueError):
7218 _lv_raw = opt_state.get(
"last_valid_obj", float(
"inf"))
7220 _lv = float(_lv_raw)
if _lv_raw
is not None else float(
"inf")
7221 except (TypeError, ValueError):
7223 best_obj_str = f
"{_bo:.3e}" if np.isfinite(_bo)
else "inf"
7224 curr_obj_str = f
"{_lv:.3e}" if np.isfinite(_lv)
else "inf"
7225 update_progress(
"Layer 1: Optimization", progress, f
"Eval {eval_str} | Curr: {curr_obj_str} | Best: {best_obj_str}")
7226 layer1_logger.info(f
"[{int(progress*100)}%] Evaluation {eval_str} - "
7227 f
"Objective: {curr_obj_str} (Best: {best_obj_str})")
7228 for handler
in layer1_logger.handlers:
7232 x_clipped = np.clip(np.asarray(x, dtype=float), lower_bounds_eff, upper_bounds_eff)
7241 x_clipped = x_clipped.copy()
7242 for di
in layer1_integer_dims:
7243 if 0 <= di < len(x_clipped):
7244 x_clipped[di] = int(round(x_clipped[di]))
7247 config, _, _ = apply_x_to_config(x_clipped, config_base)
7251 cg = ensure_chamber_geometry(config)
7252 A_throat_check = float(cg.A_throat
or 0.0)
7253 D_chamber_inner = float(cg.chamber_diameter
or 0.0)
7254 A_chamber_check = np.pi * (D_chamber_inner / 2) ** 2
7256 geom = getattr(getattr(config,
"injector",
None),
"geometry",
None)
7257 has_pintle = bool(getattr(getattr(config,
"injector",
None),
"type",
None) ==
"pintle" and geom
is not None)
7258 has_impinging = bool(getattr(getattr(config,
"injector",
None),
"type",
None) ==
"impinging" and geom
is not None)
7260 A_lox_injector = np.nan
7261 A_fuel_injector = np.nan
7264 area_ratio_error = np.nan
7266 infeasibility_score = 0.0
7272 if A_chamber_check > 0
and A_throat_check > 0:
7273 contraction_ratio_check = A_chamber_check / A_throat_check
7275 infeasibility_score +=
max(0.0, (A_throat_check * 1.1) / A_chamber_check - 1.0) ** 2
7278 if contraction_ratio_check < 1.5:
7279 infeasibility_score += (1.5 - contraction_ratio_check) ** 2
7280 elif contraction_ratio_check > 15.0:
7281 infeasibility_score += (contraction_ratio_check - 15.0) ** 2
7284 if has_pintle
and D_chamber_inner > 0:
7285 d_pintle_tip_check = float(geom.fuel.d_pintle_tip)
7287 infeasibility_score +=
max(0.0, (d_pintle_tip_check * 1.1) / D_chamber_inner - 1.0) ** 2
7291 D_throat_check = np.sqrt(4.0 * A_throat_check / np.pi)
if A_throat_check > 0
else 0.0
7292 if D_throat_check > 0
and D_chamber_inner > 0:
7295 infeasibility_score +=
max(0.0, D_throat_check - D_chamber_inner * 0.95) ** 2 * 10.0
7298 D_exit_check = float(cg.exit_diameter
or 0.0)
7299 R_throat_check = np.sqrt(
max(0, A_throat_check / np.pi))
7300 D_throat_check = 2 * R_throat_check
7301 if D_exit_check > 0
and D_throat_check > 0:
7303 infeasibility_score +=
max(0.0, D_throat_check / D_exit_check - 1.0) ** 2
7307 if has_pintle
and A_throat_check > 0:
7308 A_lox_injector = float(geom.lox.n_orifices * np.pi * (geom.lox.d_orifice / 2) ** 2)
7309 R_inner = float(geom.fuel.d_pintle_tip / 2)
7310 R_outer = float(R_inner + geom.fuel.h_gap)
7311 A_fuel_injector = float(np.pi * (R_outer ** 2 - R_inner ** 2))
7314 if has_impinging
and A_throat_check > 0:
7317 A_lox_injector = float(oxg.n_elements * np.pi * (oxg.d_jet / 2.0) ** 2)
7318 A_fuel_injector = float(fug.n_elements * np.pi * (fug.d_jet / 2.0) ** 2)
7319 sp_O = float(oxg.spacing)
7320 sp_F = float(fug.spacing)
7321 djo = float(oxg.d_jet)
7322 djf = float(fug.d_jet)
7328 n_el_O=int(oxg.n_elements),
7329 n_el_F=int(fug.n_elements),
7330 D_chamber_inner=D_chamber_inner,
7334 P_O_psi = float(np.clip(x_clipped[idx_P_O], bounds[idx_P_O][0], bounds[idx_P_O][1]))
7335 P_F_psi = float(np.clip(x_clipped[idx_P_F], bounds[idx_P_F][0], bounds[idx_P_F][1]))
7336 P_O_test = P_O_psi * psi_to_Pa
7337 P_F_test = P_F_psi * psi_to_Pa
7338 P_O_ratio = P_O_psi / max_lox_P_psi
if max_lox_P_psi > 0
else 0.0
7339 P_F_ratio = P_F_psi / max_fuel_P_psi
if max_fuel_P_psi > 0
else 0.0
7341 skip_physics_eval =
False
7342 if has_impinging
and A_throat_check > 0
and geom
is not None:
7346 d_jet_O=float(oxg.d_jet),
7347 d_jet_F=float(fug.d_jet),
7348 sp_O=float(oxg.spacing),
7349 sp_F=float(fug.spacing),
7350 D_chamber_inner=D_chamber_inner,
7351 D_throat_check=D_throat_check,
7352 A_chamber_check=A_chamber_check,
7353 A_throat_check=A_throat_check,
7354 n_elements=float(oxg.n_elements),
7355 angle_O_deg=float(oxg.impingement_angle),
7356 angle_F_deg=float(fug.impingement_angle),
7357 center_clear_dia_m=layer1_injector_center_clear_dia_m,
7358 min_web_m=layer1_injector_min_web_m,
7359 wall_clearance_m=layer1_injector_wall_clearance_m,
7360 min_face_incidence_deg=layer1_injector_min_face_incidence_deg,
7361 spray_radius_frac=layer1_injector_spray_radius_frac,
7362 spray_radius_tol=layer1_injector_spray_radius_tol,
7363 plate_thickness_m=layer1_injector_plate_thickness_m,
7364 counterbore_dia_m=layer1_injector_counterbore_dia_m,
7365 land_ld_O=layer1_injector_land_ld_O,
7366 land_ld_F=layer1_injector_land_ld_F,
7367 min_back_web_m=layer1_injector_min_back_web_m,
7368 back_keepout_dia_m=layer1_injector_back_keepout_dia_m,
7369 back_keepout_wall_m=layer1_injector_back_keepout_wall_m,
7370 face_contoured=layer1_injector_face_contoured,
7371 exit_land=layer1_injector_exit_land,
7372 back_channels=layer1_injector_back_channels,
7373 passage_ld_O=layer1_injector_land_ld_O,
7374 passage_ld_F=layer1_injector_land_ld_F,
7375 channel_width_m=layer1_injector_channel_width_m,
7376 channel_floor=layer1_injector_channel_floor,
7377 r_plate_m=layer1_injector_r_plate_m,
7378 r_port_m=layer1_injector_r_port_m,
7379 liner_thickness_m=layer1_injector_liner_thickness_m,
7381 elif infeasibility_score > 0.0:
7382 skip_physics_eval =
True
7384 eval_success =
False
7385 final_results: Dict[str, Any] = {}
7386 final_pressures = (P_O_test, P_F_test)
7387 eval_error_str: Optional[str] =
None
7393 cache_enabled =
not bool(opt_state.get(
"layer1_lbfgs_phase",
False))
7394 cache_key = _make_eval_cache_key(x_clipped)
if cache_enabled
else None
7395 if not skip_physics_eval:
7396 cached = eval_cache.get(cache_key)
if cache_enabled
and cache_key
is not None else None
7398 if cached
is not None:
7400 if isinstance(cached, dict):
7402 if "results" in cached:
7403 eval_success = bool(cached.get(
"success",
False))
7404 final_results = copy.deepcopy(cached.get(
"results", {}))
if eval_success
else {}
7405 eval_error_str = cached.get(
"error",
None)
7411 _xs = cached.get(
"x_solved")
7412 if _xs
is not None and len(_xs) == len(x_clipped):
7413 x_clipped[:] = np.asarray(_xs, dtype=float)
7416 A_throat=float(x_clipped[0]),
7417 Lstar=float(x_clipped[1]),
7418 expansion_ratio=float(x_clipped[2]),
7420 float(x_clipped[3]), _od_increment_in,
7421 wall_m=wall_total_m, snap_target=_od_snap_target),
7422 max_nozzle_exit=max_nozzle_exit,
7423 wall_thickness_m=wall_total_m,
7437 config_runner = copy.deepcopy(config)
7438 if os.environ.get(
"ED_L1_INLINE_THERMAL_OFF",
"0") ==
"1":
7439 if hasattr(config_runner,
"ablative_cooling")
and config_runner.ablative_cooling:
7440 config_runner.ablative_cooling.enabled =
False
7441 if hasattr(config_runner,
"graphite_insert")
and config_runner.graphite_insert:
7442 config_runner.graphite_insert.enabled =
False
7444 test_runner = PintleEngineRunner(config_runner)
7448 def _evaluate_once_slow():
7453 from engine
import accel
as _accel
7454 _r = _accel.evaluate(test_runner.config, test_runner.cea_cache,
7455 P_O_test, P_F_test, target_P_exit)
7458 return test_runner.evaluate(
7459 P_O_test, P_F_test, P_ambient=target_P_exit, silent=
True)
7460 except Exception
as _e:
7461 _last_err[0] = str(_e)
7469 for _grow
in (1.0, 1.3, 1.7, 2.2):
7471 x_clipped[0] = float(
min(x_clipped[0] * (_grow / _prev_grow),
7472 derive_At_range[1]))
7474 config_runner, A_throat=float(x_clipped[0]),
7475 Lstar=float(x_clipped[1]), expansion_ratio=float(x_clipped[2]),
7477 float(x_clipped[3]), _od_increment_in,
7478 wall_m=wall_total_m, snap_target=_od_snap_target),
7479 max_nozzle_exit=max_nozzle_exit, wall_thickness_m=wall_total_m)
7481 _first = _evaluate_once_slow()
7482 _err = (_last_err[0]
or "")
7483 if (_first
is not None or not derive_throat_from_thrust
7484 or not (
"Supply < Demand" in _err
or "Insufficient mass flow" in _err)):
7487 eval_success =
False
7488 eval_error_str = _last_err[0]
or "evaluate() raised"
7493 x_clipped, config_runner, _derive_constants,
7494 _evaluate_once_slow, _first,
7495 cea_cache=getattr(test_runner,
"cea_cache",
None),
7497 final_results = _solved
if isinstance(_solved, dict)
else _first
7503 A_throat=float(x_clipped[0]),
7504 Lstar=float(x_clipped[1]),
7505 expansion_ratio=float(x_clipped[2]),
7507 float(x_clipped[3]), _od_increment_in,
7508 wall_m=wall_total_m, snap_target=_od_snap_target),
7509 max_nozzle_exit=max_nozzle_exit,
7510 wall_thickness_m=wall_total_m,
7513 if cache_enabled
and cache_key
is not None:
7514 eval_cache[cache_key] = {
7515 "success": bool(eval_success),
7516 "results": copy.deepcopy(final_results)
if eval_success
else {},
7517 "error": eval_error_str,
7518 "x_solved": x_clipped.tolist(),
7522 if eval_success
and isinstance(final_results, dict):
7523 final_results[
"eps"] = float(final_results.get(
"eps", x_clipped[2]))
7525 final_results
if eval_success
else None, x_clipped, requirements, constants_dict,
7526 cea_cache=_l1_cea_cache, eval_error=eval_error_str)
7527 obj = float(_T[
"objective"])
7529 infeasibility_score = float(_T.get(
"infeasibility_score", 1.0))
7530 inf_residual = float(_T.get(
"inf_residual", infeasibility_score))
7531 length_violation = bool(_T.get(
"length_violation",
False))
7532 raw_feasible_obj = obj
if (eval_success
and inf_residual <= 0.0
and not length_violation)
else float(
"nan")
7534 F_actual = float(_T.get(
"F", np.nan))
7535 Isp_actual = float(_T.get(
"Isp", np.nan))
7536 MR_actual = float(_T.get(
"MR", np.nan))
7537 Pc_actual = float(_T.get(
"Pc", np.nan))
7538 Cf_actual = float(_T.get(
"Cf", np.nan))
7539 thrust_error = float(_T.get(
"thrust_error", 1.0))
7540 of_error = float(_T.get(
"of_error", 1.0))
7541 stability = _T.get(
"stability_results", {})
or {}
7542 stability_state = _T.get(
"stability_state",
"unstable")
7543 stability_score = float(_T.get(
"stability_score", 0.0))
7544 chugging_margin = float(_T.get(
"chugging_margin", 0.0))
7545 acoustic_margin = float(_T.get(
"acoustic_margin", 0.0))
7546 feed_margin = float(_T.get(
"feed_margin", 0.0))
7547 state_ok = bool(_T.get(
"state_ok",
False))
7548 effective_min_score = float(_T.get(
"effective_min_score", 0.0))
7549 effective_margin = float(_T.get(
"effective_margin", 0.0))
7550 ratio_o_obj = _T.get(
"injector_dp_ratio_O")
7551 ratio_f_obj = _T.get(
"injector_dp_ratio_F")
7552 L_chamber_curr = float(_T.get(
"L_chamber", np.nan))
7553 smd_eff_um = float(_T.get(
"effective_smd_microns", np.nan))
7554 imp_angle_deg = float(_T.get(
"impingement_angle_deg", np.nan))
7555 cf_hinge = float(_T.get(
"cf_hinge", 0.0))
7558 if eval_success
and isinstance(final_results, dict)
and np.isfinite(F_actual):
7559 opt_state[
"last_good_eval_bundle"] = {
7560 "results": copy.deepcopy(final_results),
7561 "P_O_Pa": float(P_O_test),
7562 "P_F_Pa": float(P_F_test),
7563 "thrust_error": float(thrust_error),
7564 "of_error": float(of_error),
7568 thrust_tol_validation = float(layer1_gate_thrust_tol)
7569 of_tol_validation = float(layer1_of_validation_tol)
7571 _ok_o_dp = injector_dp_ratio_within_gate(ratio_o_obj, layer1_dp_o_band[0], layer1_dp_o_band[1])
7572 if _ok_o_dp
is not None:
7573 dp_gate_obj &= bool(_ok_o_dp)
7574 _ok_f_dp = injector_dp_ratio_within_gate(ratio_f_obj, layer1_dp_f_band[0], layer1_dp_f_band[1])
7575 if _ok_f_dp
is not None:
7576 dp_gate_obj &= bool(_ok_f_dp)
7577 momentum_gate_obj =
True
7578 _R_m_gate = float(_T.get(
"momentum_ratio_R", np.nan))
7579 if (has_impinging
and np.isfinite(_R_m_gate)
7580 and layer1_impinging_R_mom_lo
is not None and layer1_impinging_R_mom_hi
is not None):
7581 momentum_gate_obj = (
7582 float(layer1_impinging_R_mom_lo) <= _R_m_gate <= float(layer1_impinging_R_mom_hi))
7583 errors_acceptable = (
7586 and (thrust_error <= thrust_tol_validation + 1e-12)
7587 and (of_error <= of_tol_validation + 1e-12)
7588 and (stability_score >= effective_min_score * 0.8)
7590 and momentum_gate_obj
7592 if errors_acceptable:
7593 if final_pressures
is not None:
7594 opt_state[
"best_pressures"] = final_pressures
7595 opt_state[
"best_results_for_validation"] = {
7596 "F": F_actual,
"MR": MR_actual,
7597 "thrust_error": thrust_error,
"of_error": of_error,
7598 "stability_score": stability_score,
"stability_state": stability_state,
7599 "chugging_margin": chugging_margin,
"acoustic_margin": acoustic_margin,
7600 "feed_margin": feed_margin,
"stability_results": stability,
7602 if not opt_state.get(
'acceptable_found_logged',
False):
7603 log_status(
"Layer 1", f
"✓ Acceptable solution found! Obj={obj:.6e}, Thrust err: {thrust_error*100:.2f}%, O/F err: {of_error*100:.2f}% (continuing optimization...)")
7604 opt_state[
'acceptable_found_logged'] =
True
7605 if np.isfinite(obj)
and obj <= obj_tolerance:
7606 opt_state[
"objective_satisfied"] =
True
7607 if obj < float(opt_state.get(
"satisfied_obj", float(
"inf"))):
7608 opt_state[
"satisfied_obj"] = float(obj)
7609 opt_state[
"satisfied_eval_count"] = int(opt_state.get(
"satisfied_eval_count", 0)) + 1
7612 if eval_success
and np.isfinite(obj):
7613 opt_state[
'consecutive_failures'] = 0
7614 opt_state[
'last_valid_obj'] = obj
7616 opt_state[
'consecutive_failures'] += 1
7619 def _finite_or_none(v: Any) -> Optional[float]:
7624 return vv
if np.isfinite(vv)
else None
7629 _at_clip = derive_At_range
if derive_throat_from_thrust
else bounds[0]
7630 _eps_clip = derive_eps_range
if derive_expansion_ratio
else bounds[2]
7631 A_throat_curr = float(np.clip(x_clipped[0], _at_clip[0], _at_clip[1]))
7632 Lstar_curr = float(np.clip(x_clipped[1], bounds[1][0], bounds[1][1]))
7633 expansion_ratio_curr = float(np.clip(x_clipped[2], _eps_clip[0], _eps_clip[1]))
7638 float(np.clip(x_clipped[3], bounds[3][0], bounds[3][1])),
7639 _od_increment_in, wall_m=wall_total_m, snap_target=_od_snap_target,
7641 if l1_injector_type ==
"impinging":
7642 n_doublets_curr = int(round(np.clip(x_clipped[4], bounds[4][0], bounds[4][1])))
7643 n_elements_O_curr = n_doublets_curr
7644 d_jet_O_curr = float(np.clip(x_clipped[5], bounds[5][0], bounds[5][1]))
7645 imp_ang_O_curr = float(np.clip(x_clipped[6], bounds[6][0], bounds[6][1]))
7646 spacing_O_curr = float(np.clip(x_clipped[7], bounds[7][0], bounds[7][1]))
7647 n_elements_F_curr = n_doublets_curr
7648 d_jet_F_curr = float(np.clip(x_clipped[8], bounds[8][0], bounds[8][1]))
7649 imp_ang_F_curr = float(np.clip(x_clipped[9], bounds[9][0], bounds[9][1]))
7650 spacing_F_curr = float(np.clip(x_clipped[10], bounds[10][0], bounds[10][1]))
7651 d_pintle_tip_curr = float(
"nan")
7652 h_gap_curr = float(
"nan")
7653 n_orifices_curr = -1
7654 d_orifice_curr = float(
"nan")
7656 d_pintle_tip_curr = float(np.clip(x_clipped[4], bounds[4][0], bounds[4][1]))
7657 h_gap_curr = float(np.clip(x_clipped[5], bounds[5][0], bounds[5][1]))
7658 n_orifices_curr = int(round(np.clip(x_clipped[6], bounds[6][0], bounds[6][1])))
7659 d_orifice_curr = float(np.clip(x_clipped[7], bounds[7][0], bounds[7][1]))
7660 n_elements_O_curr = -1
7661 n_elements_F_curr = -1
7662 d_jet_O_curr = float(
"nan")
7663 d_jet_F_curr = float(
"nan")
7664 imp_ang_O_curr = float(
"nan")
7665 imp_ang_F_curr = float(
"nan")
7666 spacing_O_curr = float(
"nan")
7667 spacing_F_curr = float(
"nan")
7668 P_O_start_psi_hist = float(np.clip(x_clipped[idx_P_O], bounds[idx_P_O][0], bounds[idx_P_O][1]))
7669 P_F_start_psi_hist = float(np.clip(x_clipped[idx_P_F], bounds[idx_P_F][0], bounds[idx_P_F][1]))
7671 D_inner_curr = D_outer_curr - wall_total_m
7672 if D_inner_curr <= 0:
7673 D_inner_curr =
max(D_outer_curr * 0.3, 0.01)
7675 lox_start_ratio_hist = P_O_start_psi_hist / max_lox_P_psi
if max_lox_P_psi > 0
else 0.7
7676 fuel_start_ratio_hist = P_F_start_psi_hist / max_fuel_P_psi
if max_fuel_P_psi > 0
else 0.7
7677 combined_stability_margin =
min(chugging_margin, acoustic_margin, feed_margin)
if eval_success
else 0.0
7678 L_chamber_hist = L_chamber_curr
if (np.isfinite(L_chamber_curr))
else None
7679 _R_hist = float(
"nan")
7680 if has_impinging
and eval_success:
7681 _Rv_hist = final_results.get(
"diagnostics", {}).get(
"momentum_ratio_R")
7682 if _Rv_hist
is not None:
7684 _R_hist = float(_Rv_hist)
7685 except (TypeError, ValueError):
7686 _R_hist = float(
"nan")
7687 _lstar_lo_hist = float(requirements.get(
"min_Lstar", 0.76))
7688 _lstar_hi_hist = float(requirements.get(
"max_Lstar", 1.5))
7689 _exp_lo_hist = float(exp_ratio_lo)
7690 _exp_hi_hist = float(exp_ratio_hi)
7692 _candidate_new_best = (
7695 and float(obj) < float(opt_state.get(
"best_objective", float(
"inf")))
7698 bool(opt_state.get(
"layer1_lbfgs_phase"))
7699 and not _candidate_new_best
7703 or float(inf_residual) > 0.0
7706 if not _skip_history:
7707 opt_state[
"history"].append({
7708 "iteration": iteration,
7709 "x": x_clipped.copy(),
7711 "A_throat": A_throat_curr,
7712 "Lstar": Lstar_curr,
7713 "expansion_ratio": expansion_ratio_curr,
7714 "D_chamber_outer": D_outer_curr,
7715 "D_chamber_inner": D_inner_curr,
7716 "L_chamber": L_chamber_hist,
7717 "d_pintle_tip": d_pintle_tip_curr,
7718 "h_gap": h_gap_curr,
7719 "n_orifices": n_orifices_curr,
7720 "d_orifice": d_orifice_curr,
7721 "n_elements_O": n_elements_O_curr,
7722 "n_elements_F": n_elements_F_curr,
7723 "n_doublets": float(n_elements_O_curr)
if l1_injector_type ==
"impinging" else float(
"nan"),
7724 "d_jet_O": d_jet_O_curr,
7725 "impingement_angle_O": imp_ang_O_curr,
7726 "spacing_O": spacing_O_curr,
7727 "d_jet_F": d_jet_F_curr,
7728 "impingement_angle_F": imp_ang_F_curr,
7729 "spacing_F": spacing_F_curr,
7730 "effective_smd_microns": _finite_or_none(smd_eff_um),
7731 "impingement_angle_deg_effective": _finite_or_none(imp_angle_deg),
7732 "P_O_start_psi": P_O_start_psi_hist,
7733 "P_F_start_psi": P_F_start_psi_hist,
7735 "thrust": _finite_or_none(F_actual),
7736 "thrust_error": thrust_error,
7737 "of_error": of_error,
7738 "Isp": _finite_or_none(Isp_actual),
7739 "MR": _finite_or_none(MR_actual),
7740 "Pc": _finite_or_none(Pc_actual),
7741 "Cf": _finite_or_none(Cf_actual),
7742 "Cf_error": float(np.sqrt(cf_hinge))
if np.isfinite(cf_hinge)
else 1.0,
7743 "Cf_penalty": float(_tw.get(
"cf_penalty", 0.0)),
7745 "stability_margin": combined_stability_margin,
7746 "stability_state": stability_state,
7747 "stability_score": stability_score,
7748 "chugging_margin": chugging_margin,
7749 "acoustic_margin": acoustic_margin,
7750 "feed_margin": feed_margin,
7752 "lox_end_ratio": lox_start_ratio_hist,
7753 "fuel_end_ratio": fuel_start_ratio_hist,
7754 "lox_start_ratio": lox_start_ratio_hist,
7755 "fuel_start_ratio": fuel_start_ratio_hist,
7757 "infeasibility_score": float(infeasibility_score),
7758 "eval_success": bool(eval_success),
7759 "eval_error": eval_error_str,
7761 "injector_dp_ratio_O": _finite_or_none(ratio_o_obj),
7762 "injector_dp_ratio_F": _finite_or_none(ratio_f_obj),
7763 "momentum_ratio_R": _finite_or_none(_R_hist),
7764 "raw_feasible_objective": _finite_or_none(raw_feasible_obj),
7765 "penalty_thrust": float(_tw.get(
"thrust_penalty", 0.0)),
7766 "penalty_pc_target": float(_tw.get(
"pc_target_penalty", 0.0)),
7767 "penalty_of": float(_tw.get(
"of_penalty", 0.0)),
7768 "penalty_cf": float(_tw.get(
"cf_penalty", 0.0)),
7769 "penalty_exit": float(_tw.get(
"exit_pressure_penalty", 0.0)),
7770 "penalty_injector_dp": float(_tw.get(
"injector_dp_penalty", 0.0)),
7771 "penalty_length": float(_tw.get(
"length_penalty", 0.0)),
7772 "penalty_chamber_shape": float(_tw.get(
"chamber_shape_penalty", 0.0)),
7773 "penalty_chamber_mass": float(_tw.get(
"chamber_mass_penalty", 0.0)),
7774 "penalty_isp": float(_tw.get(
"isp_penalty", 0.0)),
7775 "Isp_merit": _finite_or_none(_T.get(
"Isp_merit_s")),
7776 "penalty_lstar": float(_tw.get(
"lstar_penalty", 0.0)),
7777 "chamber_mass_kg": float(_T.get(
"chamber_mass_kg", float(
"nan"))),
7778 "penalty_momentum": float(_tw.get(
"momentum_balance_penalty", 0.0)),
7779 "penalty_impingement_angle": float(_tw.get(
"impingement_angle_penalty", 0.0)),
7780 "penalty_jet_asymmetry": float(_tw.get(
"jet_angle_asymmetry_penalty", 0.0)),
7781 "penalty_smd": float(_tw.get(
"smd_penalty", 0.0)),
7782 "penalty_tank_equal": float(_tw.get(
"tank_pressure_equal_penalty", 0.0)),
7783 "penalty_geom_fit": float(_tw.get(
"impinging_geometry_fit_penalty", 0.0)),
7784 "penalty_geom_ao_af": float(_tw.get(
"geom_ao_af_momentum_penalty", 0.0)),
7785 "band_violation_Lstar_sq": float(
7787 )
if np.isfinite(Lstar_curr)
else None,
7788 "band_violation_expansion_sq": float(
7791 "band_violation_dp_O_sq": float(
7794 r_band_lo=float(layer1_dp_o_aim[0]),
7795 r_band_hi=float(layer1_dp_o_aim[1]),
7797 )
if ratio_o_obj
is not None and np.isfinite(float(ratio_o_obj))
else None,
7798 "band_violation_dp_F_sq": float(
7801 r_band_lo=float(layer1_dp_f_aim[0]),
7802 r_band_hi=float(layer1_dp_f_aim[1]),
7804 )
if ratio_f_obj
is not None and np.isfinite(float(ratio_f_obj))
else None,
7805 "band_violation_momentum_R_sq": float(
7808 r_band_lo=layer1_R_mom_aim_lo,
7809 r_band_hi=layer1_R_mom_aim_hi,
7813 and np.isfinite(_R_hist)
7814 and layer1_impinging_R_mom_lo
is not None
7815 and layer1_impinging_R_mom_hi
is not None
7817 "band_violation_angle_sq": float(
7820 angle_band_lo_deg=layer1_impinging_angle_deg_lo,
7821 angle_band_hi_deg=layer1_impinging_angle_deg_hi,
7823 )
if layer1_W_impinging_angle > 0.0
else None,
7832 and float(obj) < float(opt_state.get(
"best_objective", float(
"inf")))
7835 opt_state[
"best_objective"] = obj
7836 opt_state[
"best_x"] = x_clipped.copy()
7837 opt_state[
"last_best_eval"] = opt_state[
"function_evaluations"]
7839 opt_state[
"last_eval_config"] = copy.deepcopy(config)
7840 opt_state[
"last_eval_config_x"] = x_clipped.copy()
7841 opt_state[
"last_eval_validation_tank_pa"] = (float(P_O_test), float(P_F_test))
7852 _breakdown_terms = dict(_tw)
7853 _penalty_sum = float(sum(_breakdown_terms.values()))
7854 opt_state[
"best_objective_breakdown"] = {
7855 "objective": float(obj),
7857 "penalty_sum": _penalty_sum,
7860 "penalty_unaccounted": float(obj) - _penalty_sum,
7862 "injector_dp_ratio_O": ratio_o_obj,
7863 "injector_dp_ratio_F": ratio_f_obj,
7864 "chamber_D_over_Dt": _finite_or_none(_T.get(
"chamber_D_over_Dt")),
7865 "chamber_L_over_D": _finite_or_none(_T.get(
"chamber_L_over_D")),
7866 "Lstar_m": _finite_or_none(Lstar_curr),
7867 "Lstar_target_m": _finite_or_none(_T.get(
"Lstar_target_m")),
7868 "chamber_mass_kg": _finite_or_none(_T.get(
"chamber_mass_kg")),
7869 "Isp_s": _finite_or_none(Isp_actual),
7870 "Isp_merit_s": _finite_or_none(_T.get(
"Isp_merit_s")),
7871 "Isp_ref_s": _finite_or_none(_T.get(
"Isp_ref_s")),
7872 "propellant_mass_kg": _finite_or_none(_T.get(
"propellant_mass_kg")),
7873 "L_face_to_throat_m": _finite_or_none(_T.get(
"L_face_to_throat")),
7874 "L_engine_m": _finite_or_none(_T.get(
"L_engine")),
7875 "momentum_ratio_R": _finite_or_none(_R_hist),
7876 "resultant_tilt_deg": _finite_or_none(_T.get(
"resultant_tilt_deg")),
7877 "resultant_tilt_breakeven_deg": _finite_or_none(_T.get(
"resultant_tilt_breakeven_deg")),
7878 "effective_smd_microns": _finite_or_none(smd_eff_um),
7879 "impingement_angle_deg": _finite_or_none(imp_angle_deg),
7880 "geom_ao_af_momentum_scale": _finite_or_none(_T.get(
"geom_ao_af_momentum_scale")),
7881 "geom_ao_af": _finite_or_none(_T.get(
"geom_ao_af")),
7882 "expected_ao_af_for_R1": _finite_or_none(_T.get(
"expected_ao_af_for_R1")),
7883 "infeasibility_penalty": float(obj - _penalty_sum)
if (inf_residual > 0.0
or length_violation
or not eval_success)
else 0.0,
7884 "length_violation": bool(length_violation),
7886 infeasibility_score, requirements, length_violation=length_violation
7891 if opt_state.get(
"valley_escape_tier", 0) > 0:
7892 layer1_logger.info(f
" *** Improvement found: Exiting Valley Escape Tier {opt_state['valley_escape_tier']} ***")
7893 opt_state[
"valley_escape_tier"] = 0
7896 opt_state[
"best_config"] = copy.deepcopy(config)
7897 opt_state[
"best_config_x"] = x_clipped.copy()
7898 opt_state[
"best_validation_tank_pa"] = (float(P_O_test), float(P_F_test))
7899 opt_state[
"best_lox_end_ratio"] = lox_start_ratio_hist
7900 opt_state[
"best_fuel_end_ratio"] = fuel_start_ratio_hist
7904 and not length_violation
7905 and final_pressures
is not None
7907 opt_state[
"best_pressures"] = final_pressures
7908 opt_state[
"best_results_for_validation"] = {
7911 "thrust_error": thrust_error,
7912 "of_error": of_error,
7913 "stability_score": stability_score,
7914 "stability_state": stability_state,
7915 "chugging_margin": chugging_margin,
7916 "acoustic_margin": acoustic_margin,
7917 "feed_margin": feed_margin,
7918 "stability_results": stability,
7921 f
" ✓ New best objective: {obj:.6f} "
7922 f
"(thrust_err: {thrust_error*100:.2f}%, O/F_err: {of_error*100:.2f}%, "
7923 f
"Cf: {Cf_actual:.3f}, stability: {stability_state}, score: {stability_score:.3f})"
7925 for handler
in layer1_logger.handlers:
7929 opt_state[
"objective_buffer"].append({
7930 "iteration": int(iteration),
7931 "objective": float(obj),
7932 "best_objective": float(opt_state.get(
"best_objective", obj)),
7937 should_report = (iteration % report_every_n == 0)
or opt_state.get(
'objective_satisfied',
False)
7938 if objective_callback
is not None and should_report:
7941 for buffered_entry
in opt_state[
"objective_buffer"]:
7945 buffered_entry[
"iteration"],
7946 buffered_entry[
"objective"],
7947 buffered_entry[
"best_objective"],
7950 opt_state[
"objective_buffer"] = []
7956 stability_acceptable = (
7958 (stability_score >= effective_min_score * 0.6)
and
7959 (chugging_margin >= effective_margin * 0.5)
and
7960 (acoustic_margin >= effective_margin * 0.5)
and
7961 (feed_margin >= effective_margin * 0.5)
7964 convergence_thrust_tol = thrust_tol * 2.0
7965 convergence_of_tol = 0.30
7969 thrust_error < convergence_thrust_tol
7970 and of_error < convergence_of_tol
7971 and stability_acceptable
7972 and np.isfinite(float(opt_state.get(
"best_objective", float(
"inf"))))
7973 and float(opt_state.get(
"best_objective", float(
"inf"))) <= obj_tolerance
7976 opt_state[
"converged"] = bool(opt_state.get(
"converged",
False)
or converged_now)
7981 layer1_logger.info(
"")
7982 layer1_logger.info(
"Starting optimization...")
7983 layer1_logger.info(
"")
7984 opt_state[
"iteration"] = 0
7985 opt_state[
"function_evaluations"] = 0
7987 class _ResultWrapper:
7988 def __init__(self, x, fun, success=True):
7989 self.x = np.asarray(x, dtype=float)
7990 self.fun = float(fun)
7991 self.success = success
8000 "CMA-ES (cma package) is required for Layer 1 optimization but not installed.\n"
8001 "Install with: pip install cma"
8005 layer1_logger.info(
"Using CMA-ES with random restart strategy for robust optimization.")
8006 update_progress(
"Layer 1: CMA-ES", 0.45,
"Running CMA-ES global solver...")
8009 target_fraction_of_range = 0.15
8013 sigma0 = float(np.percentile(span, 25) * target_fraction_of_range)
8014 if not np.isfinite(sigma0)
or sigma0 <= 0:
8019 cma_stds = np.ones_like(span)
8020 for i
in range(len(span)):
8023 desired_step = span[i] * target_fraction_of_range
8026 cma_stds[i] =
max(0.1, desired_step / sigma0)
if sigma0 > 0
else 1.0
8031 for idx_disc
in layer1_integer_dims:
8032 if 0 <= idx_disc < len(cma_stds):
8033 cma_stds[idx_disc] =
max(cma_stds[idx_disc], target_step_n / sigma0)
8038 layer1_logger.info(f
"Population size: {popsize}")
8050 requirements,
"layer1_generations_per_restart", 50.0)))
8051 if layer1_cma_restarts
is not None:
8052 num_restarts =
max(1, int(layer1_cma_restarts))
8054 num_restarts =
max(1, int(round(max_iterations / float(_gens_per_restart))))
8055 if layer1_smoke
and layer1_cma_restarts
is None:
8058 "Restart plan: %d generations total / %d per restart = %d restarts (%d evals at popsize %d)",
8059 max_iterations, _gens_per_restart, num_restarts, max_iterations * popsize, popsize,
8061 layer1_logger.info(f
"Using {num_restarts} CMA restart(s)")
8064 total_eval_budget = max_iterations
8069 opt_state[
"eval_budget"] =
max(1, int(max_iterations) * int(popsize) * int(num_restarts))
8071 best_x_global = x0_refined
8072 best_f_global = float(
'inf')
8079 debug_strict = getattr(config_obj.optimizer,
'debug_strict',
False)
if hasattr(config_obj,
'optimizer')
else False
8081 layer1_logger.info(f
"Using {num_workers} worker processes for parallel evaluation")
8091 _worker_cfg_src = copy.deepcopy(config_base)
8096 if os.environ.get(
"ED_L1_WORKER_THERMAL_OFF",
"0") ==
"1":
8097 if getattr(_worker_cfg_src,
"ablative_cooling",
None):
8098 _worker_cfg_src.ablative_cooling.enabled =
False
8099 if getattr(_worker_cfg_src,
"graphite_insert",
None):
8100 _worker_cfg_src.graphite_insert.enabled =
False
8102 bounds_array = np.column_stack([lower_bounds, upper_bounds])
8109 'layer1_contraction_half_angle_deg': requirements.get(
'layer1_contraction_half_angle_deg'),
8110 'layer1_min_Lcyl_over_D': requirements.get(
'layer1_min_Lcyl_over_D'),
8111 'layer1_max_element_pitch_m': requirements.get(
'layer1_max_element_pitch_m'),
8112 'target_thrust': target_thrust,
8113 'optimal_of': optimal_of,
8114 'P_ambient': target_P_exit,
8115 'max_lox_P_psi': max_lox_P_psi,
8116 'max_fuel_P_psi': max_fuel_P_psi,
8117 'TOTAL_WALL_THICKNESS_M': wall_total_m,
8118 'max_nozzle_exit': max_nozzle_exit,
8119 'max_chamber_od': max_chamber_od,
8121 'layer1_derive_expansion_ratio': derive_expansion_ratio,
8122 'layer1_derive_throat_from_thrust': derive_throat_from_thrust,
8123 'layer1_derive_max_iters': derive_max_iters,
8124 'layer1_derive_thrust_tol_rel': derive_thrust_tol,
8128 'derive_eps_min': derive_eps_range[0],
8129 'derive_eps_max': derive_eps_range[1],
8130 'derive_At_min': derive_At_range[0],
8131 'derive_At_max': derive_At_range[1],
8132 'derive_dF_min': derive_dF_range[0],
8133 'derive_dF_max': derive_dF_range[1],
8134 'layer1_derive_fuel_jet_from_of': derive_fuel_jet_from_of,
8135 'layer1_derive_tank_from_dp_ratio': derive_tank_from_dp,
8136 'layer1_dp_ratio_target': dp_ratio_target,
8137 'lock_P_T_min': lock_P_T_range[0],
8138 'lock_P_T_max': lock_P_T_range[1],
8139 'injector_type': l1_injector_type,
8142 'W_MOM': layer1_W_MOM,
8143 'W_DP': layer1_W_DP,
8144 'W_DP_CENTER': layer1_W_DP_CENTER,
8145 'W_DP_O': layer1_W_DP_O,
8146 'W_DP_F': layer1_W_DP_F,
8147 'W_DP_HIGH': layer1_W_DP_HIGH,
8150 'injector_dp_ratio_O_min': layer1_dp_o_aim[0],
8151 'injector_dp_ratio_O_max': layer1_dp_o_aim[1],
8152 'injector_dp_ratio_F_min': layer1_dp_f_aim[0],
8153 'injector_dp_ratio_F_max': layer1_dp_f_aim[1],
8154 'W_geom_ao_af_momentum': layer1_W_geom_ao_af,
8155 'rho_oxidizer': float(config_base.fluids[
"oxidizer"].density),
8156 'rho_fuel': float(config_base.fluids[
"fuel"].density),
8160 'impinging_momentum_R_min': layer1_R_mom_aim_lo,
8161 'impinging_momentum_R_max': layer1_R_mom_aim_hi,
8162 'layer1_momentum_log_deadband_rel': layer1_momentum_log_deadband_rel,
8163 'W_IMPINGING_ANGLE': layer1_W_impinging_angle,
8164 'W_IMPINGING_JET_ASYM': layer1_W_impinging_jet_asym,
8165 'layer1_impinging_jet_angle_max_asym_deg': layer1_impinging_jet_max_asym_deg,
8166 'layer1_exit_pressure_inside_quad_scale': layer1_exit_pressure_inside_quad,
8167 'layer1_impinging_angle_deg_min': layer1_impinging_angle_deg_lo,
8168 'layer1_impinging_angle_deg_max': layer1_impinging_angle_deg_hi,
8169 'W_CHAMBER_SHAPE': layer1_W_chamber_shape,
8170 'layer1_chamber_dt_ratio_min': layer1_chamber_dt_ratio_min,
8171 'layer1_chamber_dt_ratio_max': layer1_chamber_dt_ratio_max,
8172 'layer1_chamber_ld_ratio_min': layer1_chamber_ld_ratio_min,
8173 'layer1_chamber_ld_ratio_max': layer1_chamber_ld_ratio_max,
8174 'injector_dp_ratio_O_soft_floor': layer1_dp_o_soft_floor,
8175 'W_DP_O_FLOOR': layer1_W_DP_O_FLOOR,
8176 'W_SMD': layer1_W_SMD,
8177 'target_smd_microns': layer1_target_smd_microns,
8178 'layer1_smd_rel_tol': layer1_smd_rel_tol,
8179 'W_TANK_EQUAL': layer1_W_TANK_EQUAL,
8180 'layer1_tank_equal_scale_psi': layer1_tank_equal_scale_psi,
8181 'layer1_tank_equal_tol_psi': layer1_tank_equal_tol_psi,
8182 'layer1_of_deadband_rel': layer1_of_deadband_rel,
8183 'layer1_chamber_od_increment_in': layer1_chamber_od_increment_in,
8184 'layer1_chamber_od_snap_target': layer1_chamber_od_snap_target,
8185 'layer1_tank_equal_inband_frac': layer1_tank_equal_inband_frac,
8186 'layer1_lock_tank_pressures': layer1_lock_tank_pressures,
8187 'lock_P_F_min': lock_P_F_range[0],
8188 'lock_P_F_max': lock_P_F_range[1],
8189 'layer1_enforce_ring_geometry': layer1_enforce_ring_geometry,
8190 'layer1_injector_center_clear_dia_m': layer1_injector_center_clear_dia_m,
8191 'layer1_injector_min_web_m': layer1_injector_min_web_m,
8192 'layer1_injector_wall_clearance_m': layer1_injector_wall_clearance_m,
8193 'layer1_injector_spray_radius_frac': layer1_injector_spray_radius_frac,
8194 'layer1_injector_spray_radius_tol': layer1_injector_spray_radius_tol,
8195 'layer1_injector_min_face_incidence_deg': layer1_injector_min_face_incidence_deg,
8196 'layer1_injector_plate_thickness_m': layer1_injector_plate_thickness_m,
8197 'layer1_injector_counterbore_dia_m': layer1_injector_counterbore_dia_m,
8198 'layer1_injector_min_back_web_m': layer1_injector_min_back_web_m,
8199 'layer1_injector_land_ld_O': layer1_injector_land_ld_O,
8200 'layer1_injector_land_ld_F': layer1_injector_land_ld_F,
8201 'layer1_injector_back_keepout_dia_m': layer1_injector_back_keepout_dia_m,
8202 'layer1_injector_back_keepout_wall_m': layer1_injector_back_keepout_wall_m,
8203 'layer1_injector_face_contoured': layer1_injector_face_contoured,
8204 'layer1_injector_back_channels': layer1_injector_back_channels,
8205 'layer1_injector_exit_land': layer1_injector_exit_land,
8206 'layer1_injector_channel_width_m': layer1_injector_channel_width_m,
8207 'layer1_injector_channel_floor': layer1_injector_channel_floor,
8208 'layer1_injector_channel_inlets': layer1_injector_channel_inlets,
8209 'layer1_injector_plate_support': layer1_injector_plate_support,
8210 'layer1_injector_r_plate_m': layer1_injector_r_plate_m,
8211 'layer1_injector_r_port_m': layer1_injector_r_port_m,
8212 'layer1_injector_liner_thickness_m': layer1_injector_liner_thickness_m,
8213 'layer1_resultant_tilt_from_reach': layer1_resultant_tilt_from_reach,
8214 'layer1_resultant_tilt_reach_margin': layer1_resultant_tilt_reach_margin,
8215 'layer1_ring_order_fuel_outboard': layer1_ring_order_fuel_outboard,
8216 'layer1_momentum_wall_side_multiplier': layer1_momentum_wall_side_multiplier,
8217 'layer1_resultant_tilt_max_deg': layer1_resultant_tilt_max_deg,
8218 'layer1_resultant_tilt_scale_deg': layer1_resultant_tilt_scale_deg,
8219 'layer1_W_TILT': layer1_W_TILT,
8220 'layer1_resultant_tilt_outward_multiplier': layer1_resultant_tilt_outward_multiplier,
8221 'rho_O': float(config_base.fluids[
"oxidizer"].density),
8222 'rho_F': float(config_base.fluids[
"fuel"].density),
8223 'layer1_momentum_low_side_multiplier': layer1_momentum_low_side_multiplier,
8224 'layer1_momentum_band_width': layer1_momentum_band_width,
8225 'layer1_momentum_scale': layer1_momentum_scale,
8226 'layer1_derive_impingement_spacing': layer1_derive_impingement_spacing,
8227 'layer1_impingement_Ld_target': _ld_target,
8228 'derive_spacing_F_min': derive_spacing_F_range[0],
8229 'derive_spacing_F_max': derive_spacing_F_range[1],
8230 'layer1_impingement_Ld_min': layer1_impingement_Ld_min,
8231 'layer1_impingement_Ld_max': layer1_impingement_Ld_max,
8232 'W_IMP_GEOM': layer1_W_IMP_GEOM,
8233 'layer1_W_EXIT': layer1_W_EXIT,
8234 'layer1_exit_pressure_deadband_rel': layer1_exit_pressure_deadband_rel,
8235 'layer1_W_MASS': layer1_W_MASS,
8236 'layer1_W_ISP': layer1_W_ISP,
8237 'layer1_Isp_ref_s': layer1_Isp_ref_s,
8238 'target_burn_time': float(target_burn_time
or 0.0),
8239 'min_Lstar': min_Lstar,
8240 'max_Lstar': max_Lstar,
8241 'layer1_W_LSTAR': layer1_W_LSTAR,
8242 'layer1_Lstar_target_m': layer1_Lstar_target_m,
8243 'layer1_chamber_wall_density_kg_m3': layer1_chamber_wall_density,
8244 'layer1_chamber_mass_ref_kg': layer1_chamber_mass_ref_kg,
8245 'layer1_W_THRUST': layer1_W_THRUST_obj,
8246 'layer1_W_OF': layer1_W_OF_obj,
8247 'layer1_W_OF_low_MR_scale': layer1_W_OF_low_MR_scale,
8248 'layer1_W_OF_high_MR_scale': layer1_W_OF_high_MR_scale,
8249 'layer1_of_validation_tol': layer1_of_validation_tol,
8256 _pc_target_psi = requirements.get(
"target_chamber_pressure_psi")
8257 if _pc_target_psi
is not None and float(_pc_target_psi) > 0:
8258 constants_dict[
'layer1_target_Pc_pa'] = float(_pc_target_psi) * 6894.757293168
8264 constants_dict[
'layer1_W_PC'] =
_requirement_float(requirements,
"layer1_W_PC", 1.0e4)
8266 "Layer 1 chamber-pressure target: %.0f psi (%.3f MPa), W_PC=%g (exact L1). Throat sizes "
8267 "for thrust; tank pressure drives Pc to target.",
8268 float(_pc_target_psi), constants_dict[
'layer1_target_Pc_pa'] / 1e6,
8269 constants_dict[
'layer1_W_PC'],
8272 integer_dims = list(layer1_integer_dims)
8275 _l1_od_step_m = float(_od_increment_in) * 0.0254
if _od_increment_in > 0
else 0.0
8277 len(bounds), integer_dims, od_index=3, od_step_m=_l1_od_step_m, bounds=bounds)
8278 if _l1_int_vars
or _l1_od_step_m > 0:
8280 "Layer 1 discrete DOFs handed to CMA: integer_variables=%s%s",
8282 (
" + chamber-OD std floor %.4f m (%.2f in step)"
8283 % (_l1_minstd[3], _od_increment_in))
if _l1_od_step_m > 0
else "",
8287 optimizer_mode =
"cma"
8288 if hasattr(config_obj,
"optimizer")
and config_obj.optimizer:
8289 optimizer_mode = config_obj.optimizer.mode
8290 effective_mode = optimizer_mode
8292 effective_mode =
"cma"
8294 "Layer 1 smoke: forcing pure CMA-ES (config optimizer.mode=%r ignored for hybrid).",
8307 from engine
import accel
as _accel
8308 if _accel.enabled():
8310 except Exception
as e:
8311 layer1_logger.warning(
8312 "Accelerator warmup failed (%s); Layer-1 falls back to Python evaluation.", e)
8318 executor_context =
None
8321 executor_context = ProcessPoolExecutor(
8322 max_workers=num_workers,
8323 initializer=_init_worker,
8324 initargs=(config_dict, bounds_array, requirements, constants_dict, debug_strict),
8326 except (PermissionError, OSError)
as e:
8327 layer1_logger.warning(
8328 "ProcessPool unavailable (%s). Falling back to serial Layer-1 evaluation.",
8333 if executor_context
is None:
8334 _init_worker(config_dict, bounds_array, requirements, constants_dict, debug_strict)
8337 with executor_context
as executor:
8342 list(executor.map(_eval_candidate, [dummy_candidate] * num_workers))
8349 best_x_global = np.asarray(x0_refined, dtype=float).copy()
8350 best_f_global = float(
"inf")
8351 layer1_warmstart_applied =
False
8352 _ws_n_raw = requirements.get(
"layer1_cma_warmstart_trials")
8353 _ws_n = 16
if _ws_n_raw
is None else int(_ws_n_raw)
8355 _ws_n =
min(_ws_n, 4)
8360 lower_bounds=lower_bounds,
8361 upper_bounds=upper_bounds,
8363 integer_dims=integer_dims,
8368 eval_cache=eval_cache,
8369 make_cache_key_fn=_make_eval_cache_key,
8370 opt_state=opt_state,
8371 logger=layer1_logger,
8373 layer1_warmstart_applied =
True
8374 best_x_global = np.asarray(x0_refined, dtype=float).copy()
8375 if np.isfinite(_ws_best_f):
8376 best_f_global = float(_ws_best_f)
8378 opt_state, best_x_global, float(_ws_best_f), parent_objective=objective
8383 int(opt_state.get(
"function_evaluations", 0)),
8385 float(opt_state.get(
"best_objective", _ws_best_f)),
8389 _l1_iso_area = (4, (5, 8), (7, 10))
if l1_injector_type ==
"impinging" else None
8390 if effective_mode ==
"hybrid_cma_blocks" and hasattr(config_obj.optimizer,
"hybrid"):
8391 layer1_logger.info(
"Using Hybrid CMA + Block Re-optimization mode.")
8392 hybrid_config = config_obj.optimizer.hybrid
8402 total_budget_evals =
max(int(popsize), int(max_iterations) * int(popsize))
8403 opt_state[
"eval_budget"] = int(total_budget_evals)
8405 "Hybrid evaluation budget: %s (max(%s, max_iterations=%s x popsize=%s))",
8413 num_tracks = hybrid_config.num_tracks
8415 best_x_global = x0_refined
8416 best_f_global = float(
'inf')
8419 layer1_logger.info(f
"Running {num_tracks} independent hybrid tracks...")
8430 for track_i
in range(num_tracks):
8432 layer1_logger.info(f
"--- Track {track_i+1}/{num_tracks} ---")
8436 x0_track = x0_refined + rng.standard_normal(len(x0_refined)) * (0.1 * span)
8437 x0_track = np.clip(x0_track, lower_bounds, upper_bounds)
8439 x0_track = x0_refined
8443 list(zip(lower_bounds, upper_bounds)),
8446 total_budget=total_budget_evals // num_tracks,
8447 logger=layer1_logger,
8448 log_status_fn=log_status,
8449 update_progress_fn=update_progress,
8450 valley_escape_tracker=opt_state,
8453 integer_dims=integer_dims,
8455 od_step_m=_l1_od_step_m,
8456 fixed_variables=_l1_fixed_vars,
8457 eval_cache=eval_cache,
8458 make_cache_key_fn=_make_eval_cache_key,
8459 stop_event=stop_event,
8460 seed=int((layer1_seed_base + track_i * 7919) % (2 ** 31)),
8462 iso_area=_l1_iso_area,
8465 if t_f < best_f_global:
8468 layer1_logger.info(f
"Track {track_i+1} found new global best: {best_f_global:.5f}")
8474 list(zip(lower_bounds, upper_bounds)),
8477 total_budget=total_budget_evals,
8478 logger=layer1_logger,
8479 log_status_fn=log_status,
8480 update_progress_fn=update_progress,
8481 valley_escape_tracker=opt_state,
8484 integer_dims=integer_dims,
8486 od_step_m=_l1_od_step_m,
8487 fixed_variables=_l1_fixed_vars,
8488 eval_cache=eval_cache,
8489 make_cache_key_fn=_make_eval_cache_key,
8490 stop_event=stop_event,
8491 seed=int(layer1_seed_base),
8493 iso_area=_l1_iso_area,
8498 best_x_global = np.asarray(x0_refined, dtype=float).copy()
8499 if not layer1_warmstart_applied:
8500 best_f_global = float(
"inf")
8502 for restart_idx
in range(num_restarts):
8503 restart_name = f
"Restart {restart_idx + 1}/{num_restarts}" if num_restarts > 1
else "Main search"
8510 if restart_idx == 0:
8511 current_sigma_fraction = target_fraction_of_range
8512 x_start = x0_refined
8513 elif restart_idx % 2 == 1:
8520 current_sigma_fraction = target_fraction_of_range
8521 x_start = best_x_global + rng.standard_normal(len(best_x_global)) * (0.3 * span)
8522 layer1_logger.info(
" (global re-exploration: wide kick off incumbent best)")
8525 current_sigma_fraction = 0.15
8526 perturbation = rng.standard_normal(len(best_x_global)) * (sigma0 * 0.3)
8527 x_start = best_x_global + perturbation
8529 x_start = np.clip(x_start, lower_bounds, upper_bounds)
8531 current_sigma0 = float(np.percentile(span, 25) * current_sigma_fraction)
8532 if not np.isfinite(current_sigma0)
or current_sigma0 <= 0:
8533 current_sigma0 = 0.05
8536 if restart_idx == 0
and layer1_warmstart_applied:
8537 _s0_scale =
_requirement_float(requirements,
"layer1_cma_restart0_sigma_scale", 0.48)
8538 current_sigma0 *= float(np.clip(_s0_scale, 0.12, 1.0))
8540 current_cma_stds = np.ones_like(span)
8541 for i
in range(len(span)):
8543 desired_step = span[i] * current_sigma_fraction
8544 current_cma_stds[i] =
max(0.1, desired_step / current_sigma0)
if current_sigma0 > 0
else 1.0
8546 if current_sigma0 > 0:
8547 for idx_disc
in layer1_integer_dims:
8548 if 0 <= idx_disc < len(current_cma_stds):
8549 current_cma_stds[idx_disc] =
max(
8550 current_cma_stds[idx_disc], 1.0 / current_sigma0
8553 layer1_logger.info(f
"")
8554 layer1_logger.info(f
"Starting {restart_name} (sigma: {current_sigma_fraction*100:.0f}% of range)...")
8559 iter_budget =
max(1, total_eval_budget // num_restarts)
8562 _cma_restart_seed = layer1_seed_base + int(restart_idx) * 1_000_003
8564 "bounds": [lower_bounds.tolist(), upper_bounds.tolist()],
8566 "maxiter": iter_budget,
8570 "integer_variables": _l1_int_vars,
8571 "minstd": _l1_minstd,
8574 "tolstagnation": 50,
8576 "seed": _cma_restart_seed,
8580 es = cma.CMAEvolutionStrategy(x_start.tolist(), current_sigma0, cma_options)
8581 while not es.stop():
8586 target_tier = _check_valley_escape_tier()
8587 if target_tier > opt_state[
"valley_escape_tier"]
and opt_state[
"function_evaluations"] > opt_state[
"cooldown_until"]:
8589 tier_names = [
"None",
"Mild",
"Medium",
"Full"]
8590 boost_factors = [1.0, 1.5, 2.0, 3.0]
8591 current_factor = boost_factors[opt_state[
"valley_escape_tier"]]
8592 target_factor = boost_factors[target_tier]
8593 relative_boost = target_factor / current_factor
8595 old_sigma = es.sigma
8596 es.sigma *= relative_boost
8599 max_sigma = 0.3 * float(np.percentile(span, 25))
8600 if es.sigma > max_sigma:
8601 es.sigma = max_sigma
8604 opt_state[
"valley_escape_tier"] = target_tier
8605 opt_state[
"cooldown_until"] = opt_state[
"function_evaluations"] + 1000
8607 layer1_logger.info(
"!" * 20)
8608 layer1_logger.info(f
"VALLEY ESCAPE TRIGGERED (Tier: {tier_names[target_tier]})")
8609 layer1_logger.info(f
"Evals: {opt_state['function_evaluations']}, Best: {opt_state['best_objective']:.2f}, Stagnation: {opt_state['function_evaluations'] - opt_state['last_best_eval']}")
8610 layer1_logger.info(f
"Sigma: {old_sigma:.6f} -> {es.sigma:.6f} (Boost: {relative_boost:.2f}x)")
8611 layer1_logger.info(
"!" * 20)
8612 for handler
in layer1_logger.handlers:
8617 candidates = es.ask()
8621 candidates_snapped = [
8628 cache_keys = [_make_eval_cache_key(c)
for c
in candidates_snapped]
8629 uncached_indices = []
8630 uncached_candidates = []
8631 values = [
None] * len(candidates)
8632 gen_meta: List[Optional[Dict[str, Any]]] = [
None] * len(candidates)
8634 for i, key
in enumerate(cache_keys):
8635 if key
in eval_cache:
8636 cached_val = eval_cache[key]
8639 if isinstance(cached_val, dict)
and 'value' in cached_val:
8640 values[i] = float(cached_val[
'value'])
8641 if cached_val.get(
"worker_payload"):
8642 gen_meta[i] = dict(cached_val[
"worker_payload"])
8643 elif isinstance(cached_val, (int, float)):
8645 values[i] = float(cached_val)
8648 uncached_indices.append(i)
8649 uncached_candidates.append(candidates_snapped[i])
8651 uncached_indices.append(i)
8652 uncached_candidates.append(candidates_snapped[i])
8655 if uncached_candidates:
8656 chunksize =
max(1, len(uncached_candidates) // (num_workers * 4))
8657 results = list(executor.map(_eval_candidate, uncached_candidates, chunksize=chunksize))
8660 for idx, res
in zip(uncached_indices, results):
8661 obj_val = res[
'value']
8662 values[idx] = obj_val
8666 eval_cache[cache_keys[idx]] = {
8668 'success': res.get(
'success',
True),
8671 'success': res.get(
'success',
True),
8672 'thrust_error': res.get(
'thrust_error'),
8673 'of_error': res.get(
'of_error'),
8674 'momentum_term': res.get(
'momentum_term'),
8675 'inf_residual': res.get(
'inf_residual'),
8677 'MR': res.get(
'MR'),
8678 'full_results': res.get(
'full_results'),
8684 opt_state[
'function_evaluations'] += 1
8687 if not res[
'success']:
8689 reason = res.get(
'error_type',
'Unknown')
8690 opt_state[
'num_failures'] = opt_state.get(
'num_failures', 0) + 1
8693 if 'fail_counts' not in opt_state:
8694 opt_state[
'fail_counts'] = {}
8695 opt_state[
'fail_counts'][reason] = opt_state[
'fail_counts'].get(reason, 0) + 1
8698 if opt_state[
'function_evaluations'] % 100 == 0:
8700 sorted_fails = sorted(opt_state[
'fail_counts'].items(), key=
lambda item: item[1], reverse=
True)
8701 top_fails = sorted_fails[:3]
8702 fail_msg =
", ".join([f
"{k}: {v}" for k, v
in top_fails])
8703 layer1_logger.warning(f
"Recent Failures (Top 3): {fail_msg}")
8705 layer1_logger.warning(f
"Eval failed: {reason}")
8709 es.tell(candidates, values)
8712 float(v)
if v
is not None and np.isfinite(float(v))
else float(
"inf")
8715 if _fv
and min(_fv) < float(
"inf"):
8716 _gen_best = float(
min(_fv))
8717 _gi = int(np.argmin(np.asarray(_fv, dtype=float)))
8720 candidates_snapped[_gi],
8722 parent_objective=objective,
8727 int(opt_state.get(
"function_evaluations", 0)),
8729 float(opt_state.get(
"best_objective", _gen_best)),
8732 iter_idx =
max(1, es.countiter)
8733 overall_progress = restart_idx / num_restarts
8734 restart_progress = iter_idx / iter_budget / num_restarts
8735 progress = 0.10 + 0.35 * (overall_progress + restart_progress)
8736 update_progress(
"Layer 1: CMA-ES", progress,
8737 f
"{restart_name} - iteration {iter_idx}/{iter_budget}")
8739 cma_result = es.result
8740 restart_best_x = np.asarray(cma_result.xbest, dtype=float)
8741 restart_best_f = float(cma_result.fbest)
8743 layer1_logger.info(f
"{restart_name} finished. Final obj: {restart_best_f:.6f}")
8745 if restart_best_f < best_f_global:
8746 best_f_global = restart_best_f
8747 best_x_global = restart_best_x
8748 layer1_logger.info(f
" ✓ New global best: {best_f_global:.6f}")
8750 if best_f_global < obj_tolerance:
8751 layer1_logger.info(f
"Found excellent solution, skipping restarts")
8754 for handler
in layer1_logger.handlers:
8757 except Exception
as e:
8758 layer1_logger.error(f
"{restart_name} failed: {e}")
8759 if restart_idx == 0:
8765 x0_refined = best_x_global
8766 best_fun = best_f_global
8767 layer1_logger.info(f
"")
8768 layer1_logger.info(f
"CMA-ES complete. Global best objective: {best_fun:.6f}")
8769 layer1_logger.info(f
"Total function evaluations: {opt_state['function_evaluations']}")
8770 log_status(
"Layer 1", f
"CMA-ES complete: {num_restarts} restart(s), obj={best_fun:.3f}, refining with L-BFGS-B...")
8771 for handler
in layer1_logger.handlers:
8775 opt_state[
"function_evaluations"] = 0
8784 maxfun_capped =
min(max_iterations * 5, 2500)
8786 opt_state[
"layer1_lbfgs_phase"] =
True
8787 update_progress(
"Layer 1: Local Refinement", 0.47, f
"Refining with L-BFGS-B (max {maxfun_capped} func evals)...")
8788 layer1_logger.info(
"Phase 2: Local refinement (L-BFGS-B)...")
8791 lbfgs_result = minimize(
8797 'maxiter': max_iterations,
8798 'maxfun': maxfun_capped,
8800 'gtol': layer1_lbfgs_gtol,
8805 if layer1_lbfgs_second_pass
and getattr(lbfgs_result,
"success",
True):
8806 _gt2 = float(
max(1e-15,
min(1e-12, layer1_lbfgs_gtol * 0.05)))
8808 lbfgs_result2 = minimize(
8810 np.asarray(lbfgs_result.x, dtype=float),
8814 "maxiter": max_iterations,
8815 "maxfun":
min(800, maxfun_capped),
8822 np.isfinite(getattr(lbfgs_result2,
"fun", float(
"nan")))
8826 or float(lbfgs_result2.fun) <= float(lbfgs_result.fun)
8829 lbfgs_result = lbfgs_result2
8831 "L-BFGS-B second pass improved feasible objective to %.6g (gtol=%.1e)",
8832 float(lbfgs_result.fun),
8835 elif np.isfinite(getattr(lbfgs_result2,
"fun", float(
"nan"))):
8837 "L-BFGS-B second pass rejected: fun=%.6g (infeasible or worse than %.6g)",
8838 float(lbfgs_result2.fun),
8839 float(lbfgs_result.fun),
8841 except Exception
as _e2:
8842 layer1_logger.debug(
"L-BFGS-B second pass skipped: %s", _e2)
8843 layer1_logger.info(
"")
8844 layer1_logger.info(
"Optimization completed")
8845 layer1_logger.info(f
"Success: {lbfgs_result.success}")
8846 _report_fun = float(getattr(lbfgs_result,
"fun", float(
"nan")))
8848 _report_fun = float(opt_state.get(
"best_objective", _report_fun))
8850 "L-BFGS-B terminal iterate was infeasible (obj≈1e6); reporting incumbent best=%.6g",
8853 layer1_logger.info(f
"Final objective value: {_report_fun:.6f}")
8854 layer1_logger.info(f
"Iterations: {lbfgs_result.nit if hasattr(lbfgs_result, 'nit') else 'N/A'}")
8855 layer1_logger.info(f
"Function evaluations: {lbfgs_result.nfev if hasattr(lbfgs_result, 'nfev') else 'N/A'}")
8856 layer1_logger.info(
"")
8857 result = lbfgs_result
8858 except Exception
as e:
8859 layer1_logger.error(f
"L-BFGS-B error: {e}")
8860 log_status(
"Layer 1 Warning", f
"L-BFGS-B error: {e}, using best result found")
8862 if 'best_fun' in locals():
8864 result = _ResultWrapper(x0_refined, best_fun)
8867 best_x = opt_state.get(
'best_x', x0)
8868 best_obj = opt_state.get(
'best_objective', float(
'inf'))
8869 result = _ResultWrapper(best_x, best_obj)
8871 opt_state.pop(
"layer1_lbfgs_phase",
None)
8883 _incumbent_obj = float(opt_state.get(
"best_objective", float(
"inf")))
8884 _incumbent_x = opt_state.get(
"best_x")
8886 _incumbent_x
is not None
8887 and np.isfinite(_incumbent_obj)
8888 and _incumbent_obj < float(getattr(result,
"fun", float(
"inf")))
8891 "Reconciling final result to incumbent best: result.fun=%.6g -> best_objective=%.6g "
8892 "(L-BFGS-B terminal iterate was worse than the tracked global best)",
8893 float(getattr(result,
"fun", float(
"inf"))),
8896 result = _ResultWrapper(np.asarray(_incumbent_x, dtype=float).copy(), _incumbent_obj)
8901 _final_inc_x = opt_state.get(
"best_x")
8902 if _final_inc_x
is not None:
8905 "Final incumbent parent re-eval (obj=%.6g) for validation snapshot and history",
8906 float(opt_state.get(
"best_objective", float(
"nan"))),
8908 objective(np.asarray(_final_inc_x, dtype=float))
8909 except Exception
as _final_exc:
8910 layer1_logger.warning(
"Final incumbent parent re-eval failed: %s", _final_exc)
8912 validation_tank_pa: Optional[Tuple[float, float]] =
None
8919 _rb_x = getattr(result,
"x",
None)
8920 _rb_src =
"result.x"
8921 if _rb_x
is None or len(_rb_x) != len(bounds)
or not np.all(np.isfinite(np.asarray(_rb_x, dtype=float))):
8922 _rb_x = opt_state.get(
"best_x"); _rb_src =
"opt_state.best_x"
8923 opt_state[
"final_rebuild_diag"] = {
"fired":
False,
"src": _rb_src}
8924 if _rb_x
is not None and len(_rb_x) == len(bounds):
8925 layer1_logger.info(
"Final rebuild: from %s (len %d)", _rb_src, len(_rb_x))
8927 np.asarray(_rb_x, dtype=float), config_base, constants_dict, logger=layer1_logger,
8928 diag=opt_state[
"final_rebuild_diag"])
8929 opt_state[
"final_rebuild_diag"][
"fired"] =
True
8930 opt_state[
"final_rebuild_diag"][
"returned"] = _rebuilt
is not None
8931 if _rebuilt
is None:
8932 layer1_logger.warning(
"Final rebuild returned None -- falling back to snapshot path")
8934 layer1_logger.warning(
8935 "Final rebuild SKIPPED: no usable vector (result.x=%s, best_x=%s, expected len %d)",
8936 None if getattr(result,
"x",
None)
is None else len(result.x),
8937 None if opt_state.get(
"best_x")
is None else len(opt_state[
"best_x"]), len(bounds))
8938 if _rebuilt
is not None:
8939 optimized_config, _rb_x_solved, _rb_result = _rebuilt
8944 _rb_result, np.asarray(_rb_x_solved, dtype=float).copy(), requirements,
8945 constants_dict, cea_cache=_l1_cea_cache)
8946 _bd = {
"objective": float(_rbT[
"objective"]), **_rbT[
"terms"]}
8947 _bd[
"penalty_sum"] = float(sum(_rbT[
"terms"].values()))
8948 _bd[
"penalty_unaccounted"] = _bd[
"objective"] - _bd[
"penalty_sum"]
8949 for _k
in (
"injector_dp_ratio_O",
"injector_dp_ratio_F",
"chamber_D_over_Dt",
8950 "chamber_L_over_D",
"Lstar_m",
"Lstar_target_m",
"chamber_mass_kg",
8951 "Isp",
"Isp_merit_s",
"Isp_ref_s",
"propellant_mass_kg",
8952 "L_face_to_throat",
"L_engine",
"momentum_ratio_R",
"resultant_tilt_deg",
8953 "resultant_tilt_breakeven_deg",
"effective_smd_microns",
8954 "impingement_angle_deg",
"infeasibility_score",
"length_violation"):
8955 _bd[_k] = _rbT.get(_k)
8956 _bd[
"source"] =
"emitted design (final rebuild)"
8957 opt_state[
"best_objective_breakdown"] = _bd
8958 except Exception
as _bd_exc:
8959 layer1_logger.warning(
"Breakdown of the emitted design failed: %s", _bd_exc)
8961 opt_state[
"best_config"] = copy.deepcopy(optimized_config)
8962 opt_state[
"best_config_x"] = np.asarray(_rb_x_solved, dtype=float).copy()
8963 opt_state[
"last_eval_config"] = copy.deepcopy(optimized_config)
8964 opt_state[
"last_eval_config_x"] = np.asarray(_rb_x_solved, dtype=float).copy()
8965 final_lox_end_ratio = opt_state.get(
"best_lox_end_ratio", 0.7)
8966 final_fuel_end_ratio = opt_state.get(
"best_fuel_end_ratio", 0.7)
8967 validation_tank_pa = (
8968 float(_rb_x_solved[idx_P_O]) * 6894.757293168,
8969 float(_rb_x_solved[idx_P_F]) * 6894.757293168,
8971 elif opt_state[
"best_config"]
is not None:
8972 optimized_config = copy.deepcopy(opt_state[
"best_config"])
8973 final_lox_end_ratio = opt_state.get(
"best_lox_end_ratio", 0.7)
8974 final_fuel_end_ratio = opt_state.get(
"best_fuel_end_ratio", 0.7)
8975 validation_tank_pa = opt_state.get(
"best_validation_tank_pa")
8976 elif opt_state.get(
"last_eval_config")
is not None:
8977 optimized_config = copy.deepcopy(opt_state[
"last_eval_config"])
8978 final_lox_end_ratio = opt_state.get(
"best_lox_end_ratio", 0.7)
8979 final_fuel_end_ratio = opt_state.get(
"best_fuel_end_ratio", 0.7)
8980 validation_tank_pa = opt_state.get(
"last_eval_validation_tank_pa")
8982 _rx = getattr(result,
"x",
None)
8983 _rx = np.asarray(_rx
if _rx
is not None else x0, dtype=float).copy()
8984 optimized_config, P_O_final_psi, P_F_final_psi = apply_x_to_config(_rx, config_base)
8985 final_lox_end_ratio = P_O_final_psi / max_lox_P_psi
if max_lox_P_psi > 0
else 0.7
8986 final_fuel_end_ratio = P_F_final_psi / max_fuel_P_psi
if max_fuel_P_psi > 0
else 0.7
8989 if optimized_config.lox_tank
is None:
8991 optimized_config.lox_tank = LOXTankConfig(lox_h=0.5, lox_radius=0.1, ox_tank_pos=1.0)
8992 if optimized_config.fuel_tank
is None:
8994 optimized_config.fuel_tank = FuelTankConfig(rp1_h=0.5, rp1_radius=0.1, fuel_tank_pos=0.5)
8997 best_x = opt_state.get(
"best_x", result.x
if hasattr(result,
'x')
else x0)
9000 x_for_tank_psi = opt_state.get(
"best_config_x")
9001 if x_for_tank_psi
is None or len(x_for_tank_psi) <= idx_P_F:
9002 x_for_tank_psi = opt_state.get(
"last_eval_config_x")
9003 if x_for_tank_psi
is None or len(x_for_tank_psi) <= idx_P_F:
9004 x_for_tank_psi = best_x
9005 if x_for_tank_psi
is not None and len(x_for_tank_psi) > idx_P_F:
9006 P_O_start_optimized_psi = float(
9007 np.clip(x_for_tank_psi[idx_P_O], bounds[idx_P_O][0], bounds[idx_P_O][1])
9013 P_F_start_optimized_psi = float(
9014 np.clip(x_for_tank_psi[idx_P_F], lower_bounds_eff[idx_P_F], upper_bounds_eff[idx_P_F])
9016 optimized_config.lox_tank.initial_pressure_psi = P_O_start_optimized_psi
9017 optimized_config.fuel_tank.initial_pressure_psi = P_F_start_optimized_psi
9022 if _rebuilt
is None:
9025 P_O_start_optimized_psi * 6894.757293168,
9026 P_F_start_optimized_psi * 6894.757293168,
9028 logger=layer1_logger,
9031 optimized_config.lox_tank.initial_pressure_psi = max_lox_P_psi * 0.8
9032 optimized_config.fuel_tank.initial_pressure_psi = max_fuel_P_psi * 0.8
9035 if hasattr(optimized_config,
'injector')
and optimized_config.injector.type ==
"pintle":
9036 if hasattr(optimized_config.injector.geometry,
'lox'):
9037 optimized_config.injector.geometry.lox.theta_orifice = 90.0
9039 iteration_history = opt_state[
"history"]
9041 def _validation_evaluate_boost(
9042 runner: PintleEngineRunner,
9045 ) -> Tuple[Dict[str, Any], float]:
9046 """Evaluate the emitted design at its OWN tank pressures.
9048 If it does not close there, find -- for the report only, never the verdict -- the
9049 smallest tank boost inside the tank limits at which it would. Any scale other than 1.0
9050 makes the design invalid (see validation_replay_ok below).
9052 po0, pf0 = float(po), float(pf)
9054 return runner.evaluate(po0, pf0, P_ambient=target_P_exit, silent=
True), 1.0
9055 except ValueError
as exc:
9056 if "Supply < Demand" not in str(exc)
and "Insufficient mass flow" not in str(exc):
9059 _cap =
min(max_lox_P_psi * psi_to_Pa /
max(po0, 1.0),
9060 max_fuel_P_psi * psi_to_Pa /
max(pf0, 1.0))
9061 for scale
in (1.03, 1.06, 1.09, 1.14, 1.20, 1.28, 1.40, 1.55, 1.72):
9065 return runner.evaluate(po0 * scale, pf0 * scale, P_ambient=target_P_exit,
9067 except ValueError
as exc:
9068 if "Supply < Demand" not in str(exc)
and "Insufficient mass flow" not in str(exc):
9072 def _validation_evaluate_or_bundle(
9073 runner: PintleEngineRunner,
9076 ) -> Tuple[Dict[str, Any], float]:
9077 """Replay validation evaluate; fall back to the last optimizer ``evaluate`` payload if the
9078 solver diverges -- for the report only; the design is then invalid.
9080 The chamber solver says it did not close in two ways: a ``ValueError`` from the root
9081 find, and a ``RuntimeError`` from its own checks on what it found ("Solution
9082 validation failed: Isp is non-positive"). Both mean the emitted design does not
9083 close at its tank pressures; only the first used to be caught, so a design whose
9084 replay landed on a non-physical root crashed the optimizer instead of being reported
9085 invalid -- on CI's runners, where the smoke search's last candidate is one of them.
9088 return _validation_evaluate_boost(runner, po, pf)
9089 except (ValueError, RuntimeError)
as exc:
9090 if isinstance(exc, RuntimeError)
and not str(exc).startswith(_CHAMBER_SOLVER_REFUSALS):
9092 b = opt_state.get(
"last_good_eval_bundle")
9093 if not isinstance(b, dict)
or not isinstance(b.get(
"results"), dict):
9095 layer1_logger.warning(
9096 "Layer 1 validation replay failed (%s); reporting the last in-loop evaluate "
9097 "payload, and the design is INVALID.", exc
9099 perf = copy.deepcopy(b[
"results"])
9100 perf[
"layer1_validation_used_last_good_bundle"] =
True
9101 perf[
"layer1_validation_replay_error"] = str(exc)[:200]
9105 update_progress(
"Layer 1: Validation", 0.52,
"Validating optimized configuration...")
9107 best_x = opt_state.get(
"best_x", result.x
if hasattr(result,
'x')
else x0)
9108 validation_pressure_scale = 1.0
9112 if validation_tank_pa
is not None and len(validation_tank_pa) == 2:
9113 P_O_initial = float(validation_tank_pa[0])
9114 P_F_initial = float(validation_tank_pa[1])
9117 P_O_initial = float(optimized_config.lox_tank.initial_pressure_psi) * psi_to_Pa
9118 P_F_initial = float(optimized_config.fuel_tank.initial_pressure_psi) * psi_to_Pa
9120 if best_x
is not None and len(best_x) > idx_P_O:
9121 P_O_initial = float(np.clip(best_x[idx_P_O], bounds[idx_P_O][0], bounds[idx_P_O][1])) * psi_to_Pa
9123 P_O_initial = max_lox_P_psi * psi_to_Pa * 0.95
9124 if best_x
is not None and len(best_x) > idx_P_F:
9125 P_F_initial = float(np.clip(best_x[idx_P_F], bounds[idx_P_F][0], bounds[idx_P_F][1])) * psi_to_Pa
9127 P_F_initial = max_fuel_P_psi * psi_to_Pa * 0.95
9133 optimized_config_runner = copy.deepcopy(optimized_config)
9134 if os.environ.get(
"ED_L1_VALIDATE_THERMAL_OFF",
"0") ==
"1":
9135 if getattr(optimized_config_runner,
"ablative_cooling",
None):
9136 optimized_config_runner.ablative_cooling.enabled =
False
9137 if getattr(optimized_config_runner,
"graphite_insert",
None):
9138 optimized_config_runner.graphite_insert.enabled =
False
9143 infeasible_reason =
None
9147 if infeasible_reason
and log_status:
9148 log_status(
"warning", infeasible_reason)
9150 infeasible_reason =
None
9153 if "best_results_for_validation" in opt_state
and opt_state[
"best_results_for_validation"]
is not None:
9154 stored_results = opt_state[
"best_results_for_validation"]
9157 initial_performance = {
9158 "F": stored_results[
"F"],
9159 "MR": stored_results[
"MR"],
9160 "stability_results": stored_results.get(
"stability_results", {}),
9162 initial_thrust_error = stored_results[
"thrust_error"]
9163 initial_MR_error = stored_results[
"of_error"]
9164 stored_stability_score = stored_results.get(
"stability_score",
None)
9165 stored_stability_state = stored_results.get(
"stability_state",
None)
9166 log_status(
"Layer 1 Validation", f
"Using stored validation results: Thrust err {initial_thrust_error*100:.2f}%, O/F err {initial_MR_error*100:.2f}%")
9171 eval_results, validation_pressure_scale = _validation_evaluate_or_bundle(
9172 optimized_runner, P_O_initial, P_F_initial
9175 if "P_exit" in eval_results:
9176 initial_performance[
"P_exit"] = eval_results[
"P_exit"]
9177 if "Cf" in eval_results
or "Cf_actual" in eval_results:
9178 initial_performance[
"Cf"] = eval_results.get(
"Cf_actual", eval_results.get(
"Cf"))
9179 initial_performance[
"Cf_actual"] = initial_performance[
"Cf"]
9181 if "Isp" in eval_results:
9182 initial_performance[
"Isp"] = eval_results[
"Isp"]
9183 if "Pc" in eval_results:
9184 initial_performance[
"Pc"] = eval_results[
"Pc"]
9186 if "mdot_total" in eval_results:
9187 initial_performance[
"mdot_total"] = eval_results[
"mdot_total"]
9188 if "mdot_O" in eval_results:
9189 initial_performance[
"mdot_O"] = eval_results[
"mdot_O"]
9190 if "mdot_F" in eval_results:
9191 initial_performance[
"mdot_F"] = eval_results[
"mdot_F"]
9192 if "eta_cstar" in eval_results:
9193 initial_performance[
"eta_cstar"] = eval_results[
"eta_cstar"]
9194 if "cstar_actual" in eval_results:
9195 initial_performance[
"cstar_actual"] = eval_results[
"cstar_actual"]
9196 if "cstar_ideal" in eval_results:
9197 initial_performance[
"cstar_ideal"] = eval_results[
"cstar_ideal"]
9199 if "F" in eval_results:
9200 initial_performance[
"F"] = eval_results[
"F"]
9201 if "MR" in eval_results:
9202 initial_performance[
"MR"] = eval_results[
"MR"]
9203 if eval_results.get(
"stability_results"):
9204 initial_performance[
"stability_results"] = eval_results[
"stability_results"]
9206 if "chamber_intrinsics" in eval_results
and eval_results[
"chamber_intrinsics"]:
9207 initial_performance[
"chamber_intrinsics"] = eval_results[
"chamber_intrinsics"]
9209 for _vk
in (
"layer1_validation_used_last_good_bundle",
"layer1_validation_replay_error"):
9210 if _vk
in eval_results:
9211 initial_performance[_vk] = eval_results[_vk]
9212 except Exception
as _replay_exc:
9213 initial_performance[
"layer1_validation_replay_error"] = str(_replay_exc)[:200]
9215 F_val = initial_performance.get(
"F", 0)
9216 Pc_val = initial_performance.get(
"Pc", 0)
9218 if A_throat_val
and A_throat_val > 0
and Pc_val > 0:
9219 Cf_calculated = F_val / (Pc_val * A_throat_val)
9220 initial_performance[
"Cf_actual"] = Cf_calculated
9221 initial_performance[
"Cf"] = Cf_calculated
9223 initial_performance, validation_pressure_scale = _validation_evaluate_or_bundle(
9224 optimized_runner, P_O_initial, P_F_initial
9226 initial_thrust_error = (
9227 abs(initial_performance.get(
"F", 0) - target_thrust) / target_thrust
if target_thrust > 0
else 1.0
9229 initial_MR_error = (
9230 abs(initial_performance.get(
"MR", 0) - optimal_of) / optimal_of
if optimal_of > 0
else 1.0
9234 if "Cf" not in initial_performance
and "Cf_actual" not in initial_performance:
9235 F_val = initial_performance.get(
"F", 0)
9236 Pc_val = initial_performance.get(
"Pc", 0)
9238 if A_throat_val
and A_throat_val > 0
and Pc_val > 0:
9239 Cf_calculated = F_val / (Pc_val * A_throat_val)
9240 initial_performance[
"Cf_actual"] = Cf_calculated
9241 initial_performance[
"Cf"] = Cf_calculated
9243 if validation_pressure_scale
is not None and np.isfinite(float(validation_pressure_scale)):
9244 sc = float(validation_pressure_scale)
9245 initial_performance[
"layer1_validation_tank_pressure_scale"] = sc
9247 _po_nom_psi =
min(float(P_O_initial) / float(psi_to_Pa), max_lox_P_psi)
9248 _pf_nom_psi =
min(float(P_F_initial) / float(psi_to_Pa), max_fuel_P_psi)
9253 optimized_config.lox_tank.initial_pressure_psi = _po_nom_psi
9254 optimized_config.fuel_tank.initial_pressure_psi = _pf_nom_psi
9255 if abs(sc - 1.0) > 1e-6:
9256 initial_performance[
"layer1_validation_replay_boosted_tanks"] =
True
9257 initial_performance[
"P_O_start_psi_validation_replay"] = float(_po_nom_psi * sc)
9258 initial_performance[
"P_F_start_psi_validation_replay"] = float(_pf_nom_psi * sc)
9261 if getattr(getattr(optimized_config,
"injector",
None),
"type",
None) ==
"impinging":
9262 _diag_vals = initial_performance.get(
"diagnostics")
9263 initial_performance[
"momentum_balance_penalty"] = 0.0
9264 if isinstance(_diag_vals, dict):
9265 _mr_promote = _diag_vals.get(
"momentum_ratio_R")
9266 if _mr_promote
is not None and np.isfinite(_mr_promote)
and float(_mr_promote) > 0:
9267 initial_performance[
"momentum_ratio_R"] = float(_mr_promote)
9268 initial_performance[
"momentum_balance_penalty"] = float(
9272 r_band_lo=layer1_R_mom_aim_lo,
9273 r_band_hi=layer1_R_mom_aim_hi,
9278 (
"v_O_bulk", float),
9279 (
"v_F_bulk", float),
9282 (
"rho_O_momentum", float),
9283 (
"rho_F_momentum", float),
9287 _rv = _diag_vals.get(_rk)
9291 if np.isfinite(_rf):
9292 initial_performance[_rk] = _rf
9293 except (TypeError, ValueError):
9295 for _ik
in (
"momentum_ratio_n_elements_O",
"momentum_ratio_n_elements_F"):
9296 _iv = _diag_vals.get(_ik)
9299 initial_performance[_ik] = int(_iv)
9300 except (TypeError, ValueError):
9303 ig = optimized_config.injector.geometry
9304 rho_o_p = float(optimized_config.fluids[
"oxidizer"].density)
9305 rho_f_p = float(optimized_config.fluids[
"fuel"].density)
9306 Aog_p = float(ig.oxidizer.n_elements * np.pi * (ig.oxidizer.d_jet / 2.0) ** 2)
9307 Afg_p = float(ig.fuel.n_elements * np.pi * (ig.fuel.d_jet / 2.0) ** 2)
9310 Aog_p, Afg_p, optimal_of, rho_o_p, rho_f_p
9312 initial_performance[
"geom_ao_af"] = float(ao_r)
9313 initial_performance[
"expected_ao_af_for_R1"] = float(exp_r)
9314 if np.isfinite(exp_r)
and float(exp_r) != 0.0:
9315 initial_performance[
"geom_ao_af_rel_error"] = float((float(ao_r) - float(exp_r)) / float(exp_r))
9316 initial_performance[
"geom_ao_af_momentum_hint_sq"] = float(_sq_g)
9321 initial_performance[
"injector_dp_out_of_range"] =
False
9322 _pc_prom = initial_performance.get(
"Pc")
9323 if _pc_prom
is not None and np.isfinite(float(_pc_prom))
and float(_pc_prom) > 0:
9324 _ro_p, _rf_p = injector_dp_ratios_from_eval_result(float(_pc_prom), initial_performance)
9325 if _ro_p
is not None:
9326 initial_performance[
"injector_dp_ratio_O"] = float(_ro_p)
9327 if _rf_p
is not None:
9328 initial_performance[
"injector_dp_ratio_F"] = float(_rf_p)
9329 initial_performance[
"injector_dp_penalty"] = float(
9330 injector_dp_ratio_penalty_weighted(
9335 o_band=layer1_dp_o_aim,
9336 f_band=layer1_dp_f_aim,
9337 w_dp_o=layer1_W_DP_O,
9338 w_dp_f=layer1_W_DP_F,
9339 o_soft_floor=layer1_dp_o_soft_floor,
9340 w_dp_o_floor=layer1_W_DP_O_FLOOR,
9343 lo_o, hi_o = layer1_dp_o_band
9344 lo_f, hi_f = layer1_dp_f_band
9346 if injector_dp_ratio_within_gate(_ro_p, float(lo_o), float(hi_o))
is False:
9348 if injector_dp_ratio_within_gate(_rf_p, float(lo_f), float(hi_f))
is False:
9350 initial_performance[
"injector_dp_out_of_range"] = bool(_oor)
9353 _Fv_gate = initial_performance.get(
"F")
9354 _MRv_gate = initial_performance.get(
"MR")
9355 if target_thrust > 0
and _Fv_gate
is not None and np.isfinite(float(_Fv_gate)):
9356 initial_thrust_error = abs(float(_Fv_gate) - float(target_thrust)) / float(target_thrust)
9358 initial_thrust_error = 1.0
9359 if optimal_of > 0
and _MRv_gate
is not None and np.isfinite(float(_MRv_gate)):
9360 initial_MR_error = abs(float(_MRv_gate) - float(optimal_of)) / float(optimal_of)
9362 initial_MR_error = 1.0
9365 Cf_final = initial_performance.get(
"Cf_actual", initial_performance.get(
"Cf"))
9366 if Cf_final
is not None and np.isfinite(Cf_final):
9367 if optimized_config.chamber_geometry
is not None:
9368 optimized_config.chamber_geometry.Cf = float(Cf_final)
9371 sr_live = initial_performance.get(
"stability_results")
9372 snap = opt_state.get(
"best_results_for_validation")
or {}
9373 if isinstance(sr_live, dict)
and sr_live:
9374 stability_results = sr_live
9375 stability_state = stability_results.get(
"stability_state",
"unstable")
9376 stability_score = float(stability_results.get(
"stability_score", 0.0))
9377 chugging_margin = float(stability_results.get(
"chugging", {}).get(
"stability_margin", 0))
9378 acoustic_margin = float(stability_results.get(
"acoustic", {}).get(
"stability_margin", 0))
9379 feed_margin = float(stability_results.get(
"feed_system", {}).get(
"stability_margin", 0))
9381 stability_results = snap.get(
"stability_results", {})
9382 stability_state = snap.get(
"stability_state",
"unstable")
9383 stability_score = float(snap.get(
"stability_score", 0.0))
9384 chugging_margin = float(snap.get(
"chugging_margin", 0))
9385 acoustic_margin = float(snap.get(
"acoustic_margin", 0))
9386 feed_margin = float(snap.get(
"feed_margin", 0))
9388 stability_results = {}
9389 stability_state =
"unstable"
9390 stability_score = 0.0
9391 chugging_margin = acoustic_margin = feed_margin = 0.0
9392 initial_stability =
min(chugging_margin, acoustic_margin, feed_margin)
9396 stability_state, stability_score, chugging_margin, acoustic_margin, feed_margin,
9397 stability_results, requirements, min_stability)
9398 require_stable_state = requirements.get(
"require_stable_state",
True)
9400 _thr_v = requirements.get(
"layer1_thrust_validation_rel_tol")
9401 thrust_tol_valid = float(_thr_v)
if _thr_v
is not None else float(thrust_tol)
9402 _of_v = requirements.get(
"layer1_of_validation_tol")
9403 of_tol_valid = float(_of_v
if _of_v
is not None else 0.15)
9404 if not np.isfinite(of_tol_valid)
or of_tol_valid <= 0.0:
9406 _gate_tol_eps = 1e-12
9407 thrust_check_passed = initial_thrust_error <= thrust_tol_valid + _gate_tol_eps
9408 of_check_passed = initial_MR_error <= of_tol_valid + _gate_tol_eps
9411 A_throat_final = float(getattr(optimized_config.chamber_geometry,
'A_throat', 0.0))
9412 D_chamber_final = float(getattr(optimized_config.chamber_geometry,
'chamber_diameter', 0.0))
9413 A_chamber_final = np.pi * (D_chamber_final / 2) ** 2
9415 geometry_check_passed =
True
9416 geometry_failure_reasons = []
9418 if A_chamber_final > 0
and A_throat_final > 0:
9419 contraction_ratio_final = A_chamber_final / A_throat_final
9420 if contraction_ratio_final < 1.1:
9421 geometry_check_passed =
False
9422 geometry_failure_reasons.append(f
"Throat area too large for chamber (contraction ratio {contraction_ratio_final:.2f} < 1.1)")
9425 getattr(optimized_config.injector,
"type",
None) ==
"pintle"
9426 and hasattr(optimized_config.injector,
'geometry')
9427 and hasattr(optimized_config.injector.geometry,
'fuel')
9429 d_pintle_final = float(optimized_config.injector.geometry.fuel.d_pintle_tip)
9430 if d_pintle_final >= D_chamber_final:
9431 geometry_check_passed =
False
9432 geometry_failure_reasons.append(f
"Pintle diameter {d_pintle_final*1000:.1f}mm >= chamber diameter {D_chamber_final*1000:.1f}mm")
9435 getattr(optimized_config.injector,
"type",
None) ==
"impinging"
9436 and hasattr(optimized_config.injector,
"geometry")
9438 ig = optimized_config.injector.geometry
9439 d_max =
max(float(ig.oxidizer.d_jet), float(ig.fuel.d_jet))
9440 if D_chamber_final > 0
and d_max * 1.15 >= D_chamber_final:
9441 geometry_check_passed =
False
9442 geometry_failure_reasons.append(
9443 f
"Impinging jet diameter (max {d_max*1000:.2f} mm) too large vs chamber inner Ø {D_chamber_final*1000:.1f} mm"
9449 _cg_len = _ensure_cg_len(optimized_config)
9450 _len_raw = getattr(_cg_len,
"length",
None)
9453 if (_len_raw
is not None and np.isfinite(float(_len_raw)))
9456 _mcl = requirements.get(
"max_chamber_length_m")
9457 max_chamber_length = float(_mcl)
if _mcl
is not None else 0.50
9459 if np.isfinite(L_chamber_final)
and max_chamber_length > 0:
9460 if L_chamber_final > max_chamber_length:
9461 geometry_check_passed =
False
9462 geometry_failure_reasons.append(f
"Chamber length {L_chamber_final*1000:.1f}mm > max allowed {max_chamber_length*1000:.1f}mm")
9465 def _as_finite_float_or_nan(v: Any) -> float:
9468 except (TypeError, ValueError):
9470 return f
if np.isfinite(f)
else float(
"nan")
9472 ratio_o_val = ratio_f_val = float(
"nan")
9473 dp_gate_passed =
True
9474 pc_val = float(initial_performance.get(
"Pc", np.nan))
9475 if np.isfinite(pc_val)
and pc_val > 0:
9476 _ro, _rf = injector_dp_ratios_from_eval_result(pc_val, initial_performance)
9477 ratio_o_val = _as_finite_float_or_nan(_ro)
9478 ratio_f_val = _as_finite_float_or_nan(_rf)
9479 _chk_o_final = injector_dp_ratio_within_gate(
9480 ratio_o_val, layer1_dp_o_band[0], layer1_dp_o_band[1]
9482 if _chk_o_final
is not None:
9483 dp_gate_passed &= bool(_chk_o_final)
9484 _chk_f_final = injector_dp_ratio_within_gate(
9485 ratio_f_val, layer1_dp_f_band[0], layer1_dp_f_band[1]
9487 if _chk_f_final
is not None:
9488 dp_gate_passed &= bool(_chk_f_final)
9491 momentum_gate_passed =
True
9492 momentum_gate_band = (layer1_impinging_R_mom_lo, layer1_impinging_R_mom_hi)
9493 momentum_r_val = float(initial_performance.get(
"momentum_ratio_R", np.nan))
9495 getattr(optimized_config.injector,
"type",
None) ==
"impinging"
9496 and layer1_impinging_R_mom_lo
is not None
9497 and layer1_impinging_R_mom_hi
is not None
9498 and np.isfinite(momentum_r_val)
9499 and momentum_r_val > 0
9508 momentum_r_val, layer1_impinging_R_mom_lo, layer1_impinging_R_mom_hi,
9509 layer1_momentum_gate_safe_slack)
9512 resultant_tilt_gate_passed =
True
9513 _tilt_final = float(
"nan")
9514 if getattr(optimized_config.injector,
"type",
None) ==
"impinging":
9515 _gf = getattr(optimized_config.injector,
"geometry",
None)
9517 initial_performance.get(
"mdot_O"), initial_performance.get(
"mdot_F"),
9518 float(optimized_config.fluids[
"oxidizer"].density),
9519 float(optimized_config.fluids[
"fuel"].density),
9520 getattr(_gf.oxidizer,
"n_elements", float(
"nan"))
if _gf
else float(
"nan"),
9521 getattr(_gf.oxidizer,
"d_jet", float(
"nan"))
if _gf
else float(
"nan"),
9522 getattr(_gf.fuel,
"d_jet", float(
"nan"))
if _gf
else float(
"nan"),
9523 getattr(_gf.oxidizer,
"impingement_angle", float(
"nan"))
if _gf
else float(
"nan"),
9524 getattr(_gf.fuel,
"impingement_angle", float(
"nan"))
if _gf
else float(
"nan"),
9525 lox_inboard=layer1_ring_order_fuel_outboard,
9526 spacing_O_m=getattr(_gf.oxidizer,
"spacing", float(
"nan"))
if _gf
else float(
"nan"),
9527 spacing_F_m=getattr(_gf.fuel,
"spacing", float(
"nan"))
if _gf
else float(
"nan"),
9529 initial_performance[
"resultant_tilt_deg"] = (
9530 float(_tilt_final)
if np.isfinite(_tilt_final)
else None)
9535 _cg_f = getattr(optimized_config,
"chamber_geometry",
None)
9536 _L_cham_f = float(
"nan")
9537 if _cg_f
is not None:
9538 _lc = getattr(_cg_f,
"length_cylindrical",
None)
9539 _ln = getattr(_cg_f,
"length_contraction",
None)
9540 if _lc
is not None and _ln
is not None:
9541 _L_cham_f = float(_lc) + float(_ln)
9542 elif getattr(_cg_f,
"length",
None)
is not None:
9543 _L_cham_f = float(_cg_f.length)
9545 n_elements=getattr(_gf.oxidizer,
"n_elements", float(
"nan"))
if _gf
else float(
"nan"),
9546 spacing_O_m=getattr(_gf.oxidizer,
"spacing", float(
"nan"))
if _gf
else float(
"nan"),
9547 spacing_F_m=getattr(_gf.fuel,
"spacing", float(
"nan"))
if _gf
else float(
"nan"),
9548 angle_O_deg=getattr(_gf.oxidizer,
"impingement_angle", float(
"nan"))
if _gf
else float(
"nan"),
9549 angle_F_deg=getattr(_gf.fuel,
"impingement_angle", float(
"nan"))
if _gf
else float(
"nan"),
9550 D_chamber_inner_m=float(getattr(_cg_f,
"chamber_diameter", float(
"nan")))
9551 if _cg_f
is not None else float(
"nan"),
9552 L_chamber_m=_L_cham_f,
9554 initial_performance[
"resultant_tilt_breakeven_deg"] = (
9555 float(_be_final)
if np.isfinite(_be_final)
else None)
9557 from_reach=layer1_resultant_tilt_from_reach,
9558 constant_deg=layer1_resultant_tilt_max_deg,
9559 breakeven_deg=_be_final,
9560 margin=layer1_resultant_tilt_reach_margin,
9562 initial_performance[
"resultant_tilt_allowed_deg"] = float(_tilt_allow_final)
9563 if np.isfinite(_be_final)
and float(_tilt_allow_final) > _be_final:
9564 layer1_logger.warning(
9565 "The permitted outward tilt, %.2f deg, exceeds the break-even tilt for this "
9566 "design's own geometry (%.2f deg): a fan at the allowed limit would reach the "
9567 "liner %.0f mm from the face, inside a %.0f mm chamber. Lower the allowance, or "
9568 "move the element ring inboard.",
9569 float(_tilt_allow_final), _be_final,
9570 1000.0 * _L_cham_f * float(np.tan(np.deg2rad(_be_final)))
9571 /
max(1e-9, float(np.tan(np.deg2rad(float(_tilt_allow_final))))),
9574 if np.isfinite(_tilt_final):
9575 resultant_tilt_gate_passed = bool(
9576 _tilt_final <= float(_tilt_allow_final)
9577 + float(layer1_resultant_tilt_gate_tol_deg))
9578 initial_performance[
"resultant_tilt_gate_passed"] = bool(resultant_tilt_gate_passed)
9584 performance=initial_performance)
9585 for _gname, (_gok, _gmsg)
in limit_gates.items():
9586 initial_performance[f
"{_gname}_gate_passed"] = bool(_gok)
9588 geometry_check_passed =
False
9589 geometry_failure_reasons.append(_gmsg)
9594 _sc_v = float(initial_performance.get(
"layer1_validation_tank_pressure_scale", 1.0)
or 1.0)
9595 validation_replay_ok =
True
9596 if abs(_sc_v - 1.0) > 1e-9:
9597 validation_replay_ok =
False
9598 failure_reasons_replay = (
9599 f
"Design does not close at its own tank pressures (needs x{_sc_v:.2f}: "
9600 f
"{float(P_O_initial) / psi_to_Pa * _sc_v:.0f}/{float(P_F_initial) / psi_to_Pa * _sc_v:.0f} psi)")
9601 elif initial_performance.get(
"layer1_validation_used_last_good_bundle"):
9602 validation_replay_ok =
False
9603 failure_reasons_replay = (
9604 "Validation replay failed ("
9605 + str(initial_performance.get(
"layer1_validation_replay_error",
"evaluate() raised"))
9606 +
"); the numbers shown are an earlier in-loop evaluation")
9607 elif initial_performance.get(
"layer1_validation_replay_error"):
9608 validation_replay_ok =
False
9609 failure_reasons_replay = (
"Validation replay failed: "
9610 + str(initial_performance[
"layer1_validation_replay_error"]))
9612 failure_reasons_replay =
""
9613 initial_performance[
"validation_replay_ok"] = bool(validation_replay_ok)
9616 pc_gate_passed =
True
9618 _pc_t_v = requirements.get(
"target_chamber_pressure_psi")
9619 _pc_v = float(initial_performance.get(
"Pc", np.nan))
9620 if _pc_t_v
is not None and float(_pc_t_v) > 0:
9621 _pc_miss = abs(_pc_v / (float(_pc_t_v) * 6894.757293168) - 1.0)
if np.isfinite(_pc_v)
else np.inf
9622 pc_gate_passed = bool(_pc_miss <= _LAYER1_PC_GATE_REL)
9623 if not pc_gate_passed:
9625 f
"Chamber pressure {_pc_v / 6894.757293168:.1f} psia is {100 * _pc_miss:.2f} % off "
9626 f
"the {float(_pc_t_v):.0f} psia target (limit {100 * _LAYER1_PC_GATE_REL:.0f} %)")
9627 initial_performance[
"pc_gate_passed"] = bool(pc_gate_passed)
9629 pressure_candidate_valid = (
9633 and resultant_tilt_gate_passed
9634 and stability_check_passed
9635 and geometry_check_passed
9637 and momentum_gate_passed
9638 and validation_replay_ok
9642 opt_state[
"converged"] = bool(pressure_candidate_valid)
9645 failure_reasons = []
9646 if not validation_replay_ok:
9647 failure_reasons.append(failure_reasons_replay)
9648 if not pc_gate_passed:
9649 failure_reasons.append(pc_gate_reason)
9650 if not thrust_check_passed:
9651 failure_reasons.append(f
"Thrust error {initial_thrust_error*100:.1f}% > {thrust_tol_valid*100:.1f}% limit")
9652 if not of_check_passed:
9653 failure_reasons.append(f
"O/F error {initial_MR_error*100:.1f}% > {of_tol_valid*100:.1f}% limit")
9654 if not geometry_check_passed:
9655 failure_reasons.extend(geometry_failure_reasons)
9656 if not stability_check_passed:
9657 failure_reasons.append(f
"Stability failed: {'; '.join(stability_parts)}")
9658 if not dp_gate_passed:
9660 if np.isfinite(ratio_o_val):
9662 f
"ΔP_O/Pc in [{layer1_dp_o_band[0]:.3f},{layer1_dp_o_band[1]:.3f}] (got {ratio_o_val:.3f})"
9664 if np.isfinite(ratio_f_val):
9666 f
"ΔP_F/Pc in [{layer1_dp_f_band[0]:.3f},{layer1_dp_f_band[1]:.3f}] (got {ratio_f_val:.3f})"
9668 failure_reasons.append(
"Injector ΔP ratio failed: " +
"; ".join(dp_parts))
9669 if not momentum_gate_passed:
9670 failure_reasons.append(
9671 f
"Momentum ratio R not in [{float(momentum_gate_band[0]):.3f},"
9672 f
"{float(momentum_gate_band[1]):.3f}] (got {momentum_r_val:.3f})"
9675 if not pressure_candidate_valid
and not failure_reasons:
9676 failure_reasons.append(
"Validation failed: no requirements met")
9679 if pressure_candidate_valid:
9680 update_progress(
"Layer 1: Validation", 0.53, f
"✓ VALID - Thrust err: {initial_thrust_error*100:.1f}%, O/F err: {initial_MR_error*100:.1f}%, Stability: {stability_state}")
9681 log_status(
"Layer 1", f
"VALID | Thrust err {initial_thrust_error*100:.1f}%, O/F err {initial_MR_error*100:.1f}%, Stability {stability_state}")
9682 layer1_logger.info(
"✓ Validation: VALID")
9684 update_progress(
"Layer 1: Validation", 0.53, f
"✗ INVALID - {'; '.join(failure_reasons)}")
9685 log_status(
"Layer 1", f
"INVALID | Reasons: {', '.join(failure_reasons)}")
9686 layer1_logger.warning(f
"✗ Validation: INVALID - {'; '.join(failure_reasons)}")
9689 final_performance = initial_performance.copy()
9690 final_performance[
"pressure_candidate_valid"] = pressure_candidate_valid
9691 final_performance[
"initial_thrust_error"] = initial_thrust_error
9692 final_performance[
"initial_MR_error"] = initial_MR_error
9693 final_performance[
"initial_stability"] = initial_stability
9694 final_performance[
"initial_stability_state"] = stability_state
9695 final_performance[
"initial_stability_score"] = stability_score
9696 final_performance[
"thrust_check_passed"] = thrust_check_passed
9697 final_performance[
"of_check_passed"] = of_check_passed
9698 final_performance[
"stability_check_passed"] = stability_check_passed
9699 final_performance[
"geometry_check_passed"] = geometry_check_passed
9700 final_performance[
"dp_gate_passed"] = dp_gate_passed
9701 final_performance[
"momentum_gate_passed"] = momentum_gate_passed
9702 final_performance[
"injector_dp_ratio_O"] = ratio_o_val
if np.isfinite(ratio_o_val)
else None
9703 final_performance[
"injector_dp_ratio_F"] = ratio_f_val
if np.isfinite(ratio_f_val)
else None
9704 final_performance[
"layer1_Lstar_search_min_m"] = float(min_Lstar)
9705 final_performance[
"layer1_Lstar_search_max_m"] = float(max_Lstar)
9706 final_performance[
"layer1_W_geom_ao_af_momentum_effective"] = float(layer1_W_geom_ao_af)
9708 if isinstance(requirements, dict):
9710 str(k): requirements.get(k)
for k
in sorted(requirements.keys(), key=str)
9712 req_used[
"W_geom_ao_af_momentum"] = float(layer1_W_geom_ao_af)
9713 final_performance[
"layer1_requirements_used"] = req_used
9715 final_performance[
"layer1_requirements_used"] = {}
9717 final_performance[
"layer1_requirements_used"] = {}
9720 _cgd = _eg_lstar_diag(optimized_config)
9721 _atd = float(_cgd.A_throat)
if _cgd.A_throat
else float(
"nan")
9722 _vold = float(_cgd.volume)
if _cgd.volume
is not None else float(
"nan")
9723 _lscd = float(_cgd.Lstar)
if _cgd.Lstar
is not None else float(
"nan")
9724 _dind = float(_cgd.chamber_diameter)
if _cgd.chamber_diameter
else float(
"nan")
9725 _acd = float(np.pi * (_dind / 2.0) ** 2)
if np.isfinite(_dind)
and _dind > 0
else float(
"nan")
9726 final_performance[
"layer1_geometry_Lstar_config_m"] = _lscd
if np.isfinite(_lscd)
else None
9727 final_performance[
"layer1_geometry_volume_over_At_m"] = (
9728 float(_vold / _atd)
if np.isfinite(_vold)
and np.isfinite(_atd)
and _atd > 0
else None
9730 final_performance[
"layer1_geometry_Ac_over_At"] = (
9731 float(_acd / _atd)
if np.isfinite(_acd)
and np.isfinite(_atd)
and _atd > 0
else None
9733 if np.isfinite(_vold)
and np.isfinite(_atd)
and _atd > 0
and np.isfinite(_lscd):
9734 final_performance[
"layer1_geometry_V_equals_Lstar_At"] = bool(
9735 np.isclose(_vold, _lscd * _atd, rtol=1e-4, atol=1e-8)
9739 if getattr(optimized_config.injector,
"type",
None) ==
"impinging":
9740 _diag_fp = final_performance.get(
"diagnostics", {})
if isinstance(final_performance, dict)
else {}
9741 _diag_map = _diag_fp
if isinstance(_diag_fp, dict)
else {}
9743 _mr_fp_raw = float(final_performance.get(
"MR", np.nan))
9744 except (TypeError, ValueError):
9745 _mr_fp_raw = float(
"nan")
9746 _mr_fp = _mr_fp_raw
if np.isfinite(_mr_fp_raw)
and _mr_fp_raw > 0
else None
9749 target_smd_microns=layer1_target_smd_microns,
9750 smd_rel_tol=layer1_smd_rel_tol,
9754 if layer1_W_SMD <= 0.0:
9757 _v_o = _diag_map.get(
"D32_O")
9758 if _v_o
is not None:
9761 except (TypeError, ValueError):
9762 _fv_o = float(
"nan")
9763 if np.isfinite(_fv_o)
and _fv_o > 0:
9765 _v_f = _diag_map.get(
"D32_F")
9766 if _v_f
is not None:
9769 except (TypeError, ValueError):
9770 _fv_f = float(
"nan")
9771 if np.isfinite(_fv_f)
and _fv_f > 0:
9773 if np.isfinite(_do)
and _do > 0
and np.isfinite(_df)
and _df > 0:
9774 _mr_b = _mr_fp
if _mr_fp
is not None else 3.0
9775 _w_o = float(_mr_b / (1.0 + _mr_b))
9776 _w_f = float(1.0 / (1.0 + _mr_b))
9777 _smd_eff_um_fp = float((_w_o * _do + _w_f * _df) * 1e6)
9778 elif np.isfinite(_do)
and _do > 0:
9779 _smd_eff_um_fp = float(_do * 1e6)
9780 elif np.isfinite(_df)
and _df > 0:
9781 _smd_eff_um_fp = float(_df * 1e6)
9782 final_performance[
"effective_smd_microns"] = _smd_eff_um_fp
if np.isfinite(_smd_eff_um_fp)
else None
9783 final_performance[
"impingement_angle_deg_effective"] = _imp_deg_fp
if np.isfinite(_imp_deg_fp)
else None
9784 final_performance[
"smd_penalty"] = float(layer1_W_SMD * _smd_term_fp)
9785 for _gk
in (
"L_imp",
"D_pitch_O",
"D_pitch_F",
"element_gap_O",
"element_gap_F",
9786 "vaporization_length_total"):
9787 _gv = _diag_map.get(_gk)
9789 _gvf = float(_gv)
if _gv
is not None else float(
"nan")
9790 except (TypeError, ValueError):
9792 final_performance[_gk] = _gvf
if np.isfinite(_gvf)
else None
9793 final_performance[
"failure_reasons"] = failure_reasons
9794 final_performance[
"layer1_warnings"] = list(layer1_warnings)
9796 final_performance[
"chugging_margin"] = chugging_margin
9797 final_performance[
"acoustic_margin"] = acoustic_margin
9798 final_performance[
"feed_margin"] = feed_margin
9802 P_O_start_optimized_psi = float(optimized_config.lox_tank.initial_pressure_psi)
9803 P_F_start_optimized_psi = float(optimized_config.fuel_tank.initial_pressure_psi)
9804 final_performance[
"P_O_start_psi"] = P_O_start_optimized_psi
9805 final_performance[
"P_F_start_psi"] = P_F_start_optimized_psi
9806 final_performance[
"P_O_start_ratio"] = P_O_start_optimized_psi / max_lox_P_psi
if max_lox_P_psi > 0
else 0.0
9807 final_performance[
"P_F_start_ratio"] = P_F_start_optimized_psi / max_fuel_P_psi
if max_fuel_P_psi > 0
else 0.0
9809 if best_x
is not None and len(best_x) > idx_P_F:
9810 P_O_start_optimized_psi = float(np.clip(best_x[idx_P_O], bounds[idx_P_O][0], bounds[idx_P_O][1]))
9811 P_F_start_optimized_psi = float(np.clip(best_x[idx_P_F], bounds[idx_P_F][0], bounds[idx_P_F][1]))
9812 final_performance[
"P_O_start_psi"] = P_O_start_optimized_psi
9813 final_performance[
"P_F_start_psi"] = P_F_start_optimized_psi
9814 final_performance[
"P_O_start_ratio"] = P_O_start_optimized_psi / max_lox_P_psi
if max_lox_P_psi > 0
else 0.0
9815 final_performance[
"P_F_start_ratio"] = P_F_start_optimized_psi / max_fuel_P_psi
if max_fuel_P_psi > 0
else 0.0
9817 final_performance[
"P_O_start_psi"] = max_lox_P_psi * 0.8
9818 final_performance[
"P_F_start_psi"] = max_fuel_P_psi * 0.8
9819 final_performance[
"P_O_start_ratio"] = 0.8
9820 final_performance[
"P_F_start_ratio"] = 0.8
9823 if "chamber_intrinsics" in initial_performance:
9824 final_performance[
"chamber_intrinsics"] = initial_performance[
"chamber_intrinsics"]
9831 cg = ensure_chamber_geometry(optimized_config)
9832 A_throat_final = float(cg.A_throat)
if cg.A_throat
and cg.A_throat > 0
else None
9834 inj = getattr(optimized_config,
'injector',
None)
9835 inj_geom = getattr(inj,
'geometry',
None)
if inj
else None
9836 is_pintle = getattr(inj,
'type',
None) ==
'pintle' if inj
else False
9837 is_impinging = getattr(inj,
'type',
None) ==
'impinging' if inj
else False
9839 if is_impinging
and inj_geom
and A_throat_final:
9840 ox = getattr(inj_geom,
'oxidizer',
None)
9841 fu = getattr(inj_geom,
'fuel',
None)
9843 n_o = float(getattr(ox,
'n_elements', 0))
9844 d_o = float(getattr(ox,
'd_jet', 0))
9845 n_f = float(getattr(fu,
'n_elements', 0))
9846 d_f = float(getattr(fu,
'd_jet', 0))
9847 A_lox_g = n_o * np.pi * (d_o / 2.0) ** 2
if d_o > 0
else 0.0
9848 A_fuel_g = n_f * np.pi * (d_f / 2.0) ** 2
if d_f > 0
else 0.0
9849 diag_ev = initial_performance.get(
"diagnostics")
or {}
9850 cd_o = diag_ev.get(
"Cd_O")
9851 cd_f = diag_ev.get(
"Cd_F")
9853 cd_o
is not None and np.isfinite(cd_o)
and float(cd_o) > 0
9854 and cd_f
is not None and np.isfinite(cd_f)
and float(cd_f) > 0
9856 A_eff = float(cd_o) * A_lox_g + float(cd_f) * A_fuel_g
9861 oxg = getattr(inj_geom,
"oxidizer",
None)
9862 fug = getattr(inj_geom,
"fuel",
None)
9863 d_o = float(getattr(oxg,
"d_jet", 0.002)
or 0.002)
9864 d_f = float(getattr(fug,
"d_jet", 0.002)
or 0.002)
9865 Cd_lox = cd_inf_from_orifice_diameter(
9866 d_o, optimized_config.discharge[
"oxidizer"]
9868 Cd_fu = cd_inf_from_orifice_diameter(
9869 d_f, optimized_config.discharge[
"fuel"]
9872 Cd_lox, Cd_fu = 0.6, 0.6
9873 A_eff = Cd_lox * A_lox_g + Cd_fu * A_fuel_g
9874 final_performance[
"effective_injector_area_ratio"] = float(A_eff / A_throat_final)
9876 if is_pintle
and inj_geom
and A_throat_final:
9878 lox_geom = getattr(inj_geom,
'lox',
None)
9879 fuel_geom = getattr(inj_geom,
'fuel',
None)
9881 if lox_geom
and fuel_geom:
9882 n_orifices = float(getattr(lox_geom,
'n_orifices', 0))
9883 d_orifice = float(getattr(lox_geom,
'd_orifice', 0))
9884 d_pintle_tip = float(getattr(fuel_geom,
'd_pintle_tip', 0))
9885 h_gap = float(getattr(fuel_geom,
'h_gap', 0))
9886 Cd_lox = float(getattr(lox_geom,
'Cd', 0.65))
9887 Cd_fuel = float(getattr(fuel_geom,
'Cd', 0.4))
9889 A_lox_injector = n_orifices * np.pi * (d_orifice / 2) ** 2
if d_orifice > 0
else 0.0
9890 R_inner = d_pintle_tip / 2
9891 R_outer = R_inner + h_gap
9892 A_fuel_injector = np.pi * (R_outer ** 2 - R_inner ** 2)
if h_gap > 0
else 0.0
9895 A_lox_effective = Cd_lox * A_lox_injector
9896 A_fuel_effective = Cd_fuel * A_fuel_injector
9899 A_injector_total_effective = A_lox_effective + A_fuel_effective
9902 effective_injector_area_ratio = A_injector_total_effective / A_throat_final
9903 final_performance[
"effective_injector_area_ratio"] = float(effective_injector_area_ratio)
9904 except Exception
as e:
9905 layer1_logger.debug(f
"Could not calculate effective_injector_area_ratio: {e}")
9909 target_thrust=target_thrust,
9910 optimal_of=optimal_of,
9911 target_P_exit=target_P_exit,
9917 "injector_type": getattr(getattr(optimized_config,
"injector",
None),
"type",
None),
9918 "performance": final_performance,
9919 "iteration_history": iteration_history,
9920 "convergence_info": {
9921 "converged": opt_state[
"converged"],
9922 "iterations": len(iteration_history),
9923 "final_change": opt_state[
"best_objective"],
9924 "best_objective": float(opt_state.get(
"best_objective", float(
"inf"))),
9925 "best_objective_breakdown": (
9926 dict(opt_state.get(
"best_objective_breakdown", {}))
9927 if isinstance(opt_state.get(
"best_objective_breakdown", {}), dict)
9930 "primary_relative_residual": _prim_rel,
9931 "infeasible_reason": infeasible_reason,
9933 "exit_pressure_targeting": {
9934 "target_P_exit": target_P_exit,
9936 "optimized_pressure_curves": {
9937 "lox_end_ratio": final_lox_end_ratio,
9938 "fuel_end_ratio": final_fuel_end_ratio,
9939 "lox_start_psi": final_performance[
"P_O_start_psi"],
9940 "fuel_start_psi": final_performance[
"P_F_start_psi"],
9943 "layer_1_pressure_candidate": pressure_candidate_valid,
9945 "optimized_parameters": extract_all_parameters(optimized_config),
9946 "final_rebuild_diag": opt_state.get(
"final_rebuild_diag"),
9950 layer1_logger.info(
"")
9951 layer1_logger.info(
"="*70)
9952 layer1_logger.info(
"Final Results Summary")
9953 layer1_logger.info(
"="*70)
9954 if "best_results_for_validation" in opt_state
and opt_state[
"best_results_for_validation"]
is not None:
9955 stored = opt_state[
"best_results_for_validation"]
9956 layer1_logger.info(f
"Thrust: {stored.get('F', 0):.1f} N (target: {target_thrust:.1f} N)")
9957 layer1_logger.info(f
"Thrust error: {stored.get('thrust_error', 0)*100:.2f}%")
9958 layer1_logger.info(f
"O/F ratio: {stored.get('MR', 0):.3f} (target: {optimal_of:.3f})")
9959 layer1_logger.info(f
"O/F error: {stored.get('of_error', 0)*100:.2f}%")
9960 layer1_logger.info(f
"Stability state: {stored.get('stability_state', 'unknown')}")
9961 layer1_logger.info(f
"Stability score: {stored.get('stability_score', 0):.3f}")
9962 layer1_logger.info(f
"Chugging margin: {stored.get('chugging_margin', 0):.3f}")
9963 layer1_logger.info(f
"Acoustic margin: {stored.get('acoustic_margin', 0):.3f}")
9964 layer1_logger.info(f
"Feed margin: {stored.get('feed_margin', 0):.3f}")
9965 if final_performance.get(
"P_O_start_psi")
is not None:
9966 layer1_logger.info(f
"LOX initial pressure: {final_performance['P_O_start_psi']:.1f} psi")
9967 if final_performance.get(
"P_F_start_psi")
is not None:
9968 layer1_logger.info(f
"Fuel initial pressure: {final_performance['P_F_start_psi']:.1f} psi")
9969 _idr_o = final_performance.get(
"injector_dp_ratio_O")
9970 _idr_f = final_performance.get(
"injector_dp_ratio_F")
9971 _idp = final_performance.get(
"injector_dp_penalty")
9972 if _idr_o
is not None and np.isfinite(_idr_o):
9973 layer1_logger.info(f
"Injector ΔP/Pc (LOX): {_idr_o:.4f}")
9974 if _idr_f
is not None and np.isfinite(_idr_f):
9975 layer1_logger.info(f
"Injector ΔP/Pc (fuel): {_idr_f:.4f}")
9976 if _idp
is not None and np.isfinite(_idp):
9977 layer1_logger.info(f
"Injector ΔP ratio penalty (weighted): {float(_idp):.6g}")
9979 getattr(getattr(optimized_config,
"injector",
None),
"type",
None) ==
"impinging"
9980 and final_performance.get(
"momentum_ratio_R")
is not None
9983 f
"Impinging jet momentum ratio R: {final_performance['momentum_ratio_R']:.3f} "
9984 f
"(sqrt momentum-flux ratio; target ~1.0)"
9986 mbp = final_performance.get(
"momentum_balance_penalty")
9987 if mbp
is not None and np.isfinite(mbp):
9989 f
"Impinging momentum-balance penalty (W_MOM × hinge): {float(mbp):.6g}"
9991 layer1_logger.info(f
"Validation: {'VALID' if pressure_candidate_valid else 'INVALID'}")
9996 target_thrust=target_thrust,
9997 optimal_of=optimal_of,
9998 target_P_exit=target_P_exit,
9999 layer1_exit_pressure_inside_quad_scale=layer1_exit_pressure_inside_quad,
10000 layer1_w_dp=layer1_W_DP,
10001 layer1_w_dp_o=layer1_W_DP_O,
10002 layer1_w_dp_f=layer1_W_DP_F,
10003 injector_dp_o_band=layer1_dp_o_band,
10004 injector_dp_f_band=layer1_dp_f_band,
10005 injector_dp_o_soft_floor=layer1_dp_o_soft_floor,
10006 layer1_w_dp_o_floor=layer1_W_DP_O_FLOOR,
10007 layer1_w_dp_center=layer1_W_DP_CENTER,
10008 layer1_w_thrust=layer1_W_THRUST_obj,
10009 layer1_w_of=layer1_W_OF_obj,
10010 layer1_w_of_low_mr_scale=layer1_W_OF_low_MR_scale,
10011 layer1_w_of_high_mr_scale=layer1_W_OF_high_MR_scale,
10013 layer1_logger.info(
"")
10014 layer1_logger.info(
10015 "Final evaluate — primary objective penalty contributions "
10016 "(W_THRUST/W_OF/W_EXIT × terms + injector_dp weighted; same structure as scalar objective)"
10018 layer1_logger.info(
10019 f
" • Thrust penalty contribution: {_primary_terms['thrust_penalty_contribution']:.6f}"
10021 layer1_logger.info(
10022 f
" • O/F penalty contribution: {_primary_terms['of_penalty_contribution']:.6f}"
10024 layer1_logger.info(
10025 f
" • Exit pressure penalty contribution: {_primary_terms['exit_pressure_penalty_contribution']:.6f}"
10027 layer1_logger.info(
10028 f
" • Injector ΔP penalty contribution: {_primary_terms['injector_dp_penalty_contribution']:.6f}"
10030 layer1_logger.info(
10031 " • Primary relative |ΔF|/F, |ΔMR|/MR, |ΔP_exit|/P_tgt: "
10032 f
"{_prim_rel['rel_thrust']:.4g}, {_prim_rel['rel_of']:.4g}, {_prim_rel['rel_P_exit']:.4g} "
10033 f
"(RMS={_prim_rel['rms_primary']:.4g}, dimensionless; not the weighted scalar)"
10035 if final_performance.get(
"injector_dp_out_of_range"):
10036 layer1_logger.warning(
10037 " injector_dp_out_of_range: True (one or both ΔP_inj/Pc streams exceed 0.50)"
10041 best_obj_value = opt_state.get(
"best_objective", float(
'inf'))
10042 if best_obj_value > 1.0:
10043 layer1_logger.info(
"")
10044 layer1_logger.info(
"-"*70)
10045 layer1_logger.info(f
"⚠ Objective Function Breakdown (objective = {best_obj_value:.6f} > 1.0)")
10046 layer1_logger.info(
"-"*70)
10047 breakdown = opt_state.get(
"best_objective_breakdown", {})
10049 if breakdown.get(
"is_infeasible",
False):
10050 layer1_logger.info(
"Solution is INFEASIBLE:")
10051 if breakdown.get(
"length_violation",
False):
10052 layer1_logger.info(f
" • Chamber length violation penalty: {breakdown.get('infeasibility_penalty', 0):.6f}")
10054 layer1_logger.info(f
" • Infeasibility penalty (constraints violated): {breakdown.get('infeasibility_penalty', 0):.6f}")
10056 layer1_logger.info(
"Objective component contributions:")
10058 (
"Thrust penalty", breakdown.get(
"thrust_penalty", 0)),
10059 (
"O/F ratio penalty", breakdown.get(
"of_penalty", 0)),
10060 (
"Cf band penalty", breakdown.get(
"cf_penalty", 0)),
10061 (
"Exit pressure penalty", breakdown.get(
"exit_pressure_penalty", 0)),
10062 (
"Injector ΔP penalty", breakdown.get(
"injector_dp_penalty", 0)),
10063 (
"Chamber length penalty", breakdown.get(
"length_penalty", 0)),
10064 (
"Chamber shape penalty", breakdown.get(
"chamber_shape_penalty", 0)),
10065 (
"Momentum balance penalty", breakdown.get(
"momentum_balance_penalty", 0)),
10066 (
"Impingement angle penalty", breakdown.get(
"impingement_angle_penalty", 0)),
10067 (
"Jet asymmetry penalty", breakdown.get(
"jet_angle_asymmetry_penalty", 0)),
10068 (
"SMD penalty", breakdown.get(
"smd_penalty", 0)),
10069 (
"Geometry fit penalty", breakdown.get(
"impinging_geometry_fit_penalty", 0)),
10070 (
"AO/AF momentum hint", breakdown.get(
"geom_ao_af_momentum_penalty", 0)),
10071 (
"Tank equal penalty", breakdown.get(
"tank_pressure_equal_penalty", 0)),
10074 components.sort(key=
lambda x: x[1], reverse=
True)
10075 total = sum(c[1]
for c
in components)
10076 for name, value
in components:
10077 if value > 0
or total > 0:
10078 pct = (value / total * 100)
if total > 0
else 0.0
10079 layer1_logger.info(f
" • {name}: {value:.6f} ({pct:.1f}%)")
10080 layer1_logger.info(f
" Total: {total:.6f}")
10082 layer1_logger.info(
" (component breakdown not available)")
10083 layer1_logger.info(
"-"*70)
10085 layer1_logger.info(
"="*70)
10086 layer1_logger.info(
"")
10087 layer1_logger.info(f
"Layer 1 optimization complete. Log saved to: {log_file_path}")
10090 layer1_logger.handlers.clear()
10092 update_progress(
"Layer 1: Complete", 1.0,
"Layer 1 optimization complete!")
10099 return optimized_config, results