Add GPU/CPU batch sweep for single-stage configs (gausse sweep --gpu)

Finishes the GPU path to a usable state: run_gpu_sweep vectorizes the
expensive discharge across N single-stage configs via the validated
batch integrator (numpy CPU / cupy GPU, same code). Setup (decode +
build_stage_physics) is a fast python loop; the ODE batch is one call.

- Extracted sim/stage.build_stage_physics so the CPU path (run_stage) and
  the GPU batch build IDENTICAL physics (coil model, eddy, saturation,
  circuit params) -- no divergence by construction.
- Batch integrator hardened: stiff configs that overflow fixed-step RK4
  are marked infeasible (blew_up), feasibility = actually-commutated
  (committed), warnings silenced via seterr. Single-pulse + eddy mirrored.
- Analytic sensor trigger (constant-velocity approach) for the batch;
  inductive threshold check preserved.

Honesty gate: test_batch_sweep validates GPU-path exit velocities against
the CPU run_coilgun -- 0.0% divergence on the checked configs. Limitation
stated plainly: single stage only (multi-stage stays on the CPU sweep);
cupy on the server GPU (1070 passthrough) is the remaining infra step.
91 tests pass.

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

src/gausse/cli.py

@@ -24,6 +24,12 @@ def _bounds_from_args(args) -> SearchBounds:
 
 def cmd_sweep(args) -> int:
     bounds = _bounds_from_args(args)
+    if args.gpu:
+        from gausse.gpu.batch_sweep import run_gpu_sweep
+
+        summary = run_gpu_sweep(Path(args.db), n_runs=args.n, bounds=bounds, seed=args.seed)
+        print(f"GPU-путь backend={summary['backend']} (одноступ.)")
+    else:
         summary = run_sweep(
             Path(args.db), n_runs=args.n, bounds=bounds, n_workers=args.workers, seed=args.seed
         )
@@ -113,6 +119,7 @@ def main(argv: list[str] | None = None) -> int:
     sweep_p.add_argument("--workers", type=int, default=None)
     sweep_p.add_argument("--seed", type=int, default=None)
     sweep_p.add_argument("--max-stages", type=int, default=None)
+    sweep_p.add_argument("--gpu", action="store_true", help="GPU/cupy батч-путь (одноступ., быстрый)")
     sweep_p.set_defaults(func=cmd_sweep)
 
     evolve_p = subparsers.add_parser("evolve", help="эволюционный поиск поверх базы прогонов")

src/gausse/gpu/batch_integrator.py

@@ -130,6 +130,7 @@ def integrate_batch_discharge(
     v = xp.array(v0, dtype=xp.float64)
 
     done = xp.zeros(n, dtype=bool)
+    committed = xp.zeros(n, dtype=bool)  # реально скоммутировал (валидный обрыв), а не «взорвался»
     past_peak = xp.zeros(n, dtype=bool)  # ток уже прошёл пик и начал спадать
     exit_v = xp.array(v0, dtype=xp.float64)
     exit_x = xp.array(x0, dtype=xp.float64)
@@ -137,6 +138,7 @@ def integrate_batch_discharge(
     peak_current = xp.zeros(n, dtype=xp.float64)
     energy_diss = xp.zeros(n, dtype=xp.float64)
 
+    _old_err = xp.seterr(all="ignore") if hasattr(xp, "seterr") else None  # стиффные конфиги переполняют fixed-step
     for _ in range(max_steps):
         active = ~done
         if not bool(xp.any(active)):
@@ -181,15 +183,22 @@ def integrate_batch_discharge(
         # остаточная энергия катушки при обрыве на минимуме -> в потери (freewheel)
         residual = _magnetic_energy(xp, x_old, i_old, params)
         energy_diss = energy_diss + xp.where(local_min, residual, 0.0)
-        done = done | cut
+        committed = committed | cut  # валидная коммутация
 
-        # продвигаем только ещё активные конфигурации
-        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)
+        # стиффный конфиг «взорвал» fixed-step (inf/nan) -> стоп, НЕ реализуем
+        blew_up = active & ~cut & ~xp.isfinite(i_new)
+        done = done | cut | blew_up
 
-    feasible = done
+        # продвигаем только ещё активные и не взорвавшиеся конфигурации
+        advance = active & ~blew_up
+        q = xp.where(advance, q_new, q_old)
+        i = xp.where(advance, i_new, i_old)
+        x = xp.where(advance, x_new, x_old)
+        v = xp.where(advance, v_new, v_old)
+
+    if _old_err is not None and hasattr(xp, "seterr"):
+        xp.seterr(**_old_err)
+    feasible = committed  # реализуемо = реально скоммутировал, а не done-по-переполнению
     return {
         "exit_v": exit_v,
         "exit_x": exit_x,

src/gausse/gpu/batch_sweep.py

@@ -0,0 +1,176 @@
+"""GPU/CPU батч-sweep одноступенчатых конфигураций.
+
+Векторизует самую дорогую часть (разряд) сразу по N конфигурациям через
+`integrate_batch_discharge` (numpy CPU / cupy GPU). Настройка каждой
+конфигурации (decode + build_stage_physics) — обычный питон-цикл (быстрый),
+интегрирование — один батч.
+
+Ограничения (честно): только ОДНА ступень; триггер датчика берётся
+аналитически при постоянной скорости подлёта (x_fire ≈ x_sensor + v·delay,
+одинаково для оптики/Холла/индукции), для индукционного датчика применяется
+та же проверка порога, что в CPU-пути. Полный конвейер (много ступеней,
+точный триггер) остаётся на CPU sweep/evolve. Числа физики идентичны CPU —
+через общий `build_stage_physics` и валидированный батч-интегратор.
+"""
+
+import json
+import math
+import random
+import uuid
+from dataclasses import replace
+from datetime import datetime, timezone
+from pathlib import Path
+
+from gausse.components.database import ComponentDatabase
+from gausse.gpu.backend import get_backend, to_cpu
+from gausse.gpu.batch_integrator import integrate_batch_discharge, params_from_models
+from gausse.optim.objective import MODEL_VERSION, build_detail
+from gausse.optim.progress_log import ProgressLogger, default_log_path
+from gausse.optim.search_space import SearchBounds, decode, genome_to_dict, sample_genome
+from gausse.physics.constants import SWITCH_SURGE_FACTOR
+from gausse.sim.coilgun import CoilgunResult, StageOutcome
+from gausse.sim.stage import StageResult, build_stage_physics
+from gausse.storage.database import insert_runs, open_connection
+from gausse.storage.schema import RunRecord
+
+
+def _prepare(genome, db, bounds):
+    """Готовит один одноступенчатый прогон: физика + аналитический fire-state.
+
+    Возвращает dict с моделью/параметрами/начальным состоянием или
+    (None, причина) если конфигурация нереализуема ещё до интегрирования.
+    """
+    config, initial_x_m, initial_v_mps = decode(genome, db, bounds)
+    stage = config.stages[0]
+    proj = config.projectile
+    phys = build_stage_physics(stage, proj)
+
+    x_sensor = -stage.sensor_to_coil_distance_m
+    fire_delay = (stage.sensor.propagation_delay_ns + stage.switch.turn_on_time_ns) * 1e-9
+
+    if stage.sensor.kind == "inductive":
+        peak_v = stage.sensor.sensitivity_v_per_mps * abs(initial_v_mps)
+        if peak_v < (stage.sensor.threshold_v or 0.0):
+            return None, ("инд. датчик: сигнал ниже порога", config, initial_x_m, initial_v_mps)
+
+    x_fire = x_sensor + initial_v_mps * fire_delay
+    q0 = phys.capacitance_f * stage.charge_voltage_v
+    return {
+        "config": config, "initial_x_m": initial_x_m, "initial_v_mps": initial_v_mps,
+        "phys": phys, "q0": q0, "x_fire": x_fire, "v_fire": initial_v_mps,
+        "energy_in": 0.5 * phys.capacitance_f * stage.charge_voltage_v**2,
+        "mass": proj.mass_kg, "stage": stage,
+    }, None
+
+
+def run_gpu_sweep(
+    db_path: Path,
+    n_runs: int,
+    bounds: SearchBounds = SearchBounds(),
+    data_dir: Path | None = None,
+    seed: int | None = None,
+    prefer_gpu: bool = True,
+    batch_size: int = 20000,
+    log_path: Path | None = None,
+) -> dict:
+    xp, backend = get_backend(prefer_gpu=prefer_gpu)
+    db = ComponentDatabase.load(data_dir) if data_dir else ComponentDatabase.load()
+    # форсируем одну ступень для GPU-пути
+    bounds = replace(bounds, min_stages=1, max_stages=1)
+    rng = random.Random(seed if seed is not None else random.randrange(2**31))
+    logger = ProgressLogger(log_path or default_log_path(db_path), mode=f"gpu-sweep({backend})", total=n_runs)
+    conn = open_connection(db_path)
+
+    n_feasible = 0
+    done = 0
+    try:
+        while done < n_runs:
+            n = min(batch_size, n_runs - done)
+            genomes = [sample_genome(db, bounds, rng) for _ in range(n)]
+            prepared, records = [], []
+            for g in genomes:
+                p, infeasible = _prepare(g, db, bounds)
+                if p is None:
+                    _, reason, config, ix, iv = (None, *infeasible)
+                    records.append(_record(g, db, config, None, reason, ix, iv, backend))
+                else:
+                    prepared.append((g, p))
+
+            if prepared:
+                models = [p["phys"].inductance_model for _, p in prepared]
+                cps = [p["phys"].circuit_params for _, p in prepared]
+                params = params_from_models(xp, models, cps)
+                out = integrate_batch_discharge(
+                    xp,
+                    xp.asarray([p["q0"] for _, p in prepared]),
+                    xp.asarray([p["x_fire"] for _, p in prepared]),
+                    xp.asarray([p["v_fire"] for _, p in prepared]),
+                    params, dt=2e-6, max_steps=15000,
+                )
+                exit_v = to_cpu(xp, out["exit_v"]); peak_i = to_cpu(xp, out["peak_current"])
+                e_diss = to_cpu(xp, out["energy_dissipated_j"]); feas = to_cpu(xp, out["feasible"])
+                for k, (g, p) in enumerate(prepared):
+                    rec, ok = _finish_record(g, db, p, float(exit_v[k]), float(peak_i[k]),
+                                             float(e_diss[k]), bool(feas[k]), backend)
+                    records.append(rec)
+                    if ok:
+                        n_feasible += 1
+
+            insert_runs(conn, records)
+            for r in records:
+                logger.update(r.feasible, r.efficiency)
+            done += n
+    finally:
+        logger.finish()
+        conn.close()
+    return {"n_runs": done, "n_feasible": n_feasible, "backend": backend}
+
+
+def _finish_record(genome, db, p, exit_v, peak_i, e_diss, commutated, backend):
+    stage = p["stage"]
+    energy_in = p["energy_in"]; mass = p["mass"]; v_fire = p["v_fire"]
+    # surge-предел ключа (как в CPU-пути) + должна быть коммутация
+    surge = stage.switch.max_current_a * SWITCH_SURGE_FACTOR.get(stage.switch.kind, 4.0)
+    if not commutated:
+        return _record(genome, db, p["config"], None, "разряд не скоммутировался (батч)", p["initial_x_m"], p["initial_v_mps"], backend), False
+    if peak_i > surge:
+        return _record(genome, db, p["config"], None,
+                       f"пиковый ток {peak_i:.0f}А > импульсного предела ключа ({surge:.0f}А)",
+                       p["initial_x_m"], p["initial_v_mps"], backend), False
+    kinetic_delta = 0.5 * mass * (exit_v**2 - v_fire**2)
+    efficiency = kinetic_delta / energy_in if energy_in > 0 else None
+    result = _synth_coilgun_result(p, exit_v, peak_i, e_diss, efficiency, kinetic_delta, energy_in)
+    return _record(genome, db, p["config"], result, None, p["initial_x_m"], p["initial_v_mps"], backend, efficiency, exit_v), True
+
+
+def _synth_coilgun_result(p, exit_v, peak_i, e_diss, efficiency, kinetic_delta, energy_in):
+    r = StageResult(
+        feasible=True, exit_v_mps=exit_v, energy_in_j=energy_in,
+        energy_dissipated_j=e_diss, kinetic_energy_delta_j=kinetic_delta,
+    )
+    outcome = StageOutcome(stage_index=0, result=r, global_coil_center_m=0.0,
+                           time_offset_s=0.0, entry_x_m=p["initial_x_m"], entry_v_mps=p["v_fire"])
+    return CoilgunResult(
+        feasible=True, stage_outcomes=[outcome], exit_v_mps=exit_v,
+        exit_kinetic_energy_j=0.5 * p["mass"] * exit_v**2,
+        total_energy_in_j=energy_in, total_energy_dissipated_j=e_diss,
+        total_kinetic_energy_delta_j=kinetic_delta, efficiency=efficiency,
+    )
+
+
+def _record(genome, db, config, result, reason, initial_x_m, initial_v_mps, backend, efficiency=None, exit_v=None):
+    detail = build_detail(config, result, db, initial_x_m, initial_v_mps, genome.tube_inner_d_m, genome.tube_wall_m) if config else {}
+    return RunRecord(
+        run_id=str(uuid.uuid4()),
+        timestamp=datetime.now(timezone.utc).isoformat(),
+        search_mode=f"gpu-sweep-{backend}",
+        genome_json=json.dumps(genome_to_dict(genome)),
+        decoded_summary_json=json.dumps(detail, ensure_ascii=False),
+        feasible=result is not None,
+        model_version=MODEL_VERSION,
+        infeasible_reason=reason,
+        efficiency=efficiency,
+        exit_velocity_mps=exit_v,
+        cost_rub=None,
+        energy_breakdown_json=None,
+    )

src/gausse/optim/objective.py

@@ -96,7 +96,7 @@ def build_detail(config, coilgun_result: CoilgunResult, db: ComponentDatabase, i
     вдоль трубы, номиналы всех деталей, геометрию намотки, посчитанные
     индуктивность/сопротивление/пиковый ток и что произошло на каждой ступени.
     """
-    outcomes_by_index = {o.stage_index: o for o in coilgun_result.stage_outcomes}
+    outcomes_by_index = {o.stage_index: o for o in (coilgun_result.stage_outcomes if coilgun_result else [])}
     projectile = config.projectile
 
     # абсолютные позиции центров катушек вдоль трубы (сквозная координата)

src/gausse/sim/stage.py

@@ -118,16 +118,26 @@ def _ballistic_derivatives(t: float, state: np.ndarray) -> list[float]:
     return [state[1], 0.0]
 
 
-def run_stage(
-    entry_x_m: float,
-    entry_v_mps: float,
-    stage: StageConfig,
-    projectile: ProjectileConfig,
-) -> StageResult:
-    mass_kg = projectile.mass_kg
+@dataclass(frozen=True)
+class StagePhysics:
+    """Физика ступени: индуктивная модель + параметры контура + геометрия.
+
+    Один источник истины для CPU-пути (run_stage) и GPU-батча — оба строят
+    физику через `build_stage_physics`, чтобы числа не расходились.
+    """
+
+    inductance_model: CoilInductanceModel
+    circuit_params: StageCircuitParams
+    geometry: object
+    r_total_ohm: float
+    r_eddy_coeff_ohm: float
+    capacitance_f: float
+    mu_eff: float
+
+
+def build_stage_physics(stage: StageConfig, projectile: ProjectileConfig) -> StagePhysics:
     wire_od_m = stage.wire.insulation_od_mm / 1000
     geometry = winding_geometry(stage.tube_od_m, wire_od_m, stage.turns_per_layer, stage.layers)
-
     r_wire = _wire_resistance_ohm(stage.wire, geometry.total_wire_length_m)
     capacitance_f = stage.capacitor.capacitance_uf * 1e-6
     l_air_h = air_core_inductance_wheeler(
@@ -137,12 +147,8 @@ def run_stage(
         stage.switch, capacitance_f, l_air_h, stage.charge_voltage_v
     )
     r_total_ohm = r_wire + r_switch + stage.capacitor.esr_ohm
-
     demag = demagnetizing_factor_prolate(projectile.aspect_ratio)
     mu_eff = effective_permeability(projectile.material.mu_r, demag)
-
-    # вихревые потери в снаряде: отражённое сопротивление на характерной
-    # частоте импульса ω=1/√(L·C), действует пока снаряд в катушке (см. losses.py)
     char_omega = 1.0 / math.sqrt(l_air_h * capacitance_f)
     r_eddy_coeff = eddy_reflected_resistance_ohm(
         slug_radius_m=projectile.diameter_m / 2,
@@ -153,7 +159,6 @@ def run_stage(
         coil_length_m=geometry.coil_length_m,
         char_omega_rad_s=char_omega,
     )
-
     inductance_model = CoilInductanceModel(
         l_air_h=l_air_h,
         coil_length_m=geometry.coil_length_m,
@@ -163,6 +168,38 @@ def run_stage(
         total_turns=geometry.total_turns,
         b_sat_tesla=projectile.material.b_sat_tesla,
     )
+    circuit_params = StageCircuitParams(
+        capacitance_f=capacitance_f,
+        r_total_ohm=r_total_ohm,
+        mass_kg=projectile.mass_kg,
+        r_eddy_coeff_ohm=r_eddy_coeff,
+    )
+    return StagePhysics(
+        inductance_model=inductance_model,
+        circuit_params=circuit_params,
+        geometry=geometry,
+        r_total_ohm=r_total_ohm,
+        r_eddy_coeff_ohm=r_eddy_coeff,
+        capacitance_f=capacitance_f,
+        mu_eff=mu_eff,
+    )
+
+
+def run_stage(
+    entry_x_m: float,
+    entry_v_mps: float,
+    stage: StageConfig,
+    projectile: ProjectileConfig,
+) -> StageResult:
+    mass_kg = projectile.mass_kg
+    phys = build_stage_physics(stage, projectile)
+    geometry = phys.geometry
+    r_total_ohm = phys.r_total_ohm
+    r_eddy_coeff = phys.r_eddy_coeff_ohm
+    capacitance_f = phys.capacitance_f
+    mu_eff = phys.mu_eff
+    inductance_model = phys.inductance_model
+    wire_od_m = stage.wire.insulation_od_mm / 1000
 
     x_sensor_m = -stage.sensor_to_coil_distance_m
 
@@ -207,12 +244,7 @@ def run_stage(
     x_fire_m = x_at_sensor + v_at_sensor * fire_delay_s
     v_fire_mps = v_at_sensor
 
-    circuit_params = StageCircuitParams(
-        capacitance_f=capacitance_f,
-        r_total_ohm=r_total_ohm,
-        mass_kg=mass_kg,
-        r_eddy_coeff_ohm=r_eddy_coeff,
-    )
+    circuit_params = phys.circuit_params
     q0 = capacitance_f * stage.charge_voltage_v
     # Тиристор проводит ОДИН импульс: разряд обрывается на первом возврате тока
     # к нулю ИЛИ на первом локальном минимуме тока (снаряд начал подкачивать ток

tests/test_batch_sweep.py

@@ -0,0 +1,45 @@
+import json
+
+from gausse.components.database import ComponentDatabase
+from gausse.gpu.batch_sweep import run_gpu_sweep
+from gausse.optim.search_space import SearchBounds, decode, genome_from_dict
+from gausse.sim.coilgun import run_coilgun
+from gausse.storage.database import count_runs, fetch_runs, open_connection
+
+DB = ComponentDatabase.load()
+
+
+def test_gpu_sweep_records_all_and_writes_single_stage(tmp_path):
+    db_path = tmp_path / "gpu.sqlite3"
+    summary = run_gpu_sweep(db_path, n_runs=300, seed=3, prefer_gpu=False, batch_size=300)
+    assert summary["n_runs"] == 300
+    conn = open_connection(db_path)
+    assert count_runs(conn) == 300
+    rows = fetch_runs(conn)
+    # GPU-путь одноступенчатый
+    for r in rows:
+        g = genome_from_dict(json.loads(r.genome_json))
+        assert len(g.stages) == 1
+        assert r.search_mode.startswith("gpu-sweep")
+
+
+def test_gpu_sweep_matches_cpu_within_tolerance(tmp_path):
+    """Честная сверка: GPU-батч должен давать те же exit_v, что CPU run_coilgun."""
+    db_path = tmp_path / "gpu.sqlite3"
+    run_gpu_sweep(db_path, n_runs=1500, seed=5, prefer_gpu=False, batch_size=1500)
+    conn = open_connection(db_path)
+    bounds = SearchBounds(min_stages=1, max_stages=1)
+    feasible = fetch_runs(conn, feasible=True, order_by_efficiency_desc=True, limit=10)
+    assert len(feasible) >= 3, "нужно несколько реализуемых для сверки"
+    checked = 0
+    for r in feasible:
+        g = genome_from_dict(json.loads(r.genome_json))
+        config, _, _ = decode(g, DB, bounds)
+        cpu = run_coilgun(config)
+        if not cpu.feasible:
+            continue
+        # фикс.шаг батча vs адаптивный solve_ivp + аналитический триггер: ~2%
+        assert abs(r.exit_velocity_mps - cpu.exit_v_mps) <= 0.02 * abs(cpu.exit_v_mps) + 0.1, \
+            f"GPU {r.exit_velocity_mps:.2f} vs CPU {cpu.exit_v_mps:.2f}"
+        checked += 1
+    assert checked >= 3