diff --git a/src/gausse/gpu/batch_sweep.py b/src/gausse/gpu/batch_sweep.py index 119dde0..0c78c82 100644 --- a/src/gausse/gpu/batch_sweep.py +++ b/src/gausse/gpu/batch_sweep.py @@ -23,7 +23,7 @@ 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.gpu.batch_integrator import NO_SATURATION_I_SAT, BatchDischargeParams, integrate_batch_discharge 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 @@ -36,13 +36,101 @@ from gausse.storage.schema import RunRecord _FEASIBLE = CoilgunResult(feasible=True, stage_outcomes=[]) +@dataclass(frozen=True) +class StagePre: + """Физика ступени, предвычисленная в ПЛОСКИЕ числа (пиклится дёшево). + + Считается один раз на геном (в воркере пула, параллельно) через тот же + `build_stage_physics` — один источник истины с CPU-путём. Главному потоку + остаётся только склейка массивов и запуск GPU-ядер: раньше он строил + физобъекты 4000 геномов сам и душил GPU (карта ждала ~половину времени). + """ + + l_air_h: float + l_iron_coeff: float + coil_length_m: float + slug_length_m: float + smoothing_width_m: float + i_sat_a: float + r_total_ohm: float + r_eddy_coeff_ohm: float + retard_const_n: float + drag_coeff_n: float + capacitance_f: float + q0: float + energy_in_j: float + sensor_to_coil_m: float + fire_delay_s: float + sensor_kind: str + sensor_sens_v_per_mps: float + sensor_threshold_v: float + pulse_limit_a: float + degenerate_reason: str | None + + +@dataclass(frozen=True) +class GenomePre: + stages: tuple # StagePre по ступеням + coil_centers: tuple # абсолютные центры катушек вдоль трубы + mass_kg: float + initial_x_m: float + initial_v_mps: float + + +def precompute_genome(genome, db: ComponentDatabase, bounds: SearchBounds) -> GenomePre: + """decode + build_stage_physics -> плоские числа (годится для пула процессов).""" + cfg, ix, iv = decode(genome, db, bounds) + centers = [0.0] + for gap in cfg.inter_stage_gaps_m: + centers.append(centers[-1] + gap) + stages = [] + for stage in cfg.stages: + phys = build_stage_physics(stage, cfg.projectile) + m = phys.inductance_model + cp_ = phys.circuit_params + tau_worst = m.l_air_h / (phys.r_total_ohm + phys.r_eddy_coeff_ohm + 1e-12) + i_sat = m.saturation_current_a + stages.append(StagePre( + l_air_h=m.l_air_h, + l_iron_coeff=m.l_iron_coeff, + coil_length_m=m.coil_length_m, + slug_length_m=m.slug_length_m, + smoothing_width_m=m.smoothing_width_m, + i_sat_a=i_sat if math.isfinite(i_sat) else NO_SATURATION_I_SAT, + r_total_ohm=phys.r_total_ohm, + r_eddy_coeff_ohm=phys.r_eddy_coeff_ohm, + retard_const_n=cp_.retard_const_n, + drag_coeff_n=cp_.drag_coeff_n_per_mps2, + capacitance_f=phys.capacitance_f, + q0=phys.capacitance_f * stage.charge_voltage_v, + energy_in_j=0.5 * phys.capacitance_f * stage.charge_voltage_v**2, + sensor_to_coil_m=stage.sensor_to_coil_distance_m, + fire_delay_s=(stage.sensor.propagation_delay_ns + stage.switch.turn_on_time_ns) * 1e-9, + sensor_kind=stage.sensor.kind, + sensor_sens_v_per_mps=getattr(stage.sensor, "sensitivity_v_per_mps", None) or 0.0, + sensor_threshold_v=getattr(stage.sensor, "threshold_v", None) or 0.0, + pulse_limit_a=switch_pulse_limit_a(stage.switch), + degenerate_reason=( + "вырожденная катушка: L/R < 100 нс (вне области модели)" if tau_worst < 1e-7 else None + ), + )) + return GenomePre( + stages=tuple(stages), coil_centers=tuple(centers), + mass_kg=cfg.projectile.mass_kg, initial_x_m=ix, initial_v_mps=iv, + ) + + +def precompute_worker(genome) -> GenomePre: + """Вариант для пула (ProcessPoolExecutor c worker_context.init_worker).""" + from gausse.optim import worker_context + + return precompute_genome(genome, worker_context.db, worker_context.bounds) + + @dataclass class _State: genome: object - config: object - initial_x_m: float - initial_v_mps: float - coil_centers: list # абсолютные центры катушек вдоль трубы + pre: GenomePre cur_x: float # текущая глобальная координата снаряда cur_v: float # текущая скорость energy_in_j: float = 0.0 @@ -69,36 +157,30 @@ def _coast_velocity(v0: float, distance_m: float, retard_const_n: float, drag_co return math.sqrt(u) if u > 0 else 0.0 -def _prepare_stage(stage, projectile, cur_x_global, cur_v, coil_center_global): - """Готовит разряд ступени: физика + аналитический fire-state из текущего (x,v).""" +def _prepare_stage(sp: StagePre, mass_kg: float, cur_x_global, cur_v, coil_center_global): + """Аналитический fire-state ступени из текущего (x,v) — только арифметика.""" 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 sp.degenerate_reason: + return None, sp.degenerate_reason + sensor_global = coil_center_global - sp.sensor_to_coil_m if sensor_global < cur_x_global - 1e-9: return None, "датчик позади снаряда (не сработает)" # подлёт до датчика — с трением о трубку и воздухом (та же физика, что на CPU) - cp_ = phys.circuit_params v_at_sensor = _coast_velocity( - cur_v, sensor_global - cur_x_global, cp_.retard_const_n, cp_.drag_coeff_n_per_mps2, projectile.mass_kg + cur_v, sensor_global - cur_x_global, sp.retard_const_n, sp.drag_coeff_n, mass_kg ) if v_at_sensor <= 1e-3: return None, "снаряд остановлен трением о трубку, не долетев до датчика" - if stage.sensor.kind == "inductive": - if stage.sensor.sensitivity_v_per_mps * v_at_sensor < (stage.sensor.threshold_v or 0.0): + if sp.sensor_kind == "inductive": + if sp.sensor_sens_v_per_mps * v_at_sensor < sp.sensor_threshold_v: return None, "инд. датчик: сигнал ниже порога" - fire_delay = (stage.sensor.propagation_delay_ns + stage.switch.turn_on_time_ns) * 1e-9 - x_fire_global = sensor_global + v_at_sensor * fire_delay + x_fire_global = sensor_global + v_at_sensor * sp.fire_delay_s x_fire_local = x_fire_global - coil_center_global - q0 = phys.capacitance_f * stage.charge_voltage_v return { - "phys": phys, "q0": q0, "x_fire": x_fire_local, "v_fire": v_at_sensor, - "energy_in": 0.5 * phys.capacitance_f * stage.charge_voltage_v**2, - "coil_center": coil_center_global, "stage": stage, "mass": projectile.mass_kg, + "sp": sp, "x_fire": x_fire_local, "v_fire": v_at_sensor, + "coil_center": coil_center_global, "mass": mass_kg, }, None @@ -106,14 +188,14 @@ def _simulate_states(xp, states: list) -> None: """Раунд-за-раундом прогоняет разряды всех живых ступеней батчами (мутирует states).""" if not states: return - max_ns = max(len(s.config.stages) for s in states) + max_ns = max(len(s.pre.stages) for s in states) for s_idx in range(max_ns): prepared = [] # (state, prep) for st in states: - if not st.alive or len(st.config.stages) <= s_idx: + if not st.alive or len(st.pre.stages) <= s_idx: continue p, reason = _prepare_stage( - st.config.stages[s_idx], st.config.projectile, st.cur_x, st.cur_v, st.coil_centers[s_idx] + st.pre.stages[s_idx], st.pre.mass_kg, st.cur_x, st.cur_v, st.pre.coil_centers[s_idx] ) if p is None: st.alive = False @@ -124,69 +206,90 @@ def _simulate_states(xp, states: list) -> None: if not prepared: continue - 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) + def col(get): + return xp.asarray([get(p["sp"], p) for _, p in prepared], dtype=xp.float64) + + params = BatchDischargeParams( + capacitance_f=col(lambda sp, p: sp.capacitance_f), + r_total_ohm=col(lambda sp, p: sp.r_total_ohm), + mass_kg=col(lambda sp, p: p["mass"]), + l_air_h=col(lambda sp, p: sp.l_air_h), + l_iron_coeff=col(lambda sp, p: sp.l_iron_coeff), + coil_length_m=col(lambda sp, p: sp.coil_length_m), + slug_length_m=col(lambda sp, p: sp.slug_length_m), + smoothing_width_m=col(lambda sp, p: sp.smoothing_width_m), + i_sat_a=col(lambda sp, p: sp.i_sat_a), + r_eddy_coeff_ohm=col(lambda sp, p: sp.r_eddy_coeff_ohm), + retard_const_n=col(lambda sp, p: sp.retard_const_n), + drag_coeff_n=col(lambda sp, p: sp.drag_coeff_n), + ) 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]), + col(lambda sp, p: sp.q0), + col(lambda sp, p: p["x_fire"]), + col(lambda sp, p: p["v_fire"]), params, dt=2e-6, max_steps=15000, ) exit_v = to_cpu(xp, out["exit_v"]); exit_x = to_cpu(xp, out["exit_x"]) peak_i = to_cpu(xp, out["peak_current"]); feas = to_cpu(xp, out["feasible"]) for k, (st, p) in enumerate(prepared): - stage = p["stage"] + sp = p["sp"] if not bool(feas[k]): st.alive = False st.reason = f"ступень {s_idx}: разряд не скоммутировался" st.failed_stage_index = s_idx continue - surge = switch_pulse_limit_a(stage.switch) - if float(peak_i[k]) > surge: + if float(peak_i[k]) > sp.pulse_limit_a: st.alive = False - st.reason = f"ступень {s_idx}: пиковый ток {float(peak_i[k]):.0f}А > импульсного предела ключа ({surge:.0f}А)" + st.reason = f"ступень {s_idx}: пиковый ток {float(peak_i[k]):.0f}А > импульсного предела ключа ({sp.pulse_limit_a:.0f}А)" st.failed_stage_index = s_idx continue ev = float(exit_v[k]) st.kinetic_delta_j += 0.5 * p["mass"] * (ev**2 - p["v_fire"] ** 2) - st.energy_in_j += p["energy_in"] + st.energy_in_j += sp.energy_in_j st.cur_v = ev st.cur_x = p["coil_center"] + float(exit_x[k]) def evaluate_genomes_gpu( - xp, genomes: list, db: ComponentDatabase, bounds: SearchBounds, build_details: bool = True + xp, genomes: list, db: ComponentDatabase, bounds: SearchBounds, + build_details: bool = True, executor=None, ) -> list: """Оценка списка геномов ОДНИМ батчем (для эволюции): та же физика и та же формула фитнеса, что в objective.evaluate (КПД либо -1+доля пройденных ступеней), но разряды всех геномов интегрируются вместе на GPU/numpy. build_details=False — не строить detail/стоимость (фитнесу они не нужны): - их параллельно собирает пул процессов (`gpu_record_worker`), иначе один - поток Python душит GPU (тот простаивал на ~2/3 времени поколения). + их параллельно собирает пул процессов (`gpu_record_worker`). + executor — пул для ПАРАЛЛЕЛЬНОГО предвычисления физики ступеней + (precompute_worker): без него один поток Python готовит физику всей + популяции и GPU простаивает, ожидая данные. """ from gausse.optim.objective import EvaluationResult, compute_cost_rub - states = [] - for g in genomes: - cfg, ix, iv = decode(g, db, bounds) - centers = [0.0] - for gap in cfg.inter_stage_gaps_m: - centers.append(centers[-1] + gap) - states.append(_State(g, cfg, ix, iv, centers, cur_x=ix, cur_v=iv)) + if executor is not None: + pres = list(executor.map(precompute_worker, genomes, chunksize=64)) + else: + pres = [precompute_genome(g, db, bounds) for g in genomes] + states = [ + _State(g, pre, cur_x=pre.initial_x_m, cur_v=pre.initial_v_mps) + for g, pre in zip(genomes, pres) + ] _simulate_states(xp, states) results = [] for st in states: - cost = compute_cost_rub(st.config, db) if build_details else None + if build_details: + cfg, _, _ = decode(st.genome, db, bounds) + cost = compute_cost_rub(cfg, db) ok = st.alive and st.energy_in_j > 0 detail = build_detail( - st.config, _FEASIBLE if ok else None, db, st.initial_x_m, st.initial_v_mps, + cfg, _FEASIBLE if ok else None, db, st.pre.initial_x_m, st.pre.initial_v_mps, st.genome.tube_inner_d_m, st.genome.tube_wall_m, ) if build_details else None + if not build_details: + cost = None if ok: eff = st.kinetic_delta_j / st.energy_in_j results.append(EvaluationResult( @@ -199,7 +302,7 @@ def evaluate_genomes_gpu( detail=detail, )) else: - n_stages = len(st.config.stages) + n_stages = len(st.pre.stages) progress = (st.failed_stage_index or 0) / n_stages if n_stages else 0.0 results.append(EvaluationResult( feasible=False, cost_rub=cost, fitness=-1.0 + progress, @@ -269,11 +372,8 @@ def run_gpu_sweep( states = [] for _ in range(n): g = sample_genome(db, bounds, rng) - cfg, ix, iv = decode(g, db, bounds) - centers = [0.0] - for gap in cfg.inter_stage_gaps_m: - centers.append(centers[-1] + gap) - states.append(_State(g, cfg, ix, iv, centers, cur_x=ix, cur_v=iv)) + pre = precompute_genome(g, db, bounds) + states.append(_State(g, pre, cur_x=pre.initial_x_m, cur_v=pre.initial_v_mps)) _simulate_states(xp, states) @@ -281,10 +381,10 @@ def run_gpu_sweep( for st in states: if st.alive and st.energy_in_j > 0: eff = st.kinetic_delta_j / st.energy_in_j - records.append(_record(st, db, backend, eff, st.cur_v)) + records.append(_record(st, db, bounds, backend, eff, st.cur_v)) n_feasible += 1 else: - records.append(_record(st, db, backend, None, None)) + records.append(_record(st, db, bounds, backend, None, None)) insert_runs(conn, records) for r in records: logger.update(r.feasible, r.efficiency) @@ -295,10 +395,11 @@ def run_gpu_sweep( return {"n_runs": done, "n_feasible": n_feasible, "backend": backend} -def _record(st: _State, db, backend, efficiency, exit_v): +def _record(st: _State, db, bounds: SearchBounds, backend, efficiency, exit_v): feasible = efficiency is not None + cfg, _, _ = decode(st.genome, db, bounds) detail = build_detail( - st.config, _FEASIBLE if feasible else None, db, st.initial_x_m, st.initial_v_mps, + cfg, _FEASIBLE if feasible else None, db, st.pre.initial_x_m, st.pre.initial_v_mps, st.genome.tube_inner_d_m, st.genome.tube_wall_m, ) return RunRecord( diff --git a/src/gausse/optim/evolutionary.py b/src/gausse/optim/evolutionary.py index 5913f05..a889da2 100644 --- a/src/gausse/optim/evolutionary.py +++ b/src/gausse/optim/evolutionary.py @@ -153,7 +153,9 @@ def run_evolution( if gpu_xp is not None: from gausse.gpu.batch_sweep import evaluate_genomes_gpu, gpu_record_worker - results = evaluate_genomes_gpu(gpu_xp, population, db, bounds, build_details=False) + results = evaluate_genomes_gpu( + gpu_xp, population, db, bounds, build_details=False, executor=executor + ) pairs = list(zip(population, results)) payload = [ (g, {