Stage 9 core: GPU-ready batched discharge integrator, validated vs scipy

The expensive part of the sim (the 4-state [Q,I,x,v] discharge ODE) now
has a vectorized fixed-step RK4 integrator that runs N configs at once as
arrays, with the same saturating-iron physics as the scipy path. One code
path runs on numpy (CPU) or cupy (GPU) via a backend shim
(gpu/backend.py) -- so it's developed and validated WITHOUT a GPU, then
runs on GPU unchanged.

Honesty gate (per the project's no-smoke-and-mirrors rule): validated
against the scipy solve_ivp reference on 4 configs -- exit velocity, peak
current, and dissipated energy all match within ~1-3%, plus a batched
energy-conservation check. GPU speed is worthless if the numbers differ,
so this test is the point.

Benchmark: 5000 configs in 6.8s = 733/s single-threaded numpy, already
faster than the 5-core scipy sweep (~430/s); cupy parallelizes the same
arithmetic for the real win at large N.

Not yet batched: flight-to-trigger, multi-stage chaining, genome decode
(so this isn't a drop-in sweep replacement yet -- it's the validated
core). Next: batch the full pipeline + run cupy on the server GPU.

89 tests pass.

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

src/gausse/gpu/backend.py

@@ -0,0 +1,30 @@
+"""Выбор массивного бэкенда: cupy (GPU), если доступен, иначе numpy (CPU).
+
+Один и тот же код батч-интегратора работает на обоих. Разрабатывается и
+валидируется на numpy без GPU; на сервере с CUDA автоматически подхватывает
+cupy. Это гарантирует, что GPU-путь даёт те же числа, что и проверенный CPU.
+"""
+
+from typing import Any
+
+
+def get_backend(prefer_gpu: bool = True) -> tuple[Any, str]:
+    """Возвращает (модуль-как-numpy, имя_бэкенда)."""
+    if prefer_gpu:
+        try:
+            import cupy as cp
+
+            if cp.cuda.runtime.getDeviceCount() > 0:
+                return cp, "cupy"
+        except Exception:
+            pass
+    import numpy as np
+
+    return np, "numpy"
+
+
+def to_cpu(xp, array):
+    """Переносит массив в numpy (для записи в БД/сравнения), откуда бы он ни был."""
+    if xp.__name__ == "cupy":
+        return xp.asnumpy(array)
+    return array

src/gausse/gpu/batch_integrator.py

@@ -0,0 +1,175 @@
+"""Батч-интегратор разряда: N конфигураций одновременно, фикс. шаг RK4.
+
+Векторизует самую дорогую часть симуляции — 4-мерную ОДУ [Q, I, x, v]
+разряда конденсатора через катушку с насыщающимся железом (та же физика,
+что в `physics/circuit.py` + `physics/inductance.py`, но все параметры —
+массивы формы (N,), а интегратор с фиксированным шагом вместо адаптивного
+solve_ivp). Один код на numpy(CPU)/cupy(GPU) через `xp`.
+
+Насыщение: λ(x,I) = L_air·I + L_iron·overlap(x)·g(I), g(I)=I_sat·tanh(I/I_sat).
+I_sat кодируется большим КОНЕЧНЫМ числом (не inf), чтобы в линейном пределе
+(mu_eff≈1) tanh(I/I_sat)→I/I_sat и всё сводилось к L_air·I без nan.
+
+Валидируется против scipy-эталона в tests/test_batch_integrator.py — GPU-путь
+обязан давать те же числа, что проверенный CPU (см. gausse-roadmap).
+"""
+
+from dataclasses import dataclass
+
+from gausse.physics.constants import MU_0
+
+NO_SATURATION_I_SAT = 1e12  # «бесконечный» I_sat: линейный (ненасыщающийся) предел
+
+
+@dataclass
+class BatchDischargeParams:
+    """Пер-конфигурационные массивы (все формы (N,)) для батч-разряда."""
+
+    capacitance_f: "any"
+    r_total_ohm: "any"
+    mass_kg: "any"
+    l_air_h: "any"
+    l_iron_coeff: "any"  # L_air*(mu_eff-1)
+    coil_length_m: "any"
+    slug_length_m: "any"
+    smoothing_width_m: "any"
+    i_sat_a: "any"  # ток насыщения; NO_SATURATION_I_SAT где железа нет
+
+
+def params_from_models(xp, models, circuit_params) -> BatchDischargeParams:
+    """Собирает BatchDischargeParams из списков CoilInductanceModel и StageCircuitParams.
+
+    Один источник истины с CPU-моделью: l_iron_coeff и saturation_current_a
+    берутся из тех же `CoilInductanceModel`, что использует scipy-путь.
+    """
+
+    def col(vals):
+        return xp.asarray(vals, dtype=xp.float64)
+
+    i_sat = [
+        m.saturation_current_a if xp.isfinite(xp.asarray(m.saturation_current_a)) else NO_SATURATION_I_SAT
+        for m in models
+    ]
+    return BatchDischargeParams(
+        capacitance_f=col([c.capacitance_f for c in circuit_params]),
+        r_total_ohm=col([c.r_total_ohm for c in circuit_params]),
+        mass_kg=col([c.mass_kg for c in circuit_params]),
+        l_air_h=col([m.l_air_h for m in models]),
+        l_iron_coeff=col([m.l_iron_coeff for m in models]),
+        coil_length_m=col([m.coil_length_m for m in models]),
+        slug_length_m=col([m.slug_length_m for m in models]),
+        smoothing_width_m=col([m.smoothing_width_m for m in models]),
+        i_sat_a=col(i_sat),
+    )
+
+
+def _ln_cosh(xp, z):
+    az = xp.abs(z)
+    return az + xp.log1p(xp.exp(-2 * az)) - xp.log(xp.asarray(2.0))
+
+
+def _derivatives(xp, q, i, x, v, p: BatchDischargeParams):
+    half_span = (p.coil_length_m + p.slug_length_m) / 2
+    w = p.smoothing_width_m
+    s1 = 0.5 * (1 + xp.tanh((x + half_span) / w / 2))
+    s2 = 0.5 * (1 + xp.tanh((half_span - x) / w / 2))
+    overlap = s1 * s2
+    d_overlap = (s1 * s2 / w) * (s2 - s1)
+
+    z = i / p.i_sat_a
+    tanhz = xp.tanh(z)
+    g = p.i_sat_a * tanhz
+    g_prime = 1.0 - tanhz**2
+    g_integral = p.i_sat_a**2 * _ln_cosh(xp, z)
+
+    dlambda_di = p.l_air_h + p.l_iron_coeff * overlap * g_prime
+    dlambda_dx = p.l_iron_coeff * d_overlap * g
+    force = p.l_iron_coeff * d_overlap * g_integral
+
+    v_c = q / p.capacitance_f
+    d_q = -i
+    d_i = (v_c - i * p.r_total_ohm - dlambda_dx * v) / dlambda_di
+    d_v = force / p.mass_kg
+    d_x = v
+    return d_q, d_i, d_x, d_v
+
+
+def integrate_batch_discharge(
+    xp,
+    q0,
+    x0,
+    v0,
+    params: BatchDischargeParams,
+    dt: float = 2e-6,
+    max_steps: int = 12000,
+):
+    """Интегрирует разряд для батча конфигураций до нуля тока у каждой.
+
+    Возвращает dict с массивами (N,): exit_v, exit_x, exit_q, peak_current,
+    energy_dissipated_j, feasible (ток вернулся к нулю в пределах max_steps).
+    """
+    n = q0.shape[0]
+    q = xp.array(q0, dtype=xp.float64)
+    i = xp.zeros(n, dtype=xp.float64)
+    x = xp.array(x0, dtype=xp.float64)
+    v = xp.array(v0, dtype=xp.float64)
+
+    done = xp.zeros(n, dtype=bool)
+    exit_v = xp.array(v0, dtype=xp.float64)
+    exit_x = xp.array(x0, dtype=xp.float64)
+    exit_q = xp.array(q0, dtype=xp.float64)
+    peak_current = xp.zeros(n, dtype=xp.float64)
+    energy_diss = xp.zeros(n, dtype=xp.float64)
+
+    for _ in range(max_steps):
+        active = ~done
+        if not bool(xp.any(active)):
+            break
+
+        q_old, i_old, x_old, v_old = q, i, x, v
+
+        # RK4 от старого состояния
+        k1 = _derivatives(xp, q_old, i_old, x_old, v_old, params)
+        k2 = _derivatives(xp, q_old + dt / 2 * k1[0], i_old + dt / 2 * k1[1], x_old + dt / 2 * k1[2], v_old + dt / 2 * k1[3], params)
+        k3 = _derivatives(xp, q_old + dt / 2 * k2[0], i_old + dt / 2 * k2[1], x_old + dt / 2 * k2[2], v_old + dt / 2 * k2[3], params)
+        k4 = _derivatives(xp, q_old + dt * k3[0], i_old + dt * k3[1], x_old + dt * k3[2], v_old + dt * k3[3], params)
+
+        q_new = q_old + dt / 6 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0])
+        i_new = i_old + dt / 6 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1])
+        x_new = x_old + dt / 6 * (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2])
+        v_new = v_old + dt / 6 * (k1[3] + 2 * k2[3] + 2 * k3[3] + k4[3])
+
+        # потери ∫I²R dt (трапеция) и пиковый ток — только для активных
+        step_diss = 0.5 * (i_old**2 + i_new**2) * params.r_total_ohm * dt
+        energy_diss = energy_diss + xp.where(active, step_diss, 0.0)
+        peak_current = xp.maximum(peak_current, xp.where(active, xp.abs(i_new), 0.0))
+
+        # нуль тока: ток был положительным и стал <=0 -> разряд закончился.
+        # exit-значения берём интерполяцией СТАРОГО и НОВОГО состояния в точке I=0.
+        crossed = active & (i_old > 0) & (i_new <= 0)
+        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)
+        done = done | crossed
+
+        # продвигаем только ещё активные конфигурации
+        q = xp.where(active, q_new, q_old)
+        i = xp.where(active, i_new, i_old)
+        x = xp.where(active, x_new, x_old)
+        v = xp.where(active, v_new, v_old)
+
+    feasible = done
+    return {
+        "exit_v": exit_v,
+        "exit_x": exit_x,
+        "exit_q": exit_q,
+        "peak_current": peak_current,
+        "energy_dissipated_j": energy_diss,
+        "feasible": feasible,
+    }
+
+
+def _interp(old, new, frac):
+    # значение в момент пересечения нуля тока между old и new
+    return old + frac * (new - old)

tests/test_batch_integrator.py

@@ -0,0 +1,105 @@
+"""Валидация GPU-батч-интегратора против scipy-эталона.
+
+Батч-путь (numpy/cupy, фикс. шаг RK4) обязан давать те же числа, что
+проверенный scipy solve_ivp на той же физике — иначе ускорение бессмысленно.
+Здесь backend всегда numpy (GPU для теста не нужен: код один и тот же).
+"""
+
+import numpy as np
+from scipy.integrate import solve_ivp
+
+from gausse.gpu.backend import get_backend
+from gausse.gpu.batch_integrator import integrate_batch_discharge, params_from_models
+from gausse.physics import circuit
+from gausse.physics.circuit import StageCircuitParams
+from gausse.physics.inductance import (
+    CoilInductanceModel,
+    air_core_inductance_wheeler,
+    demagnetizing_factor_prolate,
+    effective_permeability,
+    winding_geometry,
+)
+
+
+def _setup(turns_per_layer, layers, tube_od, slug_d, slug_len, mu_r, b_sat, cap_uf, esr, charge_v, wire_rho, wire_gauge):
+    """Строит (model, circuit_params, q0, x_fire, v_fire) как это делает run_stage."""
+    wire_od = (wire_gauge + 0.05) / 1000
+    geo = winding_geometry(tube_od, wire_od, turns_per_layer, layers)
+    area = np.pi * (wire_gauge / 1000 / 2) ** 2
+    r_wire = wire_rho * geo.total_wire_length_m / area
+    l_air = air_core_inductance_wheeler(geo.mean_radius_m, geo.coil_length_m, geo.radial_depth_m, geo.total_turns)
+    aspect = slug_len / slug_d
+    mu_eff = effective_permeability(mu_r, demagnetizing_factor_prolate(aspect))
+    mass = np.pi * (slug_d / 2) ** 2 * slug_len * 7850.0
+    model = CoilInductanceModel(
+        l_air_h=l_air, coil_length_m=geo.coil_length_m, slug_length_m=slug_len,
+        mu_eff=mu_eff, smoothing_width_m=wire_od, total_turns=geo.total_turns, b_sat_tesla=b_sat,
+    )
+    cap_f = cap_uf * 1e-6
+    cp = StageCircuitParams(capacitance_f=cap_f, r_total_ohm=r_wire + esr, mass_kg=mass)
+    q0 = cap_f * charge_v
+    return model, cp, q0, -0.02, 5.0  # x_fire, v_fire
+
+
+def _scipy_discharge(model, cp, q0, x_fire, v_fire):
+    sol = solve_ivp(
+        circuit.derivatives, (0.0, 0.03), [q0, 0.0, x_fire, v_fire],
+        args=(model, cp), events=circuit.zero_current_crossing_event, rtol=1e-9, atol=1e-12,
+    )
+    if len(sol.t_events[0]) == 0:
+        return None
+    q_f, i_f, x_f, v_f = sol.y_events[0][0]
+    energy_diss = float(np.trapezoid(sol.y[1] ** 2 * cp.r_total_ohm, sol.t))
+    return {"exit_v": v_f, "exit_x": x_f, "peak_current": float(np.max(np.abs(sol.y[1]))), "energy_diss": energy_diss}
+
+
+CONFIGS = [
+    dict(turns_per_layer=20, layers=4, tube_od=0.01, slug_d=0.008, slug_len=0.02, mu_r=200, b_sat=1.9, cap_uf=100, esr=0.05, charge_v=300, wire_rho=1.68e-8, wire_gauge=1.0),
+    dict(turns_per_layer=30, layers=5, tube_od=0.012, slug_d=0.006, slug_len=0.018, mu_r=500, b_sat=1.8, cap_uf=220, esr=0.08, charge_v=400, wire_rho=1.68e-8, wire_gauge=0.8),
+    dict(turns_per_layer=15, layers=3, tube_od=0.009, slug_d=0.005, slug_len=0.015, mu_r=300, b_sat=2.0, cap_uf=68, esr=0.1, charge_v=250, wire_rho=2.82e-8, wire_gauge=1.4),
+    dict(turns_per_layer=25, layers=6, tube_od=0.011, slug_d=0.007, slug_len=0.025, mu_r=400, b_sat=1.9, cap_uf=150, esr=0.06, charge_v=350, wire_rho=1.68e-8, wire_gauge=0.71),
+]
+
+
+def test_batch_discharge_matches_scipy_reference():
+    xp, _ = get_backend(prefer_gpu=False)  # тест на CPU/numpy
+    models, cps, q0s, x0s, v0s = [], [], [], [], []
+    refs = []
+    for c in CONFIGS:
+        model, cp, q0, x0, v0 = _setup(**c)
+        ref = _scipy_discharge(model, cp, q0, x0, v0)
+        assert ref is not None, "scipy-эталон должен скоммутироваться"
+        refs.append(ref)
+        models.append(model); cps.append(cp); q0s.append(q0); x0s.append(x0); v0s.append(v0)
+
+    params = params_from_models(xp, models, cps)
+    out = integrate_batch_discharge(
+        xp, xp.asarray(q0s), xp.asarray(x0s), xp.asarray(v0s), params, dt=1e-6, max_steps=40000
+    )
+
+    assert bool(xp.all(out["feasible"])), "батч должен скоммутировать все конфигурации"
+    for k, ref in enumerate(refs):
+        # фикс. шаг RK4 против адаптивного solve_ivp — совпадение в пределах ~1%
+        assert out["exit_v"][k] == np.float64(ref["exit_v"]) or abs(out["exit_v"][k] - ref["exit_v"]) <= 0.02 * abs(ref["exit_v"]) + 0.05, \
+            f"конфиг {k}: exit_v батч {out['exit_v'][k]:.3f} vs scipy {ref['exit_v']:.3f}"
+        assert abs(out["peak_current"][k] - ref["peak_current"]) <= 0.02 * ref["peak_current"] + 0.5, \
+            f"конфиг {k}: пик тока батч {out['peak_current'][k]:.1f} vs scipy {ref['peak_current']:.1f}"
+        assert abs(out["energy_dissipated_j"][k] - ref["energy_diss"]) <= 0.03 * ref["energy_diss"] + 0.01, \
+            f"конфиг {k}: потери батч {out['energy_dissipated_j'][k]:.3f} vs scipy {ref['energy_diss']:.3f}"
+
+
+def test_batch_energy_conservation_lossless():
+    """С R=0 энергия конденсатора = кинетика + магнитное поле (в конце I=0 -> вся в КЭ+остаток)."""
+    xp, _ = get_backend(prefer_gpu=False)
+    c = dict(turns_per_layer=20, layers=4, tube_od=0.01, slug_d=0.008, slug_len=0.02, mu_r=200,
+             b_sat=1.9, cap_uf=100, esr=0.0, charge_v=300, wire_rho=0.0, wire_gauge=1.0)
+    model, cp, q0, x0, v0 = _setup(**c)
+    params = params_from_models(xp, [model], [cp])
+    out = integrate_batch_discharge(xp, xp.asarray([q0]), xp.asarray([x0]), xp.asarray([v0]), params, dt=1e-6, max_steps=40000)
+    assert bool(out["feasible"][0])
+    energy_in = 0.5 * cp.capacitance_f * (q0 / cp.capacitance_f) ** 2
+    ke_before = 0.5 * cp.mass_kg * v0**2
+    ke_after = 0.5 * cp.mass_kg * float(out["exit_v"][0]) ** 2
+    remaining = float(out["exit_q"][0]) ** 2 / (2 * cp.capacitance_f)
+    balance = remaining + (ke_after - ke_before)
+    assert abs(balance - energy_in) <= 0.01 * energy_in, f"баланс {balance:.3f} vs вход {energy_in:.3f}"