Физика v3: коэффициент заполнения железа + фикс зависаний решателя

ЭКСПЛОЙТ (нашёлся по подозрительным 75-82% КПД на дашборде): модель
умножала ВСЮ индуктивность катушки на (mu_eff-1)*overlap, как будто железо
заполняет всё сечение. Тонкий снаряд (4мм) в толстой катушке (29мм) получал
~50x нефизичной силы бесплатно, и эволюция сгрузилась именно туда.
Фикс: iron_fill_factor = A_снаряда/A_среднего витка (классическое
λ = n*(B_fe*A_fe + B_возд*(A-A_fe))).

ЗАВИСАНИЯ (сервер стоял 35+ мин на одной конфигурации):
- сглаживание сухого трения tanh(v/0.01) вместо разрывного sign(v)
- событие тока удержания ключа (SCR закрывается при I<0.1А после пика) --
  передемпфированный хвост больше не молотит весь 20мс бюджет
- полная формула отражения вихревого контура (ωM)²R2/(R2²+(ωL2)²) вместо
  резистивного приближения (114кОм -> 2.5кОм на вырожденных катушках)
- гард: L/R < 100нс = вне области применимости модели, мгновенный отказ
Бенч 200 геномов: 608с -> 22с, худшая оценка 540с -> 2.5с.

MODEL_VERSION -> gausse-physics-v3. 104 теста.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

src/gausse/gpu/batch_integrator.py

@@ -108,7 +108,8 @@ def _derivatives(xp, q, i, x, v, p: BatchDischargeParams):
     v_c = q / p.capacitance_f
     d_q = -i
     d_i = (v_c - i * r_eff - dlambda_dx * v) / dlambda_di
-    retard = xp.sign(v) * (p.retard_const_n + p.drag_coeff_n * v * v)  # трение + воздух
+    # трение + воздух; tanh-сглаживание у v=0 — как в physics/circuit.py
+    retard = xp.tanh(v / 0.01) * (p.retard_const_n + p.drag_coeff_n * v * v)
     d_v = (force - retard) / p.mass_kg
     d_x = v
     return d_q, d_i, d_x, d_v
@@ -145,7 +146,7 @@ def _get_fused_step(cp):
         reff = rt + red * ov
         dq = -i
         di = (q / C - i * reff - liron * dov * g * v) / dl_di
-        dv = (liron * dov * G - cp.sign(v) * (fc + dc * v * v)) / m
+        dv = (liron * dov * G - cp.tanh(v * 100.0) * (fc + dc * v * v)) / m
         return dq, di, v, dv
 
     @cp.fuse()
@@ -242,22 +243,25 @@ def integrate_batch_discharge(
         peak_current = xp.maximum(peak_current, xp.where(active, xp.abs(i_new), 0.0))
         past_peak = past_peak | (active & (i_old > 1.0) & (i_new < i_old))
 
-        # ОДИН импульс: обрыв на нуле тока ИЛИ на первом локальном минимуме
+        # ОДИН импульс: обрыв на нуле тока, ИЛИ на первом локальном минимуме
-        # (снаряд начал подкачивать ток обратно). Что раньше.
+        # (снаряд начал подкачивать ток обратно), ИЛИ при падении ниже тока
+        # удержания ключа (передемпфированный хвост; см. physics/circuit.py).
         crossed = active & (i_old > 0) & (i_new <= 0)
-        local_min = active & past_peak & (i_new > i_old) & ~crossed
+        decayed = active & past_peak & (i_new < 0.1) & ~crossed  # I_hold=0.1А
-        cut = crossed | local_min
+        local_min = active & past_peak & (i_new > i_old) & ~crossed & ~decayed
+        cut = crossed | local_min | decayed
         # zero-crossing: интерполяция к I=0; local_min: берём состояние минимума (old)
         frac = xp.where(crossed, i_old / (i_old - i_new + 1e-30), 0.0)
         exit_v = xp.where(crossed, _interp(v_old, v_new, frac), exit_v)
         exit_x = xp.where(crossed, _interp(x_old, x_new, frac), exit_x)
         exit_q = xp.where(crossed, _interp(q_old, q_new, frac), exit_q)
-        exit_v = xp.where(local_min, v_old, exit_v)
+        cut_at_state = local_min | decayed  # обрыв в текущем состоянии (не интерп.)
-        exit_x = xp.where(local_min, x_old, exit_x)
+        exit_v = xp.where(cut_at_state, v_old, exit_v)
-        exit_q = xp.where(local_min, q_old, exit_q)
+        exit_x = xp.where(cut_at_state, x_old, exit_x)
-        # остаточная энергия катушки при обрыве на минимуме -> в потери (freewheel)
+        exit_q = xp.where(cut_at_state, q_old, exit_q)
+        # остаточная энергия катушки при обрыве -> в потери (freewheel)
         residual = _magnetic_energy(xp, x_old, i_old, params)
-        energy_diss = energy_diss + xp.where(local_min, residual, 0.0)
+        energy_diss = energy_diss + xp.where(cut_at_state, residual, 0.0)
         committed = committed | cut  # валидная коммутация
 
         # стиффный конфиг «взорвал» fixed-step (inf/nan) -> стоп, НЕ реализуем

src/gausse/gpu/batch_sweep.py

@@ -73,6 +73,10 @@ def _prepare_stage(stage, projectile, cur_x_global, cur_v, coil_center_global):
     if cur_v <= 1e-6:
         return None, "снаряд остановился/пошёл назад"
     phys = build_stage_physics(stage, projectile)
+    # вырожденный контур (см. sim/stage.py) — тот же честный быстрый отказ
+    tau_worst = phys.inductance_model.l_air_h / (phys.r_total_ohm + phys.r_eddy_coeff_ohm + 1e-12)
+    if tau_worst < 1e-7:
+        return None, "вырожденная катушка: L/R < 100 нс (вне области модели)"
     sensor_global = coil_center_global - stage.sensor_to_coil_distance_m
     if sensor_global < cur_x_global - 1e-9:
         return None, "датчик позади снаряда (не сработает)"

src/gausse/optim/objective.py

@@ -21,7 +21,10 @@ from gausse.physics.inductance import air_core_inductance_wheeler, winding_geome
 from gausse.sim.coilgun import CoilgunResult, run_coilgun
 from gausse.storage.schema import RunRecord
 
-MODEL_VERSION = "gausse-physics-v2"  # v2: скин/близость (Доуэлл), трение+воздух, паспортные импульсные токи ключей
+# v3: коэффициент заполнения железа (fill factor) — закрыт эксплойт «тонкий
+# снаряд в толстой катушке»; полная формула отражения вихревого контура;
+# гард вырожденных катушек. v2: Доуэлл, трение+воздух, импульсные токи ключей.
+MODEL_VERSION = "gausse-physics-v3"
 
 
 @dataclass

src/gausse/physics/circuit.py

@@ -35,9 +35,18 @@ class StageCircuitParams:
     drag_coeff_n_per_mps2: float = 0.0
 
 
+# Сглаживание сухого трения около v=0: sign(v) -> tanh(v/ε). РАЗРЫВНАЯ сила
+# (чистый sign) ломает адаптивный решатель: у v≈0 шаг дробится бесконечно
+# (chattering), одна конфигурация считается десятки минут. tanh — стандартная
+# регуляризация (модель залипания); при |v| >> ε неотличима от sign.
+FRICTION_SMOOTHING_V_MPS = 0.01
+
+
 def retarding_force_n(params: StageCircuitParams, v: float) -> float:
-    """Тормозящая сила (трение + воздух) со знаком ПРОТИВ скорости."""
+    """Тормозящая сила (трение + воздух), гладкая по v, против скорости."""
-    return float(np.sign(v)) * (params.retard_const_n + params.drag_coeff_n_per_mps2 * v * v)
+    return float(np.tanh(v / FRICTION_SMOOTHING_V_MPS)) * (
+        params.retard_const_n + params.drag_coeff_n_per_mps2 * v * v
+    )
 
 
 def effective_resistance(inductance_model: CoilInductanceModel, params: StageCircuitParams, x) -> float:
@@ -91,3 +100,28 @@ def current_local_min_event(t: float, state: np.ndarray, inductance_model, param
 
 current_local_min_event.terminal = True
 current_local_min_event.direction = 1.0
+
+# Ток удержания ключа: тиристор сам закрывается, когда ток после пика падает
+# ниже I_H (типично 30-100 мА; берём консервативно 0.1 А). Без этого события
+# передемпфированный разряд (ток -> 0 асимптотически, нуля не пересекает)
+# заставляет решатель молотить весь временной бюджет микросекундными шагами —
+# одна конфигурация считалась минутами. Обрыв при I<I_hold — физика, не хак.
+HOLDING_CURRENT_A = 0.1
+
+
+def make_current_decay_event(i_hold_a: float = HOLDING_CURRENT_A):
+    """Терминальное событие: ток после пика упал ниже тока удержания ключа."""
+    peak = 0.0
+
+    def event(t: float, state: np.ndarray, *_args) -> float:
+        nonlocal peak
+        i = float(state[I])
+        if i > peak:
+            peak = i
+        if peak < 2.0 * i_hold_a:  # импульса ещё толком не было
+            return 1.0
+        return i - i_hold_a
+
+    event.terminal = True
+    event.direction = -1.0
+    return event

src/gausse/physics/inductance.py

@@ -107,10 +107,15 @@ class CoilInductanceModel:
     smoothing_width_m: float
     total_turns: int = 0
     b_sat_tesla: float = 1e9  # по умолчанию насыщение отключено (для геом. тестов)
+    # Доля сечения катушки, занятая железом: A_снаряда/A_среднего витка (≤1).
+    # λ = n·(B_fe·A_fe + B_возд·(A−A_fe)) ⇒ L = L_air·(1+(μ−1)·fill·overlap).
+    # БЕЗ этого тонкий снаряд в толстой катушке получал усиление как будто
+    # железо заполняет ВСЁ сечение — оптимизатор эксплуатировал (КПД «82%»).
+    iron_fill_factor: float = 1.0
 
     @property
     def l_iron_coeff(self) -> float:
-        return self.l_air_h * (self.mu_eff - 1.0)
+        return self.l_air_h * (self.mu_eff - 1.0) * self.iron_fill_factor
 
     @property
     def saturation_current_a(self) -> float:

src/gausse/physics/losses.py

@@ -46,7 +46,13 @@ def eddy_reflected_resistance_ohm(
     r_eddy_loop = 2 * math.pi * slug_resistivity_ohm_m / slug_length_m
     if r_eddy_loop <= 0:
         return 0.0
-    return (char_omega_rad_s * mutual) ** 2 / r_eddy_loop
+    # собственная индуктивность вихревого контура (виток-соленоид длиной снаряда):
+    # полная формула отражения R = (ωM)²·R₂/(R₂²+(ωL₂)²) вместо резистивно-
+    # доминированного приближения (ωM)²/R₂ — иначе на высоких ω (вырожденные
+    # короткие катушки) отражённое сопротивление нефизично улетало в сотни кОм.
+    l_eddy_loop = MU_0 * a_slug / slug_length_m
+    om = char_omega_rad_s
+    return (om * mutual) ** 2 * r_eddy_loop / (r_eddy_loop**2 + (om * l_eddy_loop) ** 2)
 
 
 def hysteresis_energy_j(slug_radius_m: float, slug_length_m: float, n_cycles: float = 1.0) -> float:

src/gausse/sim/stage.py

@@ -122,10 +122,14 @@ def _switch_equivalent_resistance_ohm(
 
 
 def _coast_derivatives(t: float, state: np.ndarray, retard_const_n: float, drag_coeff: float, mass_kg: float) -> list[float]:
-    """Подлёт к датчику: не баллистика в вакууме, а с трением о трубку и воздухом."""
+    """Подлёт к датчику: не баллистика в вакууме, а с трением о трубку и воздухом.
+
+    Трение сглажено tanh'ом (см. circuit.FRICTION_SMOOTHING_V_MPS) — разрывный
+    sign(v) заставляет адаптивный решатель бесконечно дробить шаг у v≈0.
+    """
     v = state[1]
-    retard = retard_const_n + drag_coeff * v * v
+    retard = (retard_const_n + drag_coeff * v * v) * math.tanh(v / circuit.FRICTION_SMOOTHING_V_MPS)
-    return [v, -math.copysign(retard, v) / mass_kg if v != 0.0 else 0.0]
+    return [v, -retard / mass_kg]
 
 
 def _stall_event(t: float, state: np.ndarray, *_args) -> float:
@@ -185,6 +189,9 @@ def build_stage_physics(stage: StageConfig, projectile: ProjectileConfig) -> Sta
         coil_length_m=geometry.coil_length_m,
         char_omega_rad_s=char_omega,
     )
+    # доля сечения катушки, реально занятая железом (см. inductance.py):
+    # тонкий снаряд в толстой катушке не может усиливать ВЕСЬ поток
+    fill = min((projectile.diameter_m / 2) ** 2 / geometry.mean_radius_m**2, 1.0)
     inductance_model = CoilInductanceModel(
         l_air_h=l_air_h,
         coil_length_m=geometry.coil_length_m,
@@ -193,6 +200,7 @@ def build_stage_physics(stage: StageConfig, projectile: ProjectileConfig) -> Sta
         smoothing_width_m=wire_od_m,
         total_turns=geometry.total_turns,
         b_sat_tesla=projectile.material.b_sat_tesla,
+        iron_fill_factor=fill,
     )
     mass_kg = projectile.mass_kg
     frontal_area_m2 = math.pi * (projectile.diameter_m / 2) ** 2
@@ -233,6 +241,21 @@ def run_stage(
     inductance_model = phys.inductance_model
     wire_od_m = stage.wire.insulation_od_mm / 1000
 
+    # Вырожденный контур — вне области применимости модели: при постоянной
+    # времени L/R (с учётом вихревых при полном перекрытии) короче ~100 нс
+    # (например, катушка в несколько витков длиной пару мм) ОДУ становится
+    # неинтегрируемо жёсткой (решатель молотит минутами), а сама инженерная
+    # модель индуктивности/вихрей там уже не имеет смысла. Честный быстрый отказ.
+    tau_worst_s = inductance_model.l_air_h / (r_total_ohm + r_eddy_coeff + 1e-12)
+    if tau_worst_s < 1e-7:
+        return StageResult(
+            feasible=False,
+            reason=(
+                f"вырожденная катушка: постоянная времени контура {tau_worst_s*1e9:.1f} нс "
+                f"(< 100 нс) — вне области применимости модели"
+            ),
+        )
+
     x_sensor_m = -stage.sensor_to_coil_distance_m
 
     if stage.sensor.kind == "inductive":
@@ -294,13 +317,18 @@ def run_stage(
         (0.0, stage.discharge_time_budget_s),
         [q0, 0.0, x_fire_m, v_fire_mps],
         args=(inductance_model, circuit_params),
-        events=(circuit.zero_current_crossing_event, circuit.current_local_min_event),
+        events=(
+            circuit.zero_current_crossing_event,
+            circuit.current_local_min_event,
+            # ключ закрывается сам при токе ниже удержания (см. circuit.py)
+            circuit.make_current_decay_event(),
+        ),
         rtol=1e-8,
         atol=1e-11,
     )
     # какое из событий оборвало разряд первым
     term_candidates = []
-    for ev_idx in range(2):
+    for ev_idx in range(3):
         if len(discharge_sol.t_events[ev_idx]) > 0:
             term_candidates.append((discharge_sol.t_events[ev_idx][0], discharge_sol.y_events[ev_idx][0]))
     if not term_candidates:

tests/test_realism_v2.py

@@ -145,3 +145,78 @@ def test_pulse_limit_fallback_without_datasheet():
     sw = replace(DB.switches[0], pulse_current_a=None)
     # без даташита — консервативный множитель по типу ключа
     assert switch_pulse_limit_a(sw) == sw.max_current_a * {"SCR": 10.0, "MOSFET": 4.0, "IGBT": 3.0}[sw.kind]
+
+
+def test_friction_force_smooth_at_zero_velocity():
+    """Регрессия на завис решателя: сила трения обязана быть ГЛАДКОЙ в v=0.
+
+    Разрывный sign(v) заставлял solve_ivp бесконечно дробить шаг у v≈0
+    (одна конфигурация считалась 35+ минут на сервере).
+    """
+    from gausse.physics.circuit import StageCircuitParams, retarding_force_n
+
+    p = StageCircuitParams(
+        capacitance_f=1e-4, r_total_ohm=0.1, mass_kg=0.01,
+        retard_const_n=0.034, drag_coeff_n_per_mps2=2e-5,
+    )
+    # непрерывность: около нуля сила ~ 0, нечётная, без скачка
+    assert abs(retarding_force_n(p, 1e-6)) < 1e-5
+    assert abs(retarding_force_n(p, 1e-6) + retarding_force_n(p, -1e-6)) < 1e-12
+    # вдали от нуля выходит на полное трение
+    assert retarding_force_n(p, 1.0) > 0.9 * p.retard_const_n
+
+
+def test_iron_fill_factor_punishes_thin_slug_in_fat_coil():
+    """Регрессия на эксплойт «КПД 82%»: тонкий снаряд в толстой катушке
+    не может усиливать весь поток — только долю сечения, которую занимает."""
+    from gausse.physics.inductance import CoilInductanceModel
+
+    base = dict(l_air_h=1e-3, coil_length_m=0.05, slug_length_m=0.02,
+                mu_eff=40.0, smoothing_width_m=0.001, total_turns=300)
+    full = CoilInductanceModel(**base, iron_fill_factor=1.0)
+    thin = CoilInductanceModel(**base, iron_fill_factor=0.02)  # ⌀4мм в ⌀29мм
+    # сила при том же токе в 50 раз меньше — халявы больше нет
+    f_full = float(full.force_newtons(-0.03, 100.0))
+    f_thin = float(thin.force_newtons(-0.03, 100.0))
+    assert f_thin < 0.03 * f_full
+
+
+def test_degenerate_coil_fails_fast():
+    """Регрессия на завис 540с: вырожденная катушка бракуется мгновенно."""
+    import time
+    from dataclasses import replace
+
+    stage = make_stage()
+    degenerate = replace(stage, turns_per_layer=2, layers=2)  # 4 витка, ~2мм
+    t0 = time.time()
+    res = run_stage(-0.06, 3.0, degenerate, make_projectile())
+    assert time.time() - t0 < 5.0, "вырожденный конфиг должен отваливаться быстро"
+    # либо честно отбракован гардом, либо честно посчитан быстро — оба исхода ок
+    if not res.feasible:
+        assert res.reason
+
+
+def test_eddy_reflection_bounded_at_high_frequency():
+    """Полная формула отражения (R₂²+(ωL₂)² в знаменателе): на высоких ω
+    заметно ниже старого резистивно-доминированного приближения, на низких —
+    сходится к нему. Совсем вырожденные катушки добивает гард L/R в run_stage."""
+    import math as _m
+
+    from gausse.physics.constants import MU_0
+    from gausse.physics.losses import EDDY_COUPLING_FACTOR, eddy_reflected_resistance_ohm
+
+    kwargs = dict(slug_radius_m=0.002, slug_length_m=0.03,
+                  slug_resistivity_ohm_m=1.6e-7, mu_eff=259.0,
+                  total_turns=5, coil_length_m=0.002)
+
+    def resistive_approx(om):
+        a = _m.pi * kwargs["slug_radius_m"] ** 2
+        m = EDDY_COUPLING_FACTOR * MU_0 * kwargs["mu_eff"] * (5 / 0.002) * a
+        r2 = 2 * _m.pi * kwargs["slug_resistivity_ohm_m"] / kwargs["slug_length_m"]
+        return (om * m) ** 2 / r2
+
+    om_hi, om_lo = 1.25e5, 1e2
+    r_hi = eddy_reflected_resistance_ohm(**kwargs, char_omega_rad_s=om_hi)
+    r_lo = eddy_reflected_resistance_ohm(**kwargs, char_omega_rad_s=om_lo)
+    assert r_hi < 0.5 * resistive_approx(om_hi)  # высокие ω: ограничена
+    assert abs(r_lo - resistive_approx(om_lo)) < 0.05 * resistive_approx(om_lo)  # низкие: сходится