From 2046dcba10747504f1006aff1268c3e86f02e961 Mon Sep 17 00:00:00 2001 From: jze9 Date: Mon, 6 Jul 2026 19:45:32 +0500 Subject: [PATCH] Add component database layer and validated physics core - components/schema.py + database.py: typed loader for real retail parts (wire, capacitors, switches, sensors, projectile materials) - physics/inductance.py: multilayer solenoid inductance (Wheeler) with ferromagnetic-slug coupling via smooth overlap model and analytic dL/dx - physics/force.py: F = 0.5*I^2*dL/dx with a saturation clamp - physics/circuit.py: coupled [Q,I,x,v] discharge ODE - Validated against textbook RLC analytical solutions (under/over/critically damped) and energy conservation, not just spot-checked by eye - PLAN.md updated with animation (per-run field visualization) and GPU-accelerated batch sweep as explicit later stages Co-Authored-By: Claude Sonnet 5 --- .gitignore | 1 + PLAN.md | 2 + src/gausse/components/database.py | 61 +++++++++ src/gausse/components/schema.py | 61 +++++++++ src/gausse/physics/circuit.py | 66 ++++++++++ src/gausse/physics/constants.py | 6 + src/gausse/physics/force.py | 36 ++++++ src/gausse/physics/inductance.py | 113 +++++++++++++++++ tests/fixtures/data/capacitors.json | 11 ++ tests/fixtures/data/projectile_materials.json | 10 ++ tests/fixtures/data/sensors.json | 18 +++ tests/fixtures/data/switches.json | 13 ++ tests/fixtures/data/wires.json | 22 ++++ tests/test_circuit_rlc_analytical.py | 116 ++++++++++++++++++ tests/test_database_loading.py | 45 +++++++ tests/test_force_model.py | 35 ++++++ tests/test_inductance.py | 82 +++++++++++++ 17 files changed, 698 insertions(+) create mode 100644 src/gausse/components/database.py create mode 100644 src/gausse/components/schema.py create mode 100644 src/gausse/physics/circuit.py create mode 100644 src/gausse/physics/constants.py create mode 100644 src/gausse/physics/force.py create mode 100644 src/gausse/physics/inductance.py create mode 100644 tests/fixtures/data/capacitors.json create mode 100644 tests/fixtures/data/projectile_materials.json create mode 100644 tests/fixtures/data/sensors.json create mode 100644 tests/fixtures/data/switches.json create mode 100644 tests/fixtures/data/wires.json create mode 100644 tests/test_circuit_rlc_analytical.py create mode 100644 tests/test_database_loading.py create mode 100644 tests/test_force_model.py create mode 100644 tests/test_inductance.py diff --git a/.gitignore b/.gitignore index 55e7099..dd4f57d 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,4 @@ dist/ *.sqlite3-shm out/ .DS_Store +.claude/settings.local.json diff --git a/PLAN.md b/PLAN.md index a2fe970..b937303 100644 --- a/PLAN.md +++ b/PLAN.md @@ -31,4 +31,6 @@ ChipDip, Cable.ru и др.). - [ ] **Этап 5 — Хранилище результатов**: `storage/schema.py` + `database.py` — SQLite (WAL), таблица `runs` со всеми прогонами (успех/провал + честная причина), однопроцессный писатель поверх многопроцессной очереди. - [ ] **Этап 6 — Поиск и оптимизация**: `optim/search_space.py` (геном переменной длины), `objective.py` (КПД), дешёвый квазистатический предфильтр, `optim/sweep.py` (Monte Carlo/LHS, миллионы прогонов), `optim/evolutionary.py` (ГА + coordinate-descent полировка) — всё пишется в общую таблицу `runs`. - [ ] **Этап 7 — Отчётность**: `report/plots.py`, `summary.py`, `bom.py`, раздел "ограничения модели", `cli.py` (`gausse sweep/evolve/simulate/report`). + - [ ] `report/animate.py` — анимация одного прогона: положение снаряда в трубе по времени + визуализация поля/тока каждой катушки (свечение/интенсивность цвета ~ ток), сохранение в GIF/MP4 (matplotlib `FuncAnimation`). Нужна по запросу пользователя — "графика где будет показана симуляция пролёта цилиндра по трубе и электромагнитные поля в каждый момент времени". - [ ] **Этап 8 — Сквозная проверка**: резюмируемый прогон на реальной базе компонентов, проверка честной записи в SQLite, финальный отчёт (КПД, скорость, стоимость, сравнение датчиков) с разделом ограничений. +- [ ] **Этап 9 — GPU-ускорение массового sweep (сервер с GTX 1070)**: после того как CPU/`scipy.solve_ivp`-модель провалидирована тестами (Этап 2-4) — батч-версия интегратора с ФИКСИРОВАННЫМ шагом (RK4/полу-неявная схема), считающая сразу N траекторий параллельно как один тензор (`cupy`, если доступна CUDA, иначе векторизованный `numpy`/`numba`), для прогона по-настоящему миллионов конфигураций на сервере. Важно: сначала корректность на CPU, потом скорость на GPU — численные результаты GPU-пути должны сверяться с CPU-эталоном на контрольной выборке, чтобы ускорение не подменило точность честными числами "для галочки". diff --git a/src/gausse/components/database.py b/src/gausse/components/database.py new file mode 100644 index 0000000..85dcef4 --- /dev/null +++ b/src/gausse/components/database.py @@ -0,0 +1,61 @@ +"""Загрузка базы реальных компонентов из JSON и фильтры для поиска.""" + +import json +from dataclasses import dataclass, field +from pathlib import Path + +from gausse.components.schema import ( + CapacitorSpec, + ProjectileMaterialSpec, + SensorSpec, + SwitchSpec, + WireSpec, +) + +DEFAULT_DATA_DIR = Path(__file__).parent / "data" + + +def _load_json(path: Path) -> list[dict]: + if not path.exists(): + return [] + return json.loads(path.read_text(encoding="utf-8")) + + +@dataclass +class ComponentDatabase: + wires: list[WireSpec] = field(default_factory=list) + capacitors: list[CapacitorSpec] = field(default_factory=list) + switches: list[SwitchSpec] = field(default_factory=list) + sensors: list[SensorSpec] = field(default_factory=list) + projectile_materials: list[ProjectileMaterialSpec] = field(default_factory=list) + + @classmethod + def load(cls, data_dir: Path = DEFAULT_DATA_DIR) -> "ComponentDatabase": + return cls( + wires=[WireSpec(**row) for row in _load_json(data_dir / "wires.json")], + capacitors=[ + CapacitorSpec(**row) for row in _load_json(data_dir / "capacitors.json") + ], + switches=[SwitchSpec(**row) for row in _load_json(data_dir / "switches.json")], + sensors=[SensorSpec(**row) for row in _load_json(data_dir / "sensors.json")], + projectile_materials=[ + ProjectileMaterialSpec(**row) + for row in _load_json(data_dir / "projectile_materials.json") + ], + ) + + def wires_by_material(self, material: str) -> list[WireSpec]: + return [w for w in self.wires if w.material == material] + + def capacitors_rated_at_least(self, voltage_v: float) -> list[CapacitorSpec]: + return [c for c in self.capacitors if c.voltage_v >= voltage_v] + + def switches_rated_for(self, current_a: float, voltage_v: float) -> list[SwitchSpec]: + return [ + s + for s in self.switches + if s.max_current_a >= current_a and s.max_voltage_v >= voltage_v + ] + + def sensors_by_kind(self, kind: str) -> list[SensorSpec]: + return [s for s in self.sensors if s.kind == kind] diff --git a/src/gausse/components/schema.py b/src/gausse/components/schema.py new file mode 100644 index 0000000..51b663e --- /dev/null +++ b/src/gausse/components/schema.py @@ -0,0 +1,61 @@ +"""Типизированные описания реальных, продающихся в розницу компонентов.""" + +from dataclasses import dataclass +from typing import Literal + + +@dataclass(frozen=True) +class WireSpec: + part_id: str + material: Literal["copper", "aluminum"] + gauge_mm: float + insulation_od_mm: float + resistivity_ohm_m: float + max_current_a: float + price_per_m: float + source: str + + +@dataclass(frozen=True) +class CapacitorSpec: + part_number: str + capacitance_uf: float + voltage_v: float + esr_ohm: float + max_current_a: float + price: float + source: str + + +@dataclass(frozen=True) +class SwitchSpec: + part_number: str + kind: Literal["SCR", "MOSFET", "IGBT"] + max_current_a: float + max_voltage_v: float + on_resistance_ohm: float | None + on_voltage_drop_v: float | None + turn_on_time_ns: float + price: float + source: str + + +@dataclass(frozen=True) +class SensorSpec: + part_number: str + kind: Literal["inductive", "optical", "hall"] + propagation_delay_ns: float + price: float + source: str + sensitivity_v_per_mps: float | None = None + threshold_v: float | None = None + + +@dataclass(frozen=True) +class ProjectileMaterialSpec: + name: str + density_kg_m3: float + mu_r: float + b_sat_tesla: float + price_per_kg: float + source: str diff --git a/src/gausse/physics/circuit.py b/src/gausse/physics/circuit.py new file mode 100644 index 0000000..461f58c --- /dev/null +++ b/src/gausse/physics/circuit.py @@ -0,0 +1,66 @@ +"""ОДУ разряда конденсатора через ключ и катушку, связанное с механикой снаряда. + +Состояние: [Q, I, x, v] — заряд конденсатора, ток в катушке, положение и +скорость снаряда. Индуктивность зависит от x (проход снаряда), поэтому +уравнение на dI/dt выведено из закона сохранения потокосцепления +Ψ = L(x)*I: + + V_c = Q / C + dQ/dt = -I + dI/dt = (V_c - I*R_total - I*(dL/dx)*v) / L(x) + dv/dt = F(x, I) / m + dx/dt = v +""" + +from dataclasses import dataclass + +import numpy as np + +from gausse.physics.force import force_on_slug_newtons +from gausse.physics.inductance import CoilInductanceModel + +Q, I, X, V = range(4) + + +@dataclass(frozen=True) +class StageCircuitParams: + capacitance_f: float + r_total_ohm: float + mass_kg: float + mu_eff: float + total_turns: int + coil_length_m: float + b_sat_tesla: float + + +def derivatives( + t: float, + state: np.ndarray, + inductance_model: CoilInductanceModel, + params: StageCircuitParams, +) -> list[float]: + q, current, x, v = state + l_x = inductance_model.l_of_x(x) + dl_dx = inductance_model.dl_dx(x) + + v_c = q / params.capacitance_f + d_q = -current + d_i = (v_c - current * params.r_total_ohm - current * dl_dx * v) / l_x + force = force_on_slug_newtons( + current_a=current, + dl_dx_unsaturated=dl_dx, + mu_eff=params.mu_eff, + total_turns=params.total_turns, + coil_length_m=params.coil_length_m, + b_sat_tesla=params.b_sat_tesla, + ) + d_v = force / params.mass_kg + return [d_q, d_i, v, d_v] + + +def zero_current_crossing_event(t: float, state: np.ndarray, *_args) -> float: + return state[I] + + +zero_current_crossing_event.terminal = True +zero_current_crossing_event.direction = -1.0 diff --git a/src/gausse/physics/constants.py b/src/gausse/physics/constants.py new file mode 100644 index 0000000..bce4acc --- /dev/null +++ b/src/gausse/physics/constants.py @@ -0,0 +1,6 @@ +"""Физические константы и справочные удельные сопротивления.""" + +MU_0 = 4e-7 * 3.141592653589793 # Гн/м, магнитная проницаемость вакуума + +RESISTIVITY_COPPER_OHM_M = 1.68e-8 # при ~20°C +RESISTIVITY_ALUMINUM_OHM_M = 2.82e-8 # при ~20°C diff --git a/src/gausse/physics/force.py b/src/gausse/physics/force.py new file mode 100644 index 0000000..34506f2 --- /dev/null +++ b/src/gausse/physics/force.py @@ -0,0 +1,36 @@ +"""Сила, действующая на ферромагнитный снаряд: F = 0.5*I^2*dL/dx с насыщением. + +Насыщение — грубая, явно приближённая поправка (v1): оцениваем поле внутри +катушки простой соленоидной формулой и, если оно превышает B_sat материала +снаряда, пропорционально ослабляем вклад сердечника в dL/dx. Реальное +насыщение нелинейно и зависит от геометрии — здесь используется линейный +клэмп как консервативная эвристика, а не точная кривая намагничивания. +""" + +from gausse.physics.constants import MU_0 + + +def solenoid_field_estimate_tesla( + mu_eff: float, total_turns: int, coil_length_m: float, current_a: float +) -> float: + turns_per_m = total_turns / coil_length_m + return MU_0 * mu_eff * turns_per_m * abs(current_a) + + +def saturation_scale(b_estimate_tesla: float, b_sat_tesla: float) -> float: + if b_estimate_tesla <= 0: + return 1.0 + return min(1.0, b_sat_tesla / b_estimate_tesla) + + +def force_on_slug_newtons( + current_a: float, + dl_dx_unsaturated: float, + mu_eff: float, + total_turns: int, + coil_length_m: float, + b_sat_tesla: float, +) -> float: + b_estimate = solenoid_field_estimate_tesla(mu_eff, total_turns, coil_length_m, current_a) + scale = saturation_scale(b_estimate, b_sat_tesla) + return 0.5 * current_a**2 * dl_dx_unsaturated * scale diff --git a/src/gausse/physics/inductance.py b/src/gausse/physics/inductance.py new file mode 100644 index 0000000..f74854c --- /dev/null +++ b/src/gausse/physics/inductance.py @@ -0,0 +1,113 @@ +"""Индуктивность многослойного соленоида с ферромагнитным сердечником. + +Приближённая инженерная модель (v1): воздушная индуктивность по формуле +Уилера для многослойных катушек, скорректированная на присутствие +ферромагнитного снаряда через коэффициент размагничивания конечного +цилиндра (аппроксимация вытянутым сфероидом). Точная картина требует МКЭ +(FEA) — это явно вынесено как ограничение модели, см. PLAN.md. +""" + +import math +from dataclasses import dataclass + +import numpy as np + +INCH_M = 0.0254 + + +@dataclass(frozen=True) +class WindingGeometry: + coil_length_m: float + radial_depth_m: float + mean_radius_m: float + total_turns: int + total_wire_length_m: float + + +def winding_geometry( + tube_od_m: float, wire_od_m: float, turns_per_layer: int, layers: int +) -> WindingGeometry: + """Геометрия обмотки из параметров намотки, а не свободной "длины катушки".""" + coil_length_m = turns_per_layer * wire_od_m + radial_depth_m = layers * wire_od_m + inner_radius_m = tube_od_m / 2 + total_wire_length_m = 0.0 + for i in range(layers): + layer_radius_m = inner_radius_m + (i + 0.5) * wire_od_m + total_wire_length_m += 2 * math.pi * layer_radius_m * turns_per_layer + mean_radius_m = inner_radius_m + radial_depth_m / 2 + return WindingGeometry( + coil_length_m=coil_length_m, + radial_depth_m=radial_depth_m, + mean_radius_m=mean_radius_m, + total_turns=turns_per_layer * layers, + total_wire_length_m=total_wire_length_m, + ) + + +def air_core_inductance_wheeler(mean_radius_m: float, length_m: float, depth_m: float, turns: int) -> float: + """Формула Уилера для многослойного воздушного соленоида (Гн). + + Точность ~1% при a, l, c одного порядка величины (Wheeler, 1928). + """ + a_in = mean_radius_m / INCH_M + l_in = length_m / INCH_M + c_in = depth_m / INCH_M + l_uh = (0.8 * a_in**2 * turns**2) / (6 * a_in + 9 * l_in + 10 * c_in) + return l_uh * 1e-6 + + +def demagnetizing_factor_prolate(aspect_ratio: float) -> float: + """Коэффициент размагничивания вдоль длинной оси вытянутого сфероида. + + aspect_ratio = длина / диаметр снаряда. При aspect_ratio -> 1 (сфера) + возвращает 1/3; при aspect_ratio -> inf стремится к 0. + """ + m = max(aspect_ratio, 1.0 + 1e-6) + root = math.sqrt(m**2 - 1) + return (1 / (m**2 - 1)) * ((m / (2 * root)) * math.log((m + root) / (m - root)) - 1) + + +def effective_permeability(mu_r: float, demagnetizing_factor: float) -> float: + return mu_r / (1 + demagnetizing_factor * (mu_r - 1)) + + +def _sigmoid(z: np.ndarray) -> np.ndarray: + return 0.5 * (1 + np.tanh(z / 2)) + + +@dataclass(frozen=True) +class CoilInductanceModel: + """L(x) и аналитическая dL/dx для катушки с проходящим ферромагнитным снарядом. + + x — положение центра снаряда относительно центра катушки (м), катушка + считается центрированной в x=0 в своей локальной системе координат. + """ + + l_air_h: float + coil_length_m: float + slug_length_m: float + mu_eff: float + smoothing_width_m: float + + def _edges(self) -> tuple[float, float]: + half_span = (self.coil_length_m + self.slug_length_m) / 2 + return -half_span, half_span + + def overlap_fraction(self, x: np.ndarray) -> np.ndarray: + e1, e2 = self._edges() + w = self.smoothing_width_m + return _sigmoid((x - e1) / w) * _sigmoid((e2 - x) / w) + + def d_overlap_dx(self, x: np.ndarray) -> np.ndarray: + e1, e2 = self._edges() + w = self.smoothing_width_m + s1 = _sigmoid((x - e1) / w) + s2 = _sigmoid((e2 - x) / w) + return (s1 * s2 / w) * (s2 - s1) + + def l_of_x(self, x: np.ndarray) -> np.ndarray: + return self.l_air_h * (1 + (self.mu_eff - 1) * self.overlap_fraction(x)) + + def dl_dx(self, x: np.ndarray) -> np.ndarray: + return self.l_air_h * (self.mu_eff - 1) * self.d_overlap_dx(x) diff --git a/tests/fixtures/data/capacitors.json b/tests/fixtures/data/capacitors.json new file mode 100644 index 0000000..a638aeb --- /dev/null +++ b/tests/fixtures/data/capacitors.json @@ -0,0 +1,11 @@ +[ + { + "part_number": "test-cap-450v-220uf", + "capacitance_uf": 220.0, + "voltage_v": 450.0, + "esr_ohm": 0.05, + "max_current_a": 200.0, + "price": 150.0, + "source": "fixture" + } +] diff --git a/tests/fixtures/data/projectile_materials.json b/tests/fixtures/data/projectile_materials.json new file mode 100644 index 0000000..dcd26e8 --- /dev/null +++ b/tests/fixtures/data/projectile_materials.json @@ -0,0 +1,10 @@ +[ + { + "name": "test-steel", + "density_kg_m3": 7850.0, + "mu_r": 200.0, + "b_sat_tesla": 1.8, + "price_per_kg": 100.0, + "source": "fixture" + } +] diff --git a/tests/fixtures/data/sensors.json b/tests/fixtures/data/sensors.json new file mode 100644 index 0000000..aa0f084 --- /dev/null +++ b/tests/fixtures/data/sensors.json @@ -0,0 +1,18 @@ +[ + { + "part_number": "test-optical", + "kind": "optical", + "propagation_delay_ns": 500.0, + "price": 20.0, + "source": "fixture" + }, + { + "part_number": "test-inductive-pickup", + "kind": "inductive", + "propagation_delay_ns": 0.0, + "price": 5.0, + "source": "fixture", + "sensitivity_v_per_mps": 0.02, + "threshold_v": 0.1 + } +] diff --git a/tests/fixtures/data/switches.json b/tests/fixtures/data/switches.json new file mode 100644 index 0000000..ef1e974 --- /dev/null +++ b/tests/fixtures/data/switches.json @@ -0,0 +1,13 @@ +[ + { + "part_number": "test-scr-25a", + "kind": "SCR", + "max_current_a": 25.0, + "max_voltage_v": 600.0, + "on_resistance_ohm": null, + "on_voltage_drop_v": 1.5, + "turn_on_time_ns": 2000.0, + "price": 80.0, + "source": "fixture" + } +] diff --git a/tests/fixtures/data/wires.json b/tests/fixtures/data/wires.json new file mode 100644 index 0000000..89999db --- /dev/null +++ b/tests/fixtures/data/wires.json @@ -0,0 +1,22 @@ +[ + { + "part_id": "test-cu-0.5", + "material": "copper", + "gauge_mm": 0.5, + "insulation_od_mm": 0.55, + "resistivity_ohm_m": 1.72e-8, + "max_current_a": 3.0, + "price_per_m": 5.0, + "source": "fixture" + }, + { + "part_id": "test-al-0.5", + "material": "aluminum", + "gauge_mm": 0.5, + "insulation_od_mm": 0.55, + "resistivity_ohm_m": 2.82e-8, + "max_current_a": 2.2, + "price_per_m": 2.0, + "source": "fixture" + } +] diff --git a/tests/test_circuit_rlc_analytical.py b/tests/test_circuit_rlc_analytical.py new file mode 100644 index 0000000..612b450 --- /dev/null +++ b/tests/test_circuit_rlc_analytical.py @@ -0,0 +1,116 @@ +"""Проверка ОДУ разряда против аналитического решения простого RLC-контура. + +dL/dx заморожена в ноль (mu_eff=1.0), поэтому механика отвязана от контура +и [Q, I] эволюционируют как классический последовательный RLC-разряд +конденсатора. Сравниваем численное решение solve_ivp с учебным аналитическим +решением во всех трёх режимах затухания. +""" + +import math + +import numpy as np +import pytest +from scipy.integrate import solve_ivp + +from gausse.physics.circuit import Q, StageCircuitParams, derivatives +from gausse.physics.inductance import CoilInductanceModel + +L_H = 100e-6 +C_F = 200e-6 +V0 = 400.0 +Q0 = C_F * V0 + + +def _make_frozen_model() -> CoilInductanceModel: + return CoilInductanceModel( + l_air_h=L_H, coil_length_m=0.05, slug_length_m=0.03, mu_eff=1.0, smoothing_width_m=0.001 + ) + + +def _integrate(r_ohm: float, t_span: tuple[float, float], t_eval: np.ndarray): + model = _make_frozen_model() + params = StageCircuitParams( + capacitance_f=C_F, + r_total_ohm=r_ohm, + mass_kg=1.0, + mu_eff=1.0, + total_turns=100, + coil_length_m=0.05, + b_sat_tesla=1e9, # эффективно отключает клэмп насыщения для этого теста + ) + state0 = [Q0, 0.0, 0.0, 0.0] + return solve_ivp( + derivatives, + t_span, + state0, + args=(model, params), + t_eval=t_eval, + method="RK45", + rtol=1e-10, + atol=1e-14, + ) + + +def test_underdamped_matches_analytical_solution(): + r_ohm = 0.05 + alpha = r_ohm / (2 * L_H) + omega0 = 1 / math.sqrt(L_H * C_F) + omega_d = math.sqrt(omega0**2 - alpha**2) + + t_eval = np.linspace(0, 5 / alpha, 200) + sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval) + + expected_q = Q0 * np.exp(-alpha * t_eval) * ( + np.cos(omega_d * t_eval) + (alpha / omega_d) * np.sin(omega_d * t_eval) + ) + expected_i = Q0 * (omega0**2 / omega_d) * np.exp(-alpha * t_eval) * np.sin(omega_d * t_eval) + + assert np.allclose(sol.y[Q], expected_q, rtol=1e-3, atol=Q0 * 1e-4) + assert np.allclose(sol.y[1], expected_i, rtol=1e-3, atol=abs(expected_i).max() * 1e-3) + + +def test_overdamped_matches_analytical_solution(): + r_ohm = 5.0 + alpha = r_ohm / (2 * L_H) + omega0 = 1 / math.sqrt(L_H * C_F) + root = math.sqrt(alpha**2 - omega0**2) + s1, s2 = -alpha + root, -alpha - root + + t_eval = np.linspace(0, 5 / abs(s1), 200) + sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval) + + a_coef = -Q0 * s2 / (s1 - s2) + b_coef = Q0 * s1 / (s1 - s2) + expected_q = a_coef * np.exp(s1 * t_eval) + b_coef * np.exp(s2 * t_eval) + expected_i = -(a_coef * s1 * np.exp(s1 * t_eval) + b_coef * s2 * np.exp(s2 * t_eval)) + + assert np.allclose(sol.y[Q], expected_q, rtol=1e-3, atol=Q0 * 1e-4) + assert np.allclose(sol.y[1], expected_i, rtol=1e-3, atol=max(abs(expected_i).max(), 1e-9) * 1e-3) + + +def test_critically_damped_matches_analytical_solution(): + alpha = 1 / math.sqrt(L_H * C_F) # R = 2*sqrt(L/C) => alpha = omega0 + r_ohm = 2 * alpha * L_H + + t_eval = np.linspace(0, 5 / alpha, 200) + sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval) + + expected_q = Q0 * (1 + alpha * t_eval) * np.exp(-alpha * t_eval) + expected_i = alpha**2 * Q0 * t_eval * np.exp(-alpha * t_eval) + + assert np.allclose(sol.y[Q], expected_q, rtol=1e-3, atol=Q0 * 1e-4) + assert np.allclose(sol.y[1], expected_i, rtol=1e-3, atol=max(abs(expected_i).max(), 1e-9) * 1e-3) + + +def test_energy_conserved_when_lossless(): + """С R=0 энергия конденсатора должна переходить в энергию катушки без потерь.""" + r_ohm = 0.0 + t_eval = np.linspace(0, 1e-3, 500) + sol = _integrate(r_ohm, (t_eval[0], t_eval[-1]), t_eval) + + cap_energy = sol.y[Q] ** 2 / (2 * C_F) + inductor_energy = 0.5 * L_H * sol.y[1] ** 2 + total_energy = cap_energy + inductor_energy + initial_energy = Q0**2 / (2 * C_F) + + assert np.allclose(total_energy, initial_energy, rtol=1e-4) diff --git a/tests/test_database_loading.py b/tests/test_database_loading.py new file mode 100644 index 0000000..e3553c6 --- /dev/null +++ b/tests/test_database_loading.py @@ -0,0 +1,45 @@ +from pathlib import Path + +from gausse.components.database import ComponentDatabase + +FIXTURES = Path(__file__).parent / "fixtures" / "data" + + +def test_loads_all_component_types(): + db = ComponentDatabase.load(FIXTURES) + assert len(db.wires) == 2 + assert len(db.capacitors) == 1 + assert len(db.switches) == 1 + assert len(db.sensors) == 2 + assert len(db.projectile_materials) == 1 + + +def test_wires_by_material_filters_correctly(): + db = ComponentDatabase.load(FIXTURES) + copper = db.wires_by_material("copper") + assert len(copper) == 1 + assert copper[0].part_id == "test-cu-0.5" + + +def test_capacitors_rated_at_least_filters_by_voltage(): + db = ComponentDatabase.load(FIXTURES) + assert len(db.capacitors_rated_at_least(400.0)) == 1 + assert len(db.capacitors_rated_at_least(500.0)) == 0 + + +def test_switches_rated_for_current_and_voltage(): + db = ComponentDatabase.load(FIXTURES) + assert len(db.switches_rated_for(20.0, 500.0)) == 1 + assert len(db.switches_rated_for(30.0, 500.0)) == 0 + + +def test_sensors_by_kind(): + db = ComponentDatabase.load(FIXTURES) + assert len(db.sensors_by_kind("optical")) == 1 + assert len(db.sensors_by_kind("inductive")) == 1 + + +def test_missing_data_dir_returns_empty_database(tmp_path): + db = ComponentDatabase.load(tmp_path) + assert db.wires == [] + assert db.capacitors == [] diff --git a/tests/test_force_model.py b/tests/test_force_model.py new file mode 100644 index 0000000..52a372d --- /dev/null +++ b/tests/test_force_model.py @@ -0,0 +1,35 @@ +import pytest + +from gausse.physics.force import ( + force_on_slug_newtons, + saturation_scale, + solenoid_field_estimate_tesla, +) + + +def test_field_estimate_scales_with_current_and_turns(): + b_low = solenoid_field_estimate_tesla(mu_eff=100, total_turns=100, coil_length_m=0.05, current_a=10) + b_high = solenoid_field_estimate_tesla(mu_eff=100, total_turns=100, coil_length_m=0.05, current_a=100) + assert b_high > b_low + + +def test_saturation_scale_is_one_below_bsat(): + assert saturation_scale(b_estimate_tesla=0.5, b_sat_tesla=1.8) == pytest.approx(1.0) + + +def test_saturation_scale_clamps_above_bsat(): + scale = saturation_scale(b_estimate_tesla=3.6, b_sat_tesla=1.8) + assert scale == pytest.approx(0.5) + + +def test_force_is_reduced_once_saturated(): + common = dict(dl_dx_unsaturated=1e-3, mu_eff=200, total_turns=300, coil_length_m=0.05, b_sat_tesla=1.8) + force_unsaturated = force_on_slug_newtons(current_a=5, **common) + force_saturated = force_on_slug_newtons(current_a=500, **common) + # без клэмпа сила росла бы как I^2 (в 10000 раз); с насыщением рост должен быть намного меньше + assert force_saturated / force_unsaturated < 5000 + + +def test_force_zero_at_zero_current(): + common = dict(dl_dx_unsaturated=1e-3, mu_eff=200, total_turns=300, coil_length_m=0.05, b_sat_tesla=1.8) + assert force_on_slug_newtons(current_a=0, **common) == pytest.approx(0.0) diff --git a/tests/test_inductance.py b/tests/test_inductance.py new file mode 100644 index 0000000..259905f --- /dev/null +++ b/tests/test_inductance.py @@ -0,0 +1,82 @@ +import math + +import numpy as np +import pytest + +from gausse.physics.inductance import ( + CoilInductanceModel, + air_core_inductance_wheeler, + demagnetizing_factor_prolate, + effective_permeability, + winding_geometry, +) + + +def test_winding_geometry_scales_with_turns_and_layers(): + geo = winding_geometry(tube_od_m=0.02, wire_od_m=0.001, turns_per_layer=50, layers=3) + assert geo.total_turns == 150 + assert geo.coil_length_m == pytest.approx(0.05) + assert geo.radial_depth_m == pytest.approx(0.003) + assert geo.total_wire_length_m > 0 + # больше слоёв -> больше провода при тех же витках на слой + geo_more_layers = winding_geometry(0.02, 0.001, 50, 5) + assert geo_more_layers.total_wire_length_m > geo.total_wire_length_m + + +def test_wheeler_inductance_is_positive_and_grows_with_turns(): + l_small = air_core_inductance_wheeler(mean_radius_m=0.015, length_m=0.05, depth_m=0.003, turns=50) + l_large = air_core_inductance_wheeler(mean_radius_m=0.015, length_m=0.05, depth_m=0.003, turns=200) + assert l_small > 0 + assert l_large > l_small + + +def test_demagnetizing_factor_bounds(): + # почти сфера (aspect_ratio -> 1): классическое значение 1/3 + assert demagnetizing_factor_prolate(1.001) == pytest.approx(1 / 3, abs=0.02) + # очень вытянутый стержень: коэффициент размагничивания стремится к 0 + assert demagnetizing_factor_prolate(200.0) < 0.01 + + +def test_effective_permeability_limits(): + # без размагничивания (Nd=0) эффективная проницаемость равна мю материала + assert effective_permeability(mu_r=500.0, demagnetizing_factor=0.0) == pytest.approx(500.0) + # полное размагничивание (Nd=1) гасит эффект проницаемости до 1 + assert effective_permeability(mu_r=500.0, demagnetizing_factor=1.0) == pytest.approx(1.0) + + +@pytest.fixture +def model() -> CoilInductanceModel: + return CoilInductanceModel( + l_air_h=100e-6, + coil_length_m=0.05, + slug_length_m=0.03, + mu_eff=50.0, + smoothing_width_m=0.001, + ) + + +def test_inductance_peaks_when_slug_centered(model): + l_center = model.l_of_x(np.array(0.0)) + l_far = model.l_of_x(np.array(0.2)) + assert l_center > l_far + assert l_far == pytest.approx(model.l_air_h, rel=1e-3) + + +def test_dl_dx_vanishes_at_symmetric_center(model): + assert model.dl_dx(np.array(0.0)) == pytest.approx(0.0, abs=1e-9) + + +def test_dl_dx_sign_flips_across_center(model): + # слева от центра индуктивность растёт (dL/dx > 0), справа убывает (< 0) + assert model.dl_dx(np.array(-0.01)) > 0 + assert model.dl_dx(np.array(0.01)) < 0 + + +def test_analytic_dl_dx_matches_finite_difference(model): + xs = np.linspace(-0.06, 0.06, 25) + h = 1e-7 + analytic = np.array([model.dl_dx(x) for x in xs]) + finite_diff = np.array( + [(model.l_of_x(x + h) - model.l_of_x(x - h)) / (2 * h) for x in xs] + ) + assert np.allclose(analytic, finite_diff, rtol=1e-3, atol=1e-9)