Canonical RHS excerpt from the registered callable used for this benchmark cell. Expand it to verify the state equations; it is not a standalone runnable fixture.
Show canonical RHS excerpt
def _spot_efficiency(wander_x: float, wander_y: float,
slew_x: float, slew_y: float,
blooming: float) -> float:
"""Fraction of beam power deposited on target, accounting for
wander, slew jitter, and thermal blooming defocus."""
r2 = (wander_x + slew_x) ** 2 + (wander_y + slew_y) ** 2
gauss_loss = np.exp(-2.0 * r2 / (_W0 ** 2))
bloom_loss = 1.0 / (1.0 + blooming ** 2)
return gauss_loss * bloom_loss
def rhs_hel_thermal(t: float, y: np.ndarray) -> np.ndarray:
T_s = max(y[0], 1.0)
T_1 = max(y[1], 1.0)
T_2 = max(y[2], 1.0)
T_3 = max(y[3], 1.0)
T_melt = y[4]
bw_x, bw_y = y[6], y[7]
bloom = y[8]
slw_x, slw_y = y[9], y[10]
eff = _spot_efficiency(bw_x, bw_y, slw_x, slw_y, bloom)
P_abs = _P_BEAM * np.exp(-_ALPHA_ATM * _RANGE) * eff
rho_cp_dx = _RHO * _CP * _DX
k_dx2 = _K_COND / (_DX ** 2)
k_dx = _K_COND / _DX
d = np.empty(11)
# Thermal layers (1D finite difference, 4 layers)
irradiance = P_abs / (np.pi * _W0 ** 2)
d[0] = (irradiance
- _SIGMA * _EPS * (T_s ** 4 - _T_AMB ** 4)
- k_dx * (T_s - T_1)) / rho_cp_dx
d[1] = k_dx2 * (T_s - 2.0 * T_1 + T_2) / (_RHO * _CP)
d[2] = k_dx2 * (T_1 - 2.0 * T_2 + T_3) / (_RHO * _CP)
d[3] = k_dx2 * (T_2 - 2.0 * T_3 + _T_AMB) / (_RHO * _CP)
# Melt tracking
d[4] = 0.0 if T_s < _T_MELT_THRESH else (T_s - _T_MELT_THRESH) / _TAU_MELT
# Arrhenius damage integral
d[5] = _A_ARR * np.exp(-_E_A / (_R_GAS * T_s))
# OU surrogates for beam wander
d[6] = -bw_x / _TAU_TURB + _SIGMA_TURB * np.sin(_FREQ_WX * t)
d[7] = -bw_y / _TAU_TURB + _SIGMA_TURB * np.sin(_FREQ_WY * t)
# Thermal blooming index
d[8] = (_N2 * _P_BEAM / (np.pi * _W0 ** 2) - bloom) / _TAU_BLOOM
# OU surrogates for slew jitter
d[9] = -slw_x / _TAU_SLEW + _SIGMA_SLEW * np.sin(_FREQ_SX * t)
d[10] = -slw_y / _TAU_SLEW + _SIGMA_SLEW * np.sin(_FREQ_SY * t)
return d- Parameters
- _ALPHA_ATM = 0.0005
- _A_ARR = 1e+13
- _CP = 500
- _DX = 0.001
- _EPS = 0.9
- _E_A = 300000
- _FREQ_SX = 19.3
- _FREQ_SY = 29.7
- _FREQ_WX = 37.7
- _FREQ_WY = 53.1
- _K_COND = 50
- _N2 = 2e-19
- _P_BEAM = 30000
- _RANGE = 1000
- _RHO = 7800
- _R_GAS = 8.314
- _SIGMA = 5.67e-08
- _SIGMA_SLEW = 0.001
- _SIGMA_TURB = 0.002
- _TAU_BLOOM = 0.1
- _TAU_MELT = 1
- _TAU_SLEW = 0.05
- _TAU_TURB = 0.01
- _T_AMB = 300
- _T_MELT_THRESH = 1800
- _W0 = 0.05
- Initial condition
- y(0) = [300, 300, 300, 300, 0, 0, 0, 0, 0, 0, 0]
- Horizon
- t ∈ [0, 5]
Canonical RHS excerpt captured from the same registered callable used for the published benchmark. Frozen closure values are summarized below; helper imports and solver settings are intentionally omitted.