From 863060ee78956bbd2a9137b503547fdf5af54a74 Mon Sep 17 00:00:00 2001 From: Codex Agent Date: Sat, 22 Nov 2025 18:04:36 +0100 Subject: [PATCH] Improve reactor controls, reactions, and maintenance UI --- AGENTS.md | 1 + artifacts/last_state.json | 104 ++++++++++++++++++++++++++++++++++ src/reactor_sim/atomic.py | 66 +++++++++++++++++++++ src/reactor_sim/constants.py | 8 ++- src/reactor_sim/control.py | 3 +- src/reactor_sim/dashboard.py | 38 +++++++++++++ src/reactor_sim/fuel.py | 50 ++++++++++++++++ src/reactor_sim/neutronics.py | 14 +++-- src/reactor_sim/reactor.py | 30 +++++++++- src/reactor_sim/simulation.py | 24 ++++++++ src/reactor_sim/state.py | 17 +++++- src/reactor_sim/thermal.py | 8 ++- src/reactor_sim/turbine.py | 4 ++ tests/test_atomic.py | 54 ++++++++++++++++++ tests/test_fuel.py | 40 +++++++++++++ 15 files changed, 447 insertions(+), 14 deletions(-) create mode 100644 artifacts/last_state.json create mode 100644 tests/test_atomic.py create mode 100644 tests/test_fuel.py diff --git a/AGENTS.md b/AGENTS.md index c66e388..c558f7b 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -23,6 +23,7 @@ Pytest is the preferred framework. Place tests under `tests/` mirroring the `rea ## Commit & Pull Request Guidelines Write commits in imperative mood (`Add turbine moisture separator model`). Squash small WIP commits before review. Pull requests must describe the scenario, attach log excerpts (`snapshots.json` diff) or plots for novel behavior, and link to any tracking issues. Include validation notes: commands run, expected trends (e.g., outlet temperature increase), and outstanding risks. +- Agents should automatically create a commit immediately after code/config/text changes are finished; use a clear imperative subject and avoid piling multiple logical changes into one commit unless the task explicitly says otherwise. ## Safety & Configuration Sim parameters live in constructors; never hard-code environment-specific paths. When adding new physics, guard unstable calculations with clamps and highlight operating limits in comments. Sensitive experiments (e.g., beyond design basis) should default to disabled flags so scripted studies remain safe by default. diff --git a/artifacts/last_state.json b/artifacts/last_state.json new file mode 100644 index 0000000..ab97f0b --- /dev/null +++ b/artifacts/last_state.json @@ -0,0 +1,104 @@ +{ + "control": { + "setpoint_mw": 3000.0, + "rod_fraction": 0.3872181667930225 + }, + "plant": { + "core": { + "fuel_temperature": 633.2342060602825, + "neutron_flux": 31335131.43776027, + "reactivity_margin": 0.006241667151538859, + "power_output_mw": 2976.900110634399, + "burnup": 0.004945521266135161 + }, + "primary_loop": { + "temperature_in": 295.0, + "temperature_out": 338.7522062115579, + "pressure": 14.0, + "mass_flow_rate": 16200.0, + "steam_quality": 0.0 + }, + "secondary_loop": { + "temperature_in": 295.0, + "temperature_out": 360.05377003828596, + "pressure": 6.650537700382859, + "mass_flow_rate": 10880.0, + "steam_quality": 0.6505377003828593 + }, + "turbines": [ + { + "steam_enthalpy": 3090.3226202297155, + "shaft_power_mw": 336.9056685758782, + "electrical_output_mw": 323.4294418328431, + "condenser_temperature": 305.0, + "load_demand_mw": 0.0, + "load_supplied_mw": 0.0 + }, + { + "steam_enthalpy": 3090.3226202297155, + "shaft_power_mw": 336.9056685758782, + "electrical_output_mw": 323.4294418328431, + "condenser_temperature": 305.0, + "load_demand_mw": 0.0, + "load_supplied_mw": 0.0 + }, + { + "steam_enthalpy": 3090.3226202297155, + "shaft_power_mw": 336.9056685758782, + "electrical_output_mw": 323.4294418328431, + "condenser_temperature": 305.0, + "load_demand_mw": 0.0, + "load_supplied_mw": 0.0 + } + ], + "time_elapsed": 600.0 + }, + "metadata": { + "primary_pump_active": true, + "secondary_pump_active": true, + "turbine_active": true, + "turbine_units": [ + true, + true, + true + ], + "shutdown": false, + "consumer": { + "online": false, + "demand_mw": 800.0, + "name": "Grid" + } + }, + "health": { + "core": { + "name": "core", + "integrity": 0.9400000000000066, + "failed": false + }, + "primary_pump": { + "name": "primary_pump", + "integrity": 0.5800000000000063, + "failed": false + }, + "secondary_pump": { + "name": "secondary_pump", + "integrity": 0.1120000000000021, + "failed": false + }, + "turbine_1": { + "name": "turbine_1", + "integrity": 0.610000000000003, + "failed": false + }, + "turbine_2": { + "name": "turbine_2", + "integrity": 0.610000000000003, + "failed": false + }, + "turbine_3": { + "name": "turbine_3", + "integrity": 0.610000000000003, + "failed": false + } + } +} \ No newline at end of file diff --git a/src/reactor_sim/atomic.py b/src/reactor_sim/atomic.py index 6ad0908..fb95afa 100644 --- a/src/reactor_sim/atomic.py +++ b/src/reactor_sim/atomic.py @@ -164,6 +164,14 @@ class FissionEvent: energy_mev: float +@dataclass(frozen=True) +class ReactionEvent: + parent: Atom + products: tuple[Atom, ...] + emitted_particles: tuple[str, ...] + energy_mev: float + + def make_atom(protons: int, neutrons: int) -> Atom: return Atom(symbol=element_symbol(protons), protons=protons, neutrons=neutrons) @@ -235,6 +243,64 @@ class AtomicPhysics: ) return event + def alpha_decay(self, atom: Atom) -> ReactionEvent: + if atom.protons <= 2 or atom.neutrons <= 2: + return ReactionEvent(parent=atom, products=(atom,), emitted_particles=(), energy_mev=0.0) + daughter = make_atom(atom.protons - 2, atom.neutrons - 2) + alpha = make_atom(2, 2) + energy = 5.0 # Typical alpha decay energy + return ReactionEvent(parent=atom, products=(daughter, alpha), emitted_particles=("alpha",), energy_mev=energy) + + def beta_minus_decay(self, atom: Atom) -> ReactionEvent: + if atom.neutrons <= 0: + return ReactionEvent(parent=atom, products=(atom,), emitted_particles=(), energy_mev=0.0) + daughter = make_atom(atom.protons + 1, atom.neutrons - 1) + energy = 1.0 + return ReactionEvent( + parent=atom, + products=(daughter,), + emitted_particles=("beta-", "anti-nu_e"), + energy_mev=energy, + ) + + def beta_plus_decay(self, atom: Atom) -> ReactionEvent: + if atom.protons <= 1: + return ReactionEvent(parent=atom, products=(atom,), emitted_particles=(), energy_mev=0.0) + daughter = make_atom(atom.protons - 1, atom.neutrons + 1) + energy = 1.0 + return ReactionEvent( + parent=atom, + products=(daughter,), + emitted_particles=("beta+", "nu_e"), + energy_mev=energy, + ) + + def electron_capture(self, atom: Atom) -> ReactionEvent: + if atom.protons <= 1: + return ReactionEvent(parent=atom, products=(atom,), emitted_particles=(), energy_mev=0.0) + daughter = make_atom(atom.protons - 1, atom.neutrons + 1) + energy = 0.5 + return ReactionEvent( + parent=atom, + products=(daughter,), + emitted_particles=("nu_e", "gamma"), + energy_mev=energy, + ) + + def gamma_emission(self, atom: Atom, energy_mev: float = 0.2) -> ReactionEvent: + return ReactionEvent(parent=atom, products=(atom,), emitted_particles=("gamma",), energy_mev=energy_mev) + + def spontaneous_fission(self, atom: Atom) -> ReactionEvent: + # Reuse the generic splitter to keep mass/charge balanced. + split = self._generic_split(atom) + neutrons = ("n",) * max(0, split.emitted_neutrons) + return ReactionEvent( + parent=atom, + products=split.products, + emitted_particles=neutrons, + energy_mev=split.energy_mev, + ) + def _generic_split(self, atom: Atom) -> FissionEvent: heavy_mass = max(1, math.floor(atom.mass_number * 0.55)) light_mass = atom.mass_number - heavy_mass diff --git a/src/reactor_sim/constants.py b/src/reactor_sim/constants.py index 111c581..ac53b76 100644 --- a/src/reactor_sim/constants.py +++ b/src/reactor_sim/constants.py @@ -10,10 +10,16 @@ COOLANT_HEAT_CAPACITY = 4_200.0 # J/(kg*K) for water/steam COOLANT_DENSITY = 700.0 # kg/m^3 averaged between phases MAX_CORE_TEMPERATURE = 1_800.0 # K MAX_PRESSURE = 15.0 # MPa typical PWR primary loop limit -CONTROL_ROD_SPEED = 0.05 # fraction insertion per second +CONTROL_ROD_SPEED = 0.03 # fraction insertion per second STEAM_TURBINE_EFFICIENCY = 0.34 GENERATOR_EFFICIENCY = 0.96 ENVIRONMENT_TEMPERATURE = 295.0 # K AMU_TO_KG = 1.660_539_066_60e-27 MEV_TO_J = 1.602_176_634e-13 ELECTRON_FISSION_CROSS_SECTION = 5e-16 # cm^2, tuned for simulation scale +# Threshold inventories (event counts) for flagging common poisons in diagnostics. +KEY_POISON_THRESHOLDS = { + "Xe": 1e20, # xenon + "I": 1e20, # iodine precursor to xenon + "Sm": 5e19, # samarium +} diff --git a/src/reactor_sim/control.py b/src/reactor_sim/control.py index 7ea5d96..536eed6 100644 --- a/src/reactor_sim/control.py +++ b/src/reactor_sim/control.py @@ -27,7 +27,8 @@ class ControlSystem: if self.manual_control: return self.rod_fraction error = (state.power_output_mw - self.setpoint_mw) / self.setpoint_mw - adjustment = -error * 0.3 + # When power is low (negative error) withdraw rods; when high, insert them. + adjustment = error * 0.2 adjustment = clamp(adjustment, -constants.CONTROL_ROD_SPEED * dt, constants.CONTROL_ROD_SPEED * dt) previous = self.rod_fraction self.rod_fraction = clamp(self.rod_fraction + adjustment, 0.0, 0.95) diff --git a/src/reactor_sim/dashboard.py b/src/reactor_sim/dashboard.py index 227292e..6df9428 100644 --- a/src/reactor_sim/dashboard.py +++ b/src/reactor_sim/dashboard.py @@ -9,6 +9,7 @@ from dataclasses import dataclass from pathlib import Path from typing import Optional +from . import constants from .commands import ReactorCommand from .reactor import Reactor from .simulation import ReactorSimulation @@ -50,8 +51,12 @@ class ReactorDashboard: DashboardKey("o", "Toggle secondary pump"), DashboardKey("t", "Toggle turbine"), DashboardKey("1/2/3", "Toggle turbine units 1-3"), + DashboardKey("y/u/i", "Maintain turbine 1/2/3"), DashboardKey("c", "Toggle consumer"), DashboardKey("r", "Reset & clear state"), + DashboardKey("m", "Maintain primary pump"), + DashboardKey("n", "Maintain secondary pump"), + DashboardKey("k", "Maintain core (requires shutdown)"), DashboardKey("+/-", "Withdraw/insert rods"), DashboardKey("[/]", "Adjust consumer demand −/+50 MW"), DashboardKey("s", "Setpoint −250 MW"), @@ -142,6 +147,18 @@ class ReactorDashboard: self._queue_command(ReactorCommand(power_setpoint=self.reactor.control.setpoint_mw + 250.0)) elif ch in (ord("a"), ord("A")): self._queue_command(ReactorCommand(rod_manual=not self.reactor.control.manual_control)) + elif ch in (ord("m"), ord("M")): + self._queue_command(ReactorCommand.maintain("primary_pump")) + elif ch in (ord("n"), ord("N")): + self._queue_command(ReactorCommand.maintain("secondary_pump")) + elif ch in (ord("k"), ord("K")): + self._queue_command(ReactorCommand.maintain("core")) + elif ch in (ord("y"), ord("Y")): + self._queue_command(ReactorCommand.maintain("turbine_1")) + elif ch in (ord("u"), ord("U")): + self._queue_command(ReactorCommand.maintain("turbine_2")) + elif ch in (ord("i"), ord("I")): + self._queue_command(ReactorCommand.maintain("turbine_3")) def _queue_command(self, command: ReactorCommand) -> None: if self.pending_command is None: @@ -240,10 +257,12 @@ class ReactorDashboard: ("Core Power", f"{state.core.power_output_mw:8.1f} MW"), ("Neutron Flux", f"{state.core.neutron_flux:10.2e}"), ("Rods", f"{self.reactor.control.rod_fraction:.3f}"), + ("Rod Mode", "AUTO" if not self.reactor.control.manual_control else "MANUAL"), ("Setpoint", f"{self.reactor.control.setpoint_mw:7.0f} MW"), ("Reactivity", f"{state.core.reactivity_margin:+.4f}"), ], ) + y = self._draw_section(win, y, "Key Poisons / Emitters", self._poison_lines(state)) y = self._draw_section( win, y, @@ -300,6 +319,7 @@ class ReactorDashboard: "Start pumps before withdrawing rods.", "Bring turbine and consumer online after thermal stabilization.", "Toggle turbine units (1/2/3) for staggered maintenance.", + "Use m/n/k/y/u/i to request maintenance (stop equipment first).", "Press 'r' to reset/clear state if you want a cold start.", "Watch component health to avoid automatic trips.", ] @@ -370,6 +390,24 @@ class ReactorDashboard: def _total_load_demand(self, state: PlantState) -> float: return sum(t.load_demand_mw for t in state.turbines) + def _poison_lines(self, state: PlantState) -> list[tuple[str, str]]: + inventory = state.core.fission_product_inventory or {} + particles = state.core.emitted_particles or {} + lines: list[tuple[str, str]] = [] + def fmt(symbol: str, label: str) -> tuple[str, str]: + qty = inventory.get(symbol, 0.0) + threshold = constants.KEY_POISON_THRESHOLDS.get(symbol) + flag = " !" if threshold is not None and qty >= threshold else "" + return (f"{label}{flag}", f"{qty:9.2e}") + + lines.append(fmt("Xe", "Xe (xenon)")) + lines.append(fmt("Sm", "Sm (samarium)")) + lines.append(fmt("I", "I (iodine)")) + lines.append(("Neutrons (src)", f"{particles.get('n', 0.0):9.2e}")) + lines.append(("Gammas", f"{particles.get('gamma', 0.0):9.2e}")) + lines.append(("Alphas", f"{particles.get('alpha', 0.0):9.2e}")) + return lines + def _draw_health_bar(self, win: "curses._CursesWindow", start_y: int) -> None: height, width = win.getmaxyx() if start_y >= height - 2: diff --git a/src/reactor_sim/fuel.py b/src/reactor_sim/fuel.py index eb3f0ea..02ab0f5 100644 --- a/src/reactor_sim/fuel.py +++ b/src/reactor_sim/fuel.py @@ -29,6 +29,7 @@ class FuelAssembly: mass_kg: float fissile_atom: Atom = field(default_factory=lambda: make_atom(92, 143)) atomic_physics: AtomicPhysics = field(default_factory=AtomicPhysics) + decay_activity_base: float = 1e16 # nominal decay events per second at burnup 0 def available_energy_j(self, state: CoreState) -> float: fraction_remaining = max(0.0, 1.0 - state.burnup) @@ -57,3 +58,52 @@ class FuelAssembly: power_mw, ) return max(0.0, power_mw), event_rate, event + + def decay_reaction_effects( + self, state: CoreState + ) -> tuple[float, float, dict[str, float], dict[str, float]]: + """Return (power MW, neutron source rate, product rates, particle rates).""" + # Scale activity with burnup (fission products) and a small flux contribution. + activity = self.decay_activity_base * (0.05 + state.burnup) * (1.0 + 0.00001 * state.neutron_flux) + activity = max(0.0, min(activity, 1e20)) + + reactions = [ + self.atomic_physics.alpha_decay(self.fissile_atom), + self.atomic_physics.beta_minus_decay(self.fissile_atom), + self.atomic_physics.beta_plus_decay(self.fissile_atom), + self.atomic_physics.electron_capture(self.fissile_atom), + self.atomic_physics.gamma_emission(self.fissile_atom), + self.atomic_physics.spontaneous_fission(self.fissile_atom), + ] + weights = [1.0, 1.2, 0.5, 0.4, 1.6, 0.05] + total_weight = sum(weights) + + power_watts = 0.0 + neutron_source = 0.0 + product_rates: dict[str, float] = {} + particle_rates: dict[str, float] = {} + + for event, weight in zip(reactions, weights): + rate = activity * (weight / total_weight) + power_watts += rate * event.energy_mev * constants.MEV_TO_J + for atom in event.products: + product_rates[atom.symbol] = product_rates.get(atom.symbol, 0.0) + rate + for particle in event.emitted_particles: + particle_rates[particle] = particle_rates.get(particle, 0.0) + rate + if particle == "n": + neutron_source += rate + + power_mw = power_watts / constants.MEGAWATT + capped_power = min(power_mw, 0.1 * max(state.power_output_mw, 1.0)) + LOGGER.debug( + "Decay reactions contribute %.2f MW (raw %.2f MW) with neutron source %.2e 1/s", + capped_power, + power_mw, + neutron_source, + ) + return max(0.0, capped_power), neutron_source, product_rates, particle_rates + + def decay_reaction_power(self, state: CoreState) -> float: + """Compatibility shim: return only power contribution.""" + power_mw, _, _, _ = self.decay_reaction_effects(state) + return power_mw diff --git a/src/reactor_sim/neutronics.py b/src/reactor_sim/neutronics.py index 01692ae..9f64121 100644 --- a/src/reactor_sim/neutronics.py +++ b/src/reactor_sim/neutronics.py @@ -15,18 +15,19 @@ LOGGER = logging.getLogger(__name__) def temperature_feedback(temp: float) -> float: """Negative coefficient: higher temperature lowers reactivity.""" reference = 900.0 - coefficient = -2.5e-5 + coefficient = -1.5e-5 return coefficient * (temp - reference) def xenon_poisoning(flux: float) -> float: - return min(0.015, 2e-9 * flux) + return min(0.01, 5e-10 * flux) @dataclass class NeutronDynamics: beta_effective: float = 0.0065 delayed_neutron_fraction: float = 0.0008 + external_source_coupling: float = 1e-6 def reactivity(self, state: CoreState, control_fraction: float) -> float: rho = ( @@ -37,14 +38,15 @@ class NeutronDynamics: ) return rho - def flux_derivative(self, state: CoreState, rho: float) -> float: + def flux_derivative(self, state: CoreState, rho: float, external_source_rate: float = 0.0) -> float: generation_time = constants.NEUTRON_LIFETIME beta = self.beta_effective - return ((rho - beta) / generation_time) * state.neutron_flux + 1e5 + source_term = self.external_source_coupling * external_source_rate + return ((rho - beta) / generation_time) * state.neutron_flux + 1e5 + source_term - def step(self, state: CoreState, control_fraction: float, dt: float) -> None: + def step(self, state: CoreState, control_fraction: float, dt: float, external_source_rate: float = 0.0) -> None: rho = self.reactivity(state, control_fraction) - d_flux = self.flux_derivative(state, rho) + d_flux = self.flux_derivative(state, rho, external_source_rate) state.neutron_flux = max(0.0, state.neutron_flux + d_flux * dt) state.reactivity_margin = rho LOGGER.debug( diff --git a/src/reactor_sim/reactor.py b/src/reactor_sim/reactor.py index 837b650..c273e52 100644 --- a/src/reactor_sim/reactor.py +++ b/src/reactor_sim/reactor.py @@ -39,6 +39,7 @@ class Reactor: turbine_active: bool = True turbine_unit_active: list[bool] = field(default_factory=lambda: [True, True, True]) shutdown: bool = False + poison_alerts: set[str] = field(default_factory=set) def __post_init__(self) -> None: if not self.turbines: @@ -72,6 +73,8 @@ class Reactor: reactivity_margin=-0.02, power_output_mw=0.1, burnup=0.0, + fission_product_inventory={}, + emitted_particles={}, ) primary = CoolantLoopState( temperature_in=ambient, @@ -111,15 +114,28 @@ class Reactor: overrides = self._apply_command(command, state) rod_fraction = overrides.get("rod_fraction", rod_fraction) - self.neutronics.step(state.core, rod_fraction, dt) + decay_power, decay_neutron_source, decay_products, decay_particles = self.fuel.decay_reaction_effects( + state.core + ) + self.neutronics.step(state.core, rod_fraction, dt, external_source_rate=decay_neutron_source) prompt_power, fission_rate, fission_event = self.fuel.prompt_energy_rate( state.core.neutron_flux, rod_fraction ) - decay_power = decay_heat_fraction(state.core.burnup) * state.core.power_output_mw - total_power = prompt_power + decay_power + decay_heat = decay_heat_fraction(state.core.burnup) * state.core.power_output_mw + total_power = prompt_power + decay_heat + decay_power state.core.power_output_mw = total_power state.core.update_burnup(dt) + # Track fission products and emitted particles for diagnostics. + products: dict[str, float] = {} + for atom in fission_event.products: + products[atom.symbol] = products.get(atom.symbol, 0.0) + fission_rate * dt + for k, v in decay_products.items(): + products[k] = products.get(k, 0.0) + v * dt + particles: dict[str, float] = {k: v * dt for k, v in decay_particles.items()} + state.core.add_products(products) + state.core.add_emitted_particles(particles) + self._check_poison_alerts(state) pump_demand = overrides.get("coolant_demand", self.control.coolant_demand(state.primary_loop)) if self.primary_pump_active: @@ -371,3 +387,11 @@ class Reactor: ) ) return plant + + def _check_poison_alerts(self, state: PlantState) -> None: + inventory = state.core.fission_product_inventory or {} + for symbol, threshold in constants.KEY_POISON_THRESHOLDS.items(): + amount = inventory.get(symbol, 0.0) + if amount >= threshold and symbol not in self.poison_alerts: + self.poison_alerts.add(symbol) + LOGGER.warning("Poison level high: %s inventory %.2e exceeds %.2e", symbol, amount, threshold) diff --git a/src/reactor_sim/simulation.py b/src/reactor_sim/simulation.py index ac52d83..b289b3a 100644 --- a/src/reactor_sim/simulation.py +++ b/src/reactor_sim/simulation.py @@ -27,6 +27,29 @@ def _default_state_path() -> Path: return Path(custom) if custom else Path("artifacts/last_state.json") +def _poison_log_path() -> Path: + custom = os.getenv("FISSION_POISON_LOG") + return Path(custom) if custom else Path("artifacts/poisons.json") + + +def _write_poison_log(path: Path, state: PlantState | None, snapshots: list[dict[str, float]] | None = None) -> None: + if state is None and snapshots: + latest = snapshots[-1] + inventory = latest.get("products", {}) + particles = latest.get("particles", {}) + time_elapsed = latest.get("time_elapsed", 0.0) + elif state is not None: + inventory = state.core.fission_product_inventory + particles = state.core.emitted_particles + time_elapsed = state.time_elapsed + else: + return + path.parent.mkdir(parents=True, exist_ok=True) + payload = {"time_elapsed": time_elapsed, "products": inventory, "particles": particles} + path.write_text(json.dumps(payload, indent=2)) + LOGGER.info("Wrote poison snapshot to %s", path) + + @dataclass class ReactorSimulation: reactor: Reactor @@ -113,6 +136,7 @@ def main() -> None: snapshots = sim.log() LOGGER.info("Captured %d snapshots", len(snapshots)) print(json.dumps(snapshots[-5:], indent=2)) + _write_poison_log(_poison_log_path(), sim.last_state, snapshots) except KeyboardInterrupt: sim.stop() LOGGER.warning("Simulation interrupted by user") diff --git a/src/reactor_sim/state.py b/src/reactor_sim/state.py index b180588..a9aba1a 100644 --- a/src/reactor_sim/state.py +++ b/src/reactor_sim/state.py @@ -16,11 +16,21 @@ class CoreState: reactivity_margin: float # delta rho power_output_mw: float # MW thermal burnup: float # fraction of fuel consumed + fission_product_inventory: dict[str, float] = field(default_factory=dict) + emitted_particles: dict[str, float] = field(default_factory=dict) def update_burnup(self, dt: float) -> None: produced_energy_mwh = self.power_output_mw * (dt / 3600.0) self.burnup = clamp(self.burnup + produced_energy_mwh * 1e-5, 0.0, 0.99) + def add_products(self, products: dict[str, float]) -> None: + for element, amount in products.items(): + self.fission_product_inventory[element] = self.fission_product_inventory.get(element, 0.0) + amount + + def add_emitted_particles(self, particles: dict[str, float]) -> None: + for name, amount in particles.items(): + self.emitted_particles[name] = self.emitted_particles.get(name, 0.0) + amount + @dataclass class CoolantLoopState: @@ -61,6 +71,8 @@ class PlantState: "primary_outlet_temp": self.primary_loop.temperature_out, "secondary_pressure": self.secondary_loop.pressure, "turbine_electric": self.total_electrical_output(), + "products": self.core.fission_product_inventory, + "particles": self.core.emitted_particles, } def total_electrical_output(self) -> float: @@ -71,6 +83,9 @@ class PlantState: @classmethod def from_dict(cls, data: dict) -> "PlantState": + core_blob = dict(data["core"]) + inventory = core_blob.pop("fission_product_inventory", {}) + particles = core_blob.pop("emitted_particles", {}) turbines_blob = data.get("turbines") if turbines_blob is None: # Compatibility with previous single-turbine snapshots. @@ -78,7 +93,7 @@ class PlantState: turbines_blob = [old_turbine] if old_turbine else [] turbines = [TurbineState(**t) for t in turbines_blob] return cls( - core=CoreState(**data["core"]), + core=CoreState(**core_blob, fission_product_inventory=inventory, emitted_particles=particles), primary_loop=CoolantLoopState(**data["primary_loop"]), secondary_loop=CoolantLoopState(**data["secondary_loop"]), turbines=turbines, diff --git a/src/reactor_sim/thermal.py b/src/reactor_sim/thermal.py index 5212e1c..bc7a305 100644 --- a/src/reactor_sim/thermal.py +++ b/src/reactor_sim/thermal.py @@ -36,8 +36,12 @@ class ThermalSolver: def step_core(self, core: CoreState, primary: CoolantLoopState, power_mw: float, dt: float) -> None: temp_rise = temperature_rise(power_mw, primary.mass_flow_rate) primary.temperature_out = primary.temperature_in + temp_rise - core.fuel_temperature += 0.1 * (power_mw - temp_rise) * dt - core.fuel_temperature = min(core.fuel_temperature, constants.MAX_CORE_TEMPERATURE) + # Fuel heats from any power not immediately convected away, and cools toward the primary outlet. + heating = 0.005 * max(0.0, power_mw - temp_rise) * dt + cooling = 0.05 * max(0.0, core.fuel_temperature - primary.temperature_out) * dt + core.fuel_temperature += heating - cooling + # Keep fuel temperature bounded and never below the coolant outlet temperature. + core.fuel_temperature = min(max(primary.temperature_out, core.fuel_temperature), constants.MAX_CORE_TEMPERATURE) LOGGER.debug( "Primary loop: flow=%.0f kg/s temp_out=%.1fK core_temp=%.1fK", primary.mass_flow_rate, diff --git a/src/reactor_sim/turbine.py b/src/reactor_sim/turbine.py index 2b19b23..266e77d 100644 --- a/src/reactor_sim/turbine.py +++ b/src/reactor_sim/turbine.py @@ -25,6 +25,7 @@ class SteamGenerator: class Turbine: generator_efficiency: float = constants.GENERATOR_EFFICIENCY mechanical_efficiency: float = constants.STEAM_TURBINE_EFFICIENCY + rated_output_mw: float = 400.0 # cap per unit electrical output def step( self, @@ -37,6 +38,9 @@ class Turbine: available_power = max(steam_power_mw, (enthalpy * mass_flow / 1_000.0) / 1_000.0) shaft_power_mw = available_power * self.mechanical_efficiency electrical = shaft_power_mw * self.generator_efficiency + if electrical > self.rated_output_mw: + electrical = self.rated_output_mw + shaft_power_mw = electrical / max(1e-6, self.generator_efficiency) condenser_temp = max(305.0, loop.temperature_in - 20.0) state.steam_enthalpy = enthalpy state.shaft_power_mw = shaft_power_mw diff --git a/tests/test_atomic.py b/tests/test_atomic.py new file mode 100644 index 0000000..6cb2509 --- /dev/null +++ b/tests/test_atomic.py @@ -0,0 +1,54 @@ +from reactor_sim.atomic import AtomicPhysics, make_atom + + +def test_alpha_decay_reduces_mass_and_charge(): + physics = AtomicPhysics(seed=1) + parent = make_atom(92, 143) + event = physics.alpha_decay(parent) + daughter, alpha = event.products + assert daughter.protons == 90 + assert daughter.neutrons == 141 + assert alpha.protons == 2 and alpha.neutrons == 2 + assert "alpha" in event.emitted_particles + assert event.energy_mev > 0 + + +def test_beta_minus_conserves_mass_number(): + physics = AtomicPhysics(seed=2) + parent = make_atom(26, 30) + event = physics.beta_minus_decay(parent) + (daughter,) = event.products + assert daughter.mass_number == parent.mass_number + assert daughter.protons == parent.protons + 1 + assert "beta-" in event.emitted_particles + + +def test_beta_plus_and_electron_capture_lower_proton_count(): + physics = AtomicPhysics(seed=3) + parent = make_atom(15, 16) + beta_plus = physics.beta_plus_decay(parent) + capture = physics.electron_capture(parent) + assert beta_plus.products[0].protons == parent.protons - 1 + assert capture.products[0].protons == parent.protons - 1 + assert "beta+" in beta_plus.emitted_particles + assert "gamma" in capture.emitted_particles + + +def test_gamma_emission_keeps_nuclide_same(): + physics = AtomicPhysics(seed=4) + parent = make_atom(8, 8) + event = physics.gamma_emission(parent, energy_mev=0.3) + (product,) = event.products + assert product == parent + assert "gamma" in event.emitted_particles + assert event.energy_mev == 0.3 + + +def test_spontaneous_fission_creates_two_fragments(): + physics = AtomicPhysics(seed=5) + parent = make_atom(92, 143) + event = physics.spontaneous_fission(parent) + assert len(event.products) == 2 + total_mass = sum(atom.mass_number for atom in event.products) + assert total_mass <= parent.mass_number + assert event.energy_mev > 0 diff --git a/tests/test_fuel.py b/tests/test_fuel.py new file mode 100644 index 0000000..49ccd6a --- /dev/null +++ b/tests/test_fuel.py @@ -0,0 +1,40 @@ +from reactor_sim.atomic import AtomicPhysics +from reactor_sim.fuel import FuelAssembly +from reactor_sim.state import CoreState + + +def _core_state( + temp: float = 600.0, + flux: float = 1e6, + burnup: float = 0.05, + power: float = 100.0, +) -> CoreState: + return CoreState( + fuel_temperature=temp, + neutron_flux=flux, + reactivity_margin=0.0, + power_output_mw=power, + burnup=burnup, + ) + + +def test_decay_reaction_power_positive_and_capped(): + physics = AtomicPhysics(seed=7) + fuel = FuelAssembly(enrichment=0.045, mass_kg=80_000.0, atomic_physics=physics) + state = _core_state() + power = fuel.decay_reaction_power(state) + assert power > 0 + # Should stay under 10% of current power_output_mw. + assert power <= state.power_output_mw * 0.1 + 1e-6 + + +def test_decay_reaction_power_scales_with_burnup_and_flux(): + physics = AtomicPhysics(seed=9) + fuel = FuelAssembly(enrichment=0.045, mass_kg=80_000.0, atomic_physics=physics) + low = _core_state(burnup=0.01, flux=1e5, power=50.0) + high = _core_state(burnup=0.5, flux=5e6, power=200.0) + low_power = fuel.decay_reaction_power(low) + high_power = fuel.decay_reaction_power(high) + assert high_power > low_power + # Still capped to a reasonable fraction of the prevailing power. + assert high_power <= high.power_output_mw * 0.1 + 1e-6