diff --git a/src/gausse/cli.py b/src/gausse/cli.py index ce8ccc4..c7e1a0a 100644 --- a/src/gausse/cli.py +++ b/src/gausse/cli.py @@ -24,9 +24,15 @@ def _bounds_from_args(args) -> SearchBounds: def cmd_sweep(args) -> int: bounds = _bounds_from_args(args) - summary = run_sweep( - Path(args.db), n_runs=args.n, bounds=bounds, n_workers=args.workers, seed=args.seed - ) + 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 + ) rate = 100 * summary["n_feasible"] / summary["n_runs"] if summary["n_runs"] else 0.0 print(f"Прогнано {summary['n_runs']}, реализуемо {summary['n_feasible']} ({rate:.1f}%)") return 0 @@ -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="эволюционный поиск поверх базы прогонов") diff --git a/src/gausse/gpu/batch_integrator.py b/src/gausse/gpu/batch_integrator.py index de594cc..7e60322 100644 --- a/src/gausse/gpu/batch_integrator.py +++ b/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, diff --git a/src/gausse/gpu/batch_sweep.py b/src/gausse/gpu/batch_sweep.py new file mode 100644 index 0000000..e18a2bf --- /dev/null +++ b/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, + ) diff --git a/src/gausse/optim/objective.py b/src/gausse/optim/objective.py index 81a66fb..551614c 100644 --- a/src/gausse/optim/objective.py +++ b/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 # абсолютные позиции центров катушек вдоль трубы (сквозная координата) diff --git a/src/gausse/sim/stage.py b/src/gausse/sim/stage.py index 9c84249..6e49313 100644 --- a/src/gausse/sim/stage.py +++ b/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 # Тиристор проводит ОДИН импульс: разряд обрывается на первом возврате тока # к нулю ИЛИ на первом локальном минимуме тока (снаряд начал подкачивать ток diff --git a/tests/test_batch_sweep.py b/tests/test_batch_sweep.py new file mode 100644 index 0000000..c382c1b --- /dev/null +++ b/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