Add cooldown telemetry, RHR cooldown path, and shutdown regression test

AGENTS.md

@@ -2,6 +2,7 @@
 
 ## Project Structure & Module Organization
 All source code lives under `src/reactor_sim`. Submodules map to plant systems: `fuel.py` and `neutronics.py` govern fission power, `thermal.py` and `coolant.py` cover heat transfer and pumps, `turbine.py` drives the steam cycle, and `consumer.py` represents external electrical loads. High-level coordination happens in `reactor.py` and `simulation.py`. The convenience runner `run_simulation.py` executes the default scenario; add notebooks or scenario scripts under `experiments/` (create as needed). Keep assets such as plots or exported telemetry inside `artifacts/`.
+Feature catalog lives in `FEATURES.md`; update it whenever behavior is added or changed. Long-term realism tasks are tracked in `TODO.md`; keep it in sync as work progresses.
 
 ## Build, Test, and Development Commands
 - `.venv/bin/python -m reactor_sim.simulation` — run the default 10-minute transient and print JSON snapshots using the checked-in virtualenv.
@@ -30,3 +31,6 @@ Sim parameters live in constructors; never hard-code environment-specific paths.
 
 ## Reliability & Persistence
 Component wear is tracked via `failures.py`; stress from overheating, pump starvation, or turbine imbalance will degrade integrity and eventually disable the affected subsystem with automatic SCRAM for core damage. Plant state now persists between runs to `artifacts/last_state.json` by default (override via `FISSION_STATE_PATH` or the explicit save/load env vars). In the dashboard, press `r` to clear the saved snapshot and reboot the reactor to a cold, green-field state whenever you need a clean slate.
+
+## Session Context
+See `CONTEXT_NOTES.md` for the latest behavioral changes, controls, and assumptions carried over between sessions.

CONTEXT_NOTES.md

@@ -0,0 +1,24 @@
+# Session Context Notes
+
+- Reactor model: two-loop but tuned to RBMK-like pressures (nominal ~7 MPa for both loops). Loop pressure clamps to saturation baseline when pumps are off; pumps ramp flow/pressure over spool time when stopping or starting.
+- Feedwater & level: secondary steam drum uses shrink/swell-aware level sensing to drive a feedwater valve; makeup flow scales with steam draw toward the target level instead of instant inventory clamps.
+- Chemistry/fouling: plant tracks dissolved O2/boron/sodium; impurities plus high temp/steam draw increase HX fouling (reduces UA) and add condenser fouling/back-pressure. Oxygen degasses with steam; impurity ingress accelerates when venting.
+- Reactivity bias: boron ppm now biases shutdown reactivity; a very slow trim nudges boron toward the power setpoint after ~300s of operation.
+- Turbines: produce zero output unless steam quality is present and effective steam flow is >10 kg/s. Turbine panel shows steam availability (enthalpy × quality × mass flow) and steam enthalpy instead of loop pressure; condenser pressure/temperature/fouling shown with nominal bounds.
+- Generators: two diesel units, rated 50 MW, spool time 10s. Auto mode default `False`; manual toggles b/v. Auto stops when no load. Relief valves toggles l (primary) / ; (secondary) and displayed per loop.
+- Pumps: per-unit controls g/h (primary), j/k (secondary). Flow/pressure ramp down over spool when pumps stop. Pump status thresholds use >0.1 kg/s to show STOPPING.
+- Maintenance hotkeys: p (core, requires shutdown), m/n (primary 1/2), ,/. (secondary 1/2), B/V (generator 1/2), y/u/i (turbine 1/2/3).
+- Dashboard: two-column layout, trends panel for fuel temp and core power (delta and rate). Power Stats show aux demand/supply, generator and turbine output. Turbine panel shows steam availability/enthalpy instead of loop pressure. Core/temp/power lines include nominal/max.
+- Thermal updates: primary/secondary inlet temps now back-computed; when secondary flow is near zero, loops cool toward ambient over time. Relief venting removes mass/enthalpy with a multi-second ramp toward ~1 MPa and quenches superheat toward target-pressure saturation. Pumps cap target pressure when reliefs are open. Condenser modeled with vacuum pump drawdown, cooling-sink temperature, and fouling-driven back-pressure penalty.
+- Meltdown threshold: 2873 K. Auto rod control clears shutdown when set to auto and adjusts rods. Control rod worth/tuning currently unchanged.
+- Tests: `pytest` passing after all changes. Key regression additions include generator manual mode, turbine no-steam output, auto rod control, and passive cool-down.
+- Coolant demand fixed: demand increases when primary outlet is above target (sign was flipped before), so hot loops ramp flow instead of backing off.
+- Pump spin-down: pressure/flow ramp down over pump spool time with STOPPING->OFF threshold at 0.1 kg/s; prevents instant drop when last pump stops.
+- Steam pressure display shows 0 unless steam quality ≥0.05 and flow ≥100 kg/s to avoid showing pump head as steam pressure.
+- Passive cool-down: when secondary flow ~0, loops cool toward ambient; primary inlet/outlet back-propagated from transferred heat and ambient.
+- Relief valves: l (primary) and ; (secondary) vent with mass/enthalpy loss, ramp pressure toward ~1 MPa over several seconds, cap pump targets while open; status displayed per loop.
+- Generator behavior: starting/running only produce power when load is present; auto off by default; manual toggles b/v; auto stops with no load; base aux drops to 0 when idle/cold.
+- Pressure tying: loop pressure floors to saturation(temp) when pumps off; pump targets aim for ~7 MPa nominal RBMK-like setpoints.
+- Turbines/governor: require meaningful steam flow/quality; otherwise zero output. Steam supply pressure in turbine panel reads 0 when no steam. Throttle now biases toward load demand with a governor term, and overspeed/overload >105% rated electrical trips a unit and spools it down.
+- Rod control now supports three banks with weighted worth; xenon/iodine tracked with decay and burn-out; new UA·ΔT_lm steam-generator model and pump head/flow curves.
+- Dashboard shows heat-exchanger ΔT/efficiency and protections; pumps and HX changes documented in FEATURES.md / TODO.md.

FEATURES.md

@@ -0,0 +1,10 @@
+## C.O.R.E. feature set
+
+- **Core physics**: point-kinetics with per-bank delayed neutron precursors, temperature feedback, fuel burnup penalty, xenon/iodine buildup with decay and burn-out, and rod-bank worth curves.
+- **Rod control**: three rod banks with weighted worth; auto controller chases 3 GW setpoint with safety backoff and filtered power feedback; manual mode with staged bank motion and SCRAM; state persists across runs. Soluble boron bias contributes slow negative reactivity and trims toward the setpoint.
+- **Coolant & hydraulics**: primary/secondary pumps with head/flow curves, power draw scaling, wear tracking; pressure floors tied to saturation; auxiliary power model with generator auto-start.
+- **Heat transfer**: steam-generator UA·ΔT_lm model with a pinch cap to keep the primary outlet hotter than the secondary, coolant heating uses total fission power with fuel heating decoupled from exchanger draw, and the secondary thermal solver includes passive cool-down plus steam-drum mass/energy balance with latent heat and a shrink/swell-aware feedwater valve controller; dissolved oxygen/sodium drive HX fouling that reduces effective UA.
+- **Pressurizer & inventory**: primary pressurizer trims toward 7 MPa with level tracking, loop inventories/levels steer flow availability, secondary steam boil-off draws down level with auto makeup, and pumps reduce flow/status to `CAV` when NPSH is insufficient.
+- **Steam cycle**: three turbines with spool dynamics, throttle mapping and a simple governor with overspeed/overload trip, condenser vacuum/back-pressure with fouling and cooling sink temperature, chemistry-driven fouling/back-pressure penalties, load dispatch to consumer, steam quality gating for output, generator states with batteries/spool, and steam enthalpy-driven availability readout on the dashboard.
+- **Protections & failures**: health monitor degrading components under stress, automatic SCRAM on core or heat-sink loss plus DNB/subcool and secondary level/pressure trips, relief valves per loop with venting/mass loss and pump pressure caps, maintenance actions to restore integrity.
+- **Persistence & ops**: snapshots auto-save/load to `artifacts/last_state.json`; dashboard with live metrics, protections/warnings, heat-exchanger telemetry, component health, and control shortcuts.

TODO.md

@@ -0,0 +1,24 @@
+## Future realism upgrades
+
+- [x] Steam generator UA·ΔT_lm heat exchange and pump head/flow curves; keep validating temps under nominal load.
+- [x] Rod banks with worth curves, xenon/samarium buildup, and delayed-group kinetics per bank.
+- [x] Pressurizer behavior, primary/secondary inventory and level effects, and pump NPSH/cavitation checks.
+- [x] Model feedwater/steam-drum mass-energy balance, turbine throttle/efficiency maps, and condenser back-pressure.
+- [x] Introduce CHF/DNB margin, clad/fuel split temps, and SCRAM matrix for subcooling loss or SG level/pressure trips.
+- [x] Flesh out condenser behavior: vacuum pump limits, cooling water temperature coupling, and dynamic back-pressure with fouling.
+- [x] Dashboard polish: compact turbine/generator rows, color critical warnings (SCRAM/heat-sink), and reduce repeated log noise.
+- [ ] Dashboard multi-page view (F1/F2): retain numeric view on F1; future F2 schematic should mirror real PWR layout with ASCII art, flow/relief status, and minimal animations; add help/status hints and size checks; keep perf sane.
+- [x] Core thermal realism: add a simple radial fuel model (pellet/rim/clad temps) with burnup-driven conductivity drop and gap conductance; upgrade CHF/DNB correlation (e.g., Groeneveld/Chen) parameterized by pressure, mass flux, and quality, then calibrate margins.
+- [ ] Transient protection ladder: add dP/dt and dT/dt trips for SG overfill/depressurization and rod-run-in alarms; implement graded warn/arm/trip stages surfaced on the dashboard.
+- [x] Chemistry & fouling: track dissolved oxygen/boron and corrosion/fouling that degrade HX efficiency and condenser vacuum; let feedwater temperature/chemistry affect steam purity/back-pressure.
+- [x] Balance-of-plant dynamics: steam-drum level controller with shrink/swell, feedwater valve model, turbine throttle governor/overspeed trip, and improved load-follow tied to grid demand ramps.
+- [ ] Neutronics/feedback smoothing: add detector/measurement lag, fuel→clad→coolant transport delays, shared boron trim for fine regulation, and retune rod gains/rate limits to reduce power hunting while keeping safety margins intact.
+- [ ] Component wear & maintenance: make wear depend on duty cycle and off-nominal conditions (cavitation, high ΔP, high temp); add preventive maintenance scheduling and show next-due in the dashboard.
+- [ ] Scenarios & tooling: presets for cold start, load-follow, and fault injection (pump fail, relief stuck) with seedable randomness; snapshot diff tooling to compare saved states.
+- [ ] Incremental realism plan:
+  - [x] Add stored enthalpy for primary/secondary loops and a steam-drum mass/energy balance (sensible + latent) while keeping existing pump logic and tests passing. Target representative PWR conditions: primary 15–16 MPa, 290–320 °C inlet/320–330 °C outlet, secondary saturation ~6–7 MPa with boil at ~490–510 K.
+  - [x] Adjust HX/pressure handling to use stored energy (saturation clamp and pressure rise) and validate steam formation with both pumps at ~3 GW. Use realistic tube-side material assumptions (Inconel 690/SS cladding) and clamp steam quality to phase-equilibrium enthalpy.
+  - [x] Update turbine power mapping to consume steam enthalpy/quality and align protection trips with real steam presence; drive inlet steam around 6–7 MPa, quality/enthalpy-based flow to ~550–600 MW(e) per machine class if steam is available.
+  - [x] Add integration test: cold start → gens/pumps 2/2 → ramp to ~3 GW → confirm steam quality threshold at the secondary drum → enable all turbines and require electrical output. Include a step that tolerates one secondary pump off for a period to prove redundancy still yields steam.
+- [x] Dashboard follow-ups: replace turbine “Steam P” with a more useful steam availability signal (enthalpy × steam flow).
+- [x] Relief modeling: vent both loops gradually to ~1 MPa when reliefs are open, removing steam enthalpy/mass and capping pump targets to prevent instant repressurization.

pyproject.toml

@@ -10,6 +10,10 @@ dependencies = []
 
 [project.optional-dependencies]
 dev = ["pytest>=7.0"]
+dashboard = [
+    "rich>=13.7.0",
+    "textual>=0.50.0",
+]
 
 [build-system]
 requires = ["setuptools>=61"]

src/reactor_sim/commands.py

@@ -25,6 +25,8 @@ class ReactorCommand:
     secondary_pumps: dict[int, bool] | None = None
     generator_units: dict[int, bool] | None = None
     generator_auto: bool | None = None
+    primary_relief: bool | None = None
+    secondary_relief: bool | None = None
     maintenance_components: tuple[str, ...] = tuple()
 
     @classmethod

src/reactor_sim/constants.py

@@ -7,12 +7,18 @@ MEGAWATT = 1_000_000.0
 NEUTRON_LIFETIME = 0.1  # seconds, prompt neutron lifetime surrogate
 FUEL_ENERGY_DENSITY = 200.0 * MEGAWATT  # J/kg released as heat
 COOLANT_HEAT_CAPACITY = 4_200.0  # J/(kg*K) for water/steam
-COOLANT_DENSITY = 700.0  # kg/m^3 averaged between phases
+COOLANT_DENSITY = 720.0  # kg/m^3 averaged between phases
+STEAM_LATENT_HEAT = 2_200_000.0  # J/kg approximate latent heat of vaporization
 CORE_MELTDOWN_TEMPERATURE = 2_873.0  # K (approx 2600C) threshold for irreversible meltdown
 MAX_CORE_TEMPERATURE = CORE_MELTDOWN_TEMPERATURE  # Allow simulation to approach meltdown temperature
-MAX_PRESSURE = 15.0  # MPa typical PWR primary loop limit
+MAX_PRESSURE = 16.0  # MPa PWR primary loop limit
+CLAD_MAX_TEMPERATURE = 1_200.0  # K clad softening / DNB concern
+CHF_MASS_FLUX_REF = 1_500.0  # kg/m2-s reference mass flux surrogate
+CHF_PRESSURE_REF_MPA = 7.0  # MPa reference pressure for CHF surrogate
 CONTROL_ROD_SPEED = 0.03  # fraction insertion per second
 CONTROL_ROD_WORTH = 0.042  # delta rho contribution when fully withdrawn
+CONTROL_ROD_BANK_WEIGHTS = (0.4, 0.35, 0.25)
+ROD_MANUAL_STEP = 0.025
 STEAM_TURBINE_EFFICIENCY = 0.34
 GENERATOR_EFFICIENCY = 0.96
 ENVIRONMENT_TEMPERATURE = 295.0  # K
@@ -20,11 +26,69 @@ 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
 PUMP_SPOOL_TIME = 5.0  # seconds to reach commanded flow
+PRIMARY_PUMP_SHUTOFF_HEAD_MPA = 17.0  # approximate shutoff head for primary pumps
+SECONDARY_PUMP_SHUTOFF_HEAD_MPA = 8.0
 TURBINE_SPOOL_TIME = 12.0  # seconds to reach steady output
+
+# Turbine/condenser parameters
+TURBINE_THROTTLE_MIN = 0.1
+TURBINE_THROTTLE_MAX = 1.0
+TURBINE_THROTTLE_EFFICIENCY_DROP = 0.15  # efficiency loss when at minimum throttle
+CONDENSER_BASE_PRESSURE_MPA = 0.01
+CONDENSER_MAX_PRESSURE_MPA = 0.3
+CONDENSER_BACKPRESSURE_PENALTY = 0.35  # fractional power loss at max back-pressure
+CONDENSER_VACUUM_PUMP_RATE = 0.05  # MPa per second drawdown toward base when below max load
+CONDENSER_COOLING_WATER_TEMP_K = 295.0  # cooling sink temperature
+CONDENSER_FOULING_RATE = 0.00002  # incremental penalty per second of hot operation
+CONDENSER_FOULING_MAX_PENALTY = 0.2  # max additional backpressure penalty from fouling
+CONDENSER_CHEM_FOULING_RATE = 0.0005  # per-second fouling increment scaled by impurity ppm
+CONDENSER_CHEM_BACKPRESSURE_FACTOR = 0.0002  # MPa increase per ppm impurities toward condenser pressure
 GENERATOR_SPOOL_TIME = 10.0  # seconds to reach full output
 # Auxiliary power assumptions
 PUMP_POWER_MW = 12.0  # MW draw per pump unit
 BASE_AUX_LOAD_MW = 5.0  # control, instrumentation, misc.
+NORMAL_CORE_POWER_MW = 3_000.0
+TEST_MAX_POWER_MW = 4_000.0
+PRIMARY_OUTLET_TARGET_K = 580.0
+SECONDARY_OUTLET_TARGET_K = 520.0
+PRIMARY_NOMINAL_PRESSURE = 15.5  # MPa PWR primary pressure
+SECONDARY_NOMINAL_PRESSURE = 6.5  # MPa steam drum/steam line pressure surrogate
+STEAM_GENERATOR_UA_MW_PER_K = 150.0  # overall UA for steam generator (MW/K)
+# Loop volume / inventory assumptions
+PRIMARY_LOOP_VOLUME_M3 = 350.0
+SECONDARY_LOOP_VOLUME_M3 = 320.0
+PRIMARY_PRESSURIZER_SETPOINT_MPA = 15.5
+PRIMARY_PRESSURIZER_DEADBAND_MPA = 0.2
+PRIMARY_PRESSURIZER_HEAT_RATE_MPA_PER_S = 0.08
+PRIMARY_PRESSURIZER_SPRAY_RATE_MPA_PER_S = 0.12
+PRIMARY_PRESSURIZER_LEVEL_DRAW_PER_S = 0.002
+PRIMARY_PRESSURIZER_LEVEL_FILL_PER_S = 0.001
+LOOP_INVENTORY_CORRECTION_RATE = 0.08  # fraction of nominal mass restored per second toward target
+PRIMARY_INVENTORY_TARGET = 0.95
+SECONDARY_INVENTORY_TARGET = 0.9
+SECONDARY_STEAM_LOSS_FRACTION = 0.02  # fraction of steam mass that leaves the loop each second
+NPSH_REQUIRED_MPA = 0.25
+LOW_LEVEL_FLOW_FLOOR = 0.05
+# Chemistry & fouling
+CHEM_MAX_PPM = 5_000.0
+CHEM_OXYGEN_DEFAULT_PPM = 50.0  # deoxygenated feedwater target (ppb -> ppm surrogate)
+CHEM_BORON_DEFAULT_PPM = 500.0
+CHEM_SODIUM_DEFAULT_PPM = 5.0
+HX_FOULING_RATE = 1e-5  # fouling increment per second scaled by impurities/temp
+HX_FOULING_HEAL_RATE = 5e-6  # cleaning/settling when cool/low steam
+HX_FOULING_MAX_PENALTY = 0.25  # fractional UA loss cap
+BORON_WORTH_PER_PPM = 8e-6  # delta rho per ppm relative to baseline boron
+BORON_TRIM_RATE_PPM_PER_S = 0.04  # slow boron trim toward setpoint when near target
+# Mild thermal/measurement lags
+FUEL_TO_CLAD_TIME_CONSTANT = 0.3  # seconds (mild lag)
+CLAD_TO_COOLANT_TIME_CONSTANT = 0.2  # seconds (mild lag)
+POWER_MEASUREMENT_TIME_CONSTANT = 10.0  # seconds
+# Passive cooldown
+PASSIVE_COOL_RATE_PRIMARY = 0.05  # K/s toward ambient when low/no transfer
+PASSIVE_COOL_RATE_SECONDARY = 0.08  # K/s toward ambient when no steam/heat sink
+RHR_ACTIVE = True
+RHR_CUTOFF_POWER_MW = 5.0
+RHR_COOL_RATE = 0.2  # K/s forced cooldown when power is near zero
 # Threshold inventories (event counts) for flagging common poisons in diagnostics.
 KEY_POISON_THRESHOLDS = {
     "Xe": 1e20,  # xenon

src/reactor_sim/control.py

@@ -2,7 +2,7 @@
 
 from __future__ import annotations
 
-from dataclasses import dataclass
+from dataclasses import dataclass, field
 import json
 import logging
 from pathlib import Path
@@ -22,30 +22,103 @@ class ControlSystem:
     setpoint_mw: float = 3_000.0
     rod_fraction: float = 0.5
     manual_control: bool = False
+    rod_banks: list[float] = field(default_factory=lambda: [0.5, 0.5, 0.5])
+    rod_target: float = 0.5
+    _filtered_power_mw: float = 0.0
+    _integral_error: float = 0.0
+    _ramp_start_mw: float = 0.0
+    _ramp_progress_mw: float = 0.0
+    _last_manual: bool = False
 
     def update_rods(self, state: CoreState, dt: float) -> float:
+        if not self.rod_banks or len(self.rod_banks) != len(constants.CONTROL_ROD_BANK_WEIGHTS):
+            self.rod_banks = [self.rod_fraction] * len(constants.CONTROL_ROD_BANK_WEIGHTS)
+        # Keep manual tweaks in sync with the target.
+        self.rod_target = clamp(self.rod_target, 0.0, 0.95)
         if self.manual_control:
+            if abs(self.rod_fraction - self.effective_insertion()) > 1e-6:
+                self.rod_target = clamp(self.rod_fraction, 0.0, 0.95)
+            self._advance_banks(self.rod_target, dt)
+            self._last_manual = True
             return self.rod_fraction
-        error = (state.power_output_mw - self.setpoint_mw) / self.setpoint_mw
+
-        # When power is low (negative error) withdraw rods; when high, insert them.
+        # On transition from manual -> auto, start a gentle setpoint ramp from current power.
-        adjustment = error * 0.2
+        if self._last_manual:
-        adjustment = clamp(adjustment, -constants.CONTROL_ROD_SPEED * dt, constants.CONTROL_ROD_SPEED * dt)
+            self._ramp_start_mw = max(0.0, state.power_output_mw)
-        previous = self.rod_fraction
+            self._ramp_progress_mw = 0.0
-        self.rod_fraction = clamp(self.rod_fraction + adjustment, 0.0, 0.95)
+        self._last_manual = False
-        LOGGER.debug("Control rods %.3f -> %.3f (error=%.3f)", previous, self.rod_fraction, error)
+
+        # Setpoint ramp: close the gap at a limited MW/s.
+        ramp_rate = 5.0  # MW per second toward setpoint
+        effective_setpoint = self.setpoint_mw
+        if self._ramp_start_mw > 0.0 and self.setpoint_mw > self._ramp_start_mw:
+            self._ramp_progress_mw = min(
+                self.setpoint_mw - self._ramp_start_mw,
+                self._ramp_progress_mw + ramp_rate * dt,
+            )
+            effective_setpoint = self._ramp_start_mw + self._ramp_progress_mw
+        else:
+            effective_setpoint = self.setpoint_mw
+
+        raw_power = state.power_output_mw
+        # Begin filtering once we're in the vicinity of the setpoint to avoid chasing noise.
+        if self._filtered_power_mw <= 0.0:
+            self._filtered_power_mw = raw_power
+        if raw_power > 0.7 * self.setpoint_mw:
+            tau = constants.POWER_MEASUREMENT_TIME_CONSTANT
+            alpha = clamp(dt / max(1e-6, tau), 0.0, 1.0)
+            self._filtered_power_mw += alpha * (raw_power - self._filtered_power_mw)
+            measured_power = self._filtered_power_mw
+        else:
+            measured_power = raw_power
+
+        error = (measured_power - effective_setpoint) / max(1e-6, effective_setpoint)
+        near = measured_power > 0.9 * effective_setpoint
+        # Deadband near setpoint to prevent dithering; only apply when we're close to target.
+        if near and abs(error) < 0.01:
+            self._advance_banks(self.rod_target, dt)
+            return self.rod_fraction
+        # Integrate a bit to remove steady-state error.
+        self._integral_error = clamp(self._integral_error + error * dt, -0.05, 0.05)
+        p_gain = 0.35 if not near else 0.2
+        i_gain = 0.015 if not near else 0.01
+        adjustment = p_gain * error + i_gain * self._integral_error
+        speed = constants.CONTROL_ROD_SPEED * dt
+        # Hard high-band clamp: above 5% over setpoint, freeze withdrawals; above 10%, insert faster.
+        if measured_power > 1.1 * effective_setpoint:
+            speed *= 0.4
+            adjustment = -speed  # force insertion at capped rate
+        elif measured_power > 1.05 * effective_setpoint:
+            adjustment = min(adjustment, 0.0)
+            speed *= 0.6
+        # Soft clamp above setpoint: if we're >5% high, enforce insertion at capped rate.
+        if error > 0.1:
+            adjustment = min(adjustment, 0.0)
+            speed *= 0.5
+        elif error > 0.05:
+            adjustment = min(adjustment, adjustment)
+            speed *= 0.7
+        if near:
+            speed *= 0.5  # slow down when near setpoint to avoid overshoot
+        adjustment = clamp(adjustment, -speed, speed)
+        self.rod_target = clamp(self.rod_target + adjustment, 0.0, 0.95)
+        self._advance_banks(self.rod_target, dt)
+        LOGGER.debug("Control rod target=%.3f (error=%.3f)", self.rod_target, error)
         return self.rod_fraction
 
     def set_rods(self, fraction: float) -> float:
-        previous = self.rod_fraction
+        self.rod_target = self._quantize_manual(fraction)
-        self.rod_fraction = clamp(fraction, 0.0, 0.95)
+        self._advance_banks(self.rod_target, 0.0)
-        LOGGER.info("Manual rod set %.3f -> %.3f", previous, self.rod_fraction)
+        LOGGER.info("Manual rod target set to %.3f", self.rod_target)
-        return self.rod_fraction
+        return self.rod_target
 
     def increment_rods(self, delta: float) -> float:
         return self.set_rods(self.rod_fraction + delta)
 
     def scram(self) -> float:
-        self.rod_fraction = 0.95
+        self.rod_target = 0.95
+        self.rod_banks = [0.95 for _ in self.rod_banks]
+        self._sync_fraction()
         LOGGER.warning("SCRAM: rods fully inserted")
         return self.rod_fraction
 
@@ -57,15 +130,88 @@ class ControlSystem:
     def set_manual_mode(self, manual: bool) -> None:
         if self.manual_control != manual:
             self.manual_control = manual
-            LOGGER.info("Rod control %s", "manual" if manual else "automatic")
+        LOGGER.info("Rod control %s", "manual" if manual else "automatic")
 
-    def coolant_demand(self, primary: CoolantLoopState) -> float:
+    def safety_backoff(self, subcooling_margin: float | None, dnb_margin: float | None, dt: float) -> None:
-        desired_temp = 580.0
+        """Insert rods proactively when thermal margins are thin."""
-        error = (primary.temperature_out - desired_temp) / 100.0
+        if self.manual_control:
-        demand = clamp(0.8 - error, 0.0, 1.0)
+            return
-        LOGGER.debug("Coolant demand %.2f for outlet %.1fK", demand, primary.temperature_out)
+        severity = 0.0
+        if subcooling_margin is not None:
+            severity = max(severity, max(0.0, 3.0 - subcooling_margin) / 3.0)
+        if dnb_margin is not None:
+            severity = max(severity, max(0.0, 0.7 - dnb_margin) / 0.7)
+        if severity <= 0.0:
+            return
+        backoff = (0.001 + 0.01 * severity) * dt
+        self.rod_target = clamp(self.rod_target + backoff, 0.0, 0.95)
+        self._advance_banks(self.rod_target, dt)
+        LOGGER.debug("Safety backoff applied: target=%.3f severity=%.2f", self.rod_target, severity)
+
+    def coolant_demand(
+        self,
+        primary: CoolantLoopState,
+        core_power_mw: float | None = None,
+        electrical_output_mw: float | None = None,
+    ) -> float:
+        desired_temp = constants.PRIMARY_OUTLET_TARGET_K
+        # Increase demand when outlet is hotter than desired, reduce when cooler.
+        temp_error = (primary.temperature_out - desired_temp) / 100.0
+        demand = 0.8 + temp_error
+        # Keep a light power-proportional floor so both pumps stay spinning without flooding the loop.
+        power_floor = 0.0
+        if core_power_mw is not None:
+            power_fraction = clamp(core_power_mw / constants.NORMAL_CORE_POWER_MW, 0.0, 1.5)
+            power_floor = 0.2 + 0.25 * power_fraction
+        demand = max(demand, power_floor)
+        # At power, keep primary pumps near full speed to preserve pressure/subcooling.
+        if core_power_mw is not None and core_power_mw > 500.0:
+            demand = max(demand, 0.8)
+        elif core_power_mw is not None and core_power_mw > 100.0:
+            demand = max(demand, 0.6)
+        demand = clamp(demand, 0.0, 1.0)
+        LOGGER.debug(
+            "Coolant demand %.2f (temp_error=%.2f, power_floor=%.2f) for outlet %.1fK power %.1f MW elec %.1f MW",
+            demand,
+            temp_error,
+            power_floor,
+            primary.temperature_out,
+            core_power_mw or 0.0,
+            electrical_output_mw or 0.0,
+        )
         return demand
 
+    def effective_insertion(self) -> float:
+        if not self.rod_banks:
+            return self.rod_fraction
+        weights = constants.CONTROL_ROD_BANK_WEIGHTS
+        total = sum(weights)
+        effective = sum(w * b for w, b in zip(weights, self.rod_banks)) / total
+        return clamp(effective, 0.0, 0.95)
+
+    def _advance_banks(self, target: float, dt: float) -> None:
+        speed = constants.CONTROL_ROD_SPEED * dt
+        new_banks: list[float] = []
+        for idx, pos in enumerate(self.rod_banks):
+            direction = 1 if target > pos else -1
+            step = direction * speed
+            updated = clamp(pos + step, 0.0, 0.95)
+            # Avoid overshoot
+            if (direction > 0 and updated > target) or (direction < 0 and updated < target):
+                updated = target
+            new_banks.append(updated)
+        self.rod_banks = new_banks
+        self._sync_fraction()
+
+    def _sync_fraction(self) -> None:
+        self.rod_fraction = self.effective_insertion()
+
+    def _quantize_manual(self, fraction: float) -> float:
+        step = constants.ROD_MANUAL_STEP
+        quantized = round(fraction / step) * step
+        return clamp(quantized, 0.0, 0.95)
+
+
     def save_state(
         self,
         filepath: str,
@@ -78,6 +224,8 @@ class ControlSystem:
                 "setpoint_mw": self.setpoint_mw,
                 "rod_fraction": self.rod_fraction,
                 "manual_control": self.manual_control,
+                "rod_banks": self.rod_banks,
+                "rod_target": self.rod_target,
             },
             "plant": plant_state.to_dict(),
             "metadata": metadata or {},
@@ -96,6 +244,9 @@ class ControlSystem:
         self.setpoint_mw = control.get("setpoint_mw", self.setpoint_mw)
         self.rod_fraction = control.get("rod_fraction", self.rod_fraction)
         self.manual_control = control.get("manual_control", self.manual_control)
+        self.rod_banks = control.get("rod_banks", self.rod_banks) or self.rod_banks
+        self.rod_target = control.get("rod_target", self.rod_fraction)
+        self._sync_fraction()
         plant = PlantState.from_dict(data["plant"])
         LOGGER.info("Loaded plant state from %s", path)
         return plant, data.get("metadata", {}), data.get("health")

src/reactor_sim/coolant.py

@@ -15,12 +15,22 @@ LOGGER = logging.getLogger(__name__)
 class Pump:
     nominal_flow: float
     efficiency: float = 0.9
+    shutoff_head_mpa: float = 6.0
     spool_time: float = constants.PUMP_SPOOL_TIME
 
     def flow_rate(self, demand: float) -> float:
         demand = max(0.0, min(1.0, demand))
         return self.nominal_flow * (0.2 + 0.8 * demand) * self.efficiency
 
+    def performance(self, demand: float) -> tuple[float, float]:
+        """Return (flow_kg_s, head_mpa) at the given demand using a simple pump curve."""
+        demand = max(0.0, min(1.0, demand))
+        flow = self.flow_rate(demand)
+        flow_frac = flow / max(1e-3, self.nominal_flow)
+        # Keep a healthy head near nominal flow; fall off gently beyond the rated point.
+        head = self.shutoff_head_mpa * max(0.2, 1.0 - 0.4 * max(0.0, flow_frac))
+        return flow, head
+
     def step(self, loop: CoolantLoopState, demand: float) -> None:
         loop.mass_flow_rate = self.flow_rate(demand)
         loop.pressure = 12.0 * demand + 2.0

src/reactor_sim/dashboard.py

@@ -13,11 +13,45 @@ from . import constants
 from .commands import ReactorCommand
 from .reactor import Reactor
 from .simulation import ReactorSimulation
-from .state import PlantState
+from .state import PlantState, PumpState
 
 LOGGER = logging.getLogger(__name__)
 
 
+def _build_numpad_mapping() -> dict[int, float]:
+    # Use keypad matrix constants when available; skip missing ones to avoid import errors on some terminals.
+    mapping: dict[int, float] = {}
+    table = {
+        "KEY_C1": 0.1,  # numpad 1
+        "KEY_C2": 0.2,  # numpad 2
+        "KEY_C3": 0.3,  # numpad 3
+        "KEY_B1": 0.4,  # numpad 4
+        "KEY_B2": 0.5,  # numpad 5
+        "KEY_B3": 0.6,  # numpad 6
+        "KEY_A1": 0.7,  # numpad 7
+        "KEY_A2": 0.8,  # numpad 8
+        "KEY_A3": 0.9,  # numpad 9
+        # Common keypad aliases when NumLock is on
+        "KEY_END": 0.1,
+        "KEY_DOWN": 0.2,
+        "KEY_NPAGE": 0.3,
+        "KEY_LEFT": 0.4,
+        "KEY_B2": 0.5,  # center stays 0.5
+        "KEY_RIGHT": 0.6,
+        "KEY_HOME": 0.7,
+        "KEY_UP": 0.8,
+        "KEY_PPAGE": 0.9,
+    }
+    for name, value in table.items():
+        code = getattr(curses, name, None)
+        if code is not None:
+            mapping[code] = value
+    return mapping
+
+
+_NUMPAD_ROD_KEYS = _build_numpad_mapping()
+
+
 @dataclass
 class DashboardKey:
     key: str
@@ -42,8 +76,10 @@ class ReactorDashboard:
         self.sim: Optional[ReactorSimulation] = None
         self.quit_requested = False
         self.reset_requested = False
+        self.page = 1  # 1=metrics, 2=schematic (placeholder)
         self._last_state: Optional[PlantState] = None
-        self.log_buffer: deque[str] = deque(maxlen=4)
+        self._trend_history: deque[tuple[float, float, float]] = deque(maxlen=120)
+        self.log_buffer: deque[str] = deque(maxlen=8)
         self._log_handler: Optional[logging.Handler] = None
         self._previous_handlers: list[logging.Handler] = []
         self._logger = logging.getLogger("reactor_sim")
@@ -55,9 +91,12 @@ class ReactorDashboard:
                     DashboardKey("space", "SCRAM"),
                     DashboardKey("r", "Reset & clear state"),
                     DashboardKey("a", "Toggle auto rod control"),
+                    DashboardKey("F1/F2", "Metrics / schematic views"),
                     DashboardKey("+/-", "Withdraw/insert rods"),
+                    DashboardKey("1-9 / Numpad", "Set rods to 0.1 … 0.9 (manual)"),
                     DashboardKey("[/]", "Adjust consumer demand −/+50 MW"),
                     DashboardKey("s/d", "Setpoint −/+250 MW"),
+                    DashboardKey("p", "Maintain core (shutdown required)"),
                 ],
             ),
             (
@@ -74,6 +113,7 @@ class ReactorDashboard:
                 [
                     DashboardKey("b/v", "Toggle generator 1/2"),
                     DashboardKey("x", "Toggle generator auto"),
+                    DashboardKey("l/;", "Toggle relief primary/secondary"),
                     DashboardKey("B/V", "Maintain generator 1/2"),
                 ],
             ),
@@ -81,7 +121,7 @@ class ReactorDashboard:
                 "Turbines / Grid",
                 [
                     DashboardKey("t", "Toggle turbine bank"),
-                    DashboardKey("1/2/3", "Toggle turbine units 1-3"),
+                    DashboardKey("Shift+1/2/3", "Toggle turbine units 1-3"),
                     DashboardKey("y/u/i", "Maintain turbine 1/2/3"),
                     DashboardKey("c", "Toggle consumer"),
                 ],
@@ -100,6 +140,7 @@ class ReactorDashboard:
         curses.init_pair(3, curses.COLOR_GREEN, -1)
         curses.init_pair(4, curses.COLOR_RED, -1)
         stdscr.nodelay(True)
+        stdscr.keypad(True)
         self._install_log_capture()
         try:
             while True:
@@ -140,14 +181,16 @@ class ReactorDashboard:
             ch = stdscr.getch()
             if ch == -1:
                 break
+            keyname = None
+            try:
+                keyname = curses.keyname(ch)
+            except curses.error:
+                keyname = None
             if ch in (ord("q"), ord("Q")):
                 self.quit_requested = True
                 return
             if ch == ord(" "):
                 self._queue_command(ReactorCommand.scram_all())
-            elif ch in (ord("p"), ord("P")):
-                # Deprecated master toggles ignored.
-                continue
             elif ch in (ord("o"), ord("O")):
                 continue
             elif ch in (ord("g"), ord("G")):
@@ -158,23 +201,44 @@ class ReactorDashboard:
                 self._toggle_secondary_pump_unit(0)
             elif ch in (ord("k"), ord("K")):
                 self._toggle_secondary_pump_unit(1)
+            elif ch in (ord("p"), ord("P")):
+                self._queue_command(ReactorCommand.maintain("core"))
             elif ch in (ord("b"), ord("B")):
                 self._toggle_generator_unit(0)
             elif ch in (ord("v"), ord("V")):
                 self._toggle_generator_unit(1)
+            elif ch == ord("l"):
+                self._queue_command(ReactorCommand(primary_relief=not self.reactor.primary_relief_open))
+            elif ch == ord(";"):
+                self._queue_command(ReactorCommand(secondary_relief=not self.reactor.secondary_relief_open))
             elif ch in (ord("x"), ord("X")):
                 self._queue_command(ReactorCommand(generator_auto=not self.reactor.generator_auto))
             elif ch in (ord("t"), ord("T")):
                 self._queue_command(ReactorCommand(turbine_on=not self.reactor.turbine_active))
+            elif keyname and keyname.decode(errors="ignore") in ("!", "@", "#", '"'):
+                name = keyname.decode(errors="ignore")
+                turbine_hotkeys = {"!": 0, "@": 1, "#": 2, '"': 1}
+                self._toggle_turbine_unit(turbine_hotkeys[name])
+            elif ch in (ord("!"), ord("@"), ord("#"), ord('"')):
+                turbine_hotkeys = {ord("!"): 0, ord("@"): 1, ord("#"): 2, ord('"'): 1}
+                self._toggle_turbine_unit(turbine_hotkeys[ch])
+            elif keyname and keyname.startswith(b"KP_") and keyname[-1:] in b"123456789":
+                target = (keyname[-1] - ord("0")) / 10.0  # type: ignore[arg-type]
+                self._queue_command(ReactorCommand(rod_position=target, rod_manual=True))
             elif ord("1") <= ch <= ord("9"):
-                idx = ch - ord("1")
+                target = (ch - ord("0")) / 10.0
-                self._toggle_turbine_unit(idx)
+                self._queue_command(ReactorCommand(rod_position=target, rod_manual=True))
+            elif ch in _NUMPAD_ROD_KEYS:
+                self._queue_command(ReactorCommand(rod_position=_NUMPAD_ROD_KEYS[ch], rod_manual=True))
+            elif curses.KEY_F1 <= ch <= curses.KEY_F9:
+                target = (ch - curses.KEY_F1 + 1) / 10.0
+                self._queue_command(ReactorCommand(rod_position=target, rod_manual=True))
             elif ch in (ord("+"), ord("=")):
                 # Insert rods (increase fraction)
-                self._queue_command(ReactorCommand(rod_position=self._clamped_rod(0.05)))
+                self._queue_command(ReactorCommand(rod_position=self._clamped_rod(constants.ROD_MANUAL_STEP)))
             elif ch == ord("-"):
                 # Withdraw rods (decrease fraction)
-                self._queue_command(ReactorCommand(rod_position=self._clamped_rod(-0.05)))
+                self._queue_command(ReactorCommand(rod_position=self._clamped_rod(-constants.ROD_MANUAL_STEP)))
             elif ch == ord("["):
                 demand = self._current_demand() - 50.0
                 self._queue_command(ReactorCommand(consumer_demand=max(0.0, demand)))
@@ -196,8 +260,6 @@ class ReactorDashboard:
                 self._queue_command(ReactorCommand.maintain("primary_pump_1"))
             elif ch in (ord("n"), ord("N")):
                 self._queue_command(ReactorCommand.maintain("primary_pump_2"))
-            elif ch in (ord("k"), ord("K")):
-                self._queue_command(ReactorCommand.maintain("core"))
             elif ch == ord(","):
                 self._queue_command(ReactorCommand.maintain("secondary_pump_1"))
             elif ch == ord("."):
@@ -290,7 +352,7 @@ class ReactorDashboard:
     def _draw(self, stdscr: "curses._CursesWindow", state: PlantState) -> None:
         stdscr.erase()
         height, width = stdscr.getmaxyx()
-        min_status = 4
+        min_status = 6
         if height < min_status + 12 or width < 70:
             stdscr.addstr(
                 0,
@@ -301,20 +363,22 @@ class ReactorDashboard:
             stdscr.refresh()
             return
 
-        status_height = min_status
+        log_rows = max(2, min(len(self.log_buffer) + 2, 8))
-        data_height = height - status_height
+        status_height = min(height - 1, max(min_status, min_status + log_rows))
+        data_height = max(1, height - status_height)
+        gap = 2
         right_width = max(28, width // 3)
-        left_width = width - right_width
+        left_width = width - right_width - gap
         if left_width < 50:
-            left_width = min(50, width - 18)
+            left_width = min(50, width - (18 + gap))
-            right_width = width - left_width
+            right_width = width - left_width - gap
 
         data_height = max(1, data_height)
         left_width = max(1, left_width)
         right_width = max(1, right_width)
 
         data_win = stdscr.derwin(data_height, left_width, 0, 0)
-        help_win = stdscr.derwin(data_height, right_width, 0, left_width)
+        help_win = stdscr.derwin(data_height, right_width, 0, left_width + gap)
         status_win = stdscr.derwin(status_height, width, data_height, 0)
 
         self._draw_data_panel(data_win, state)
@@ -326,6 +390,7 @@ class ReactorDashboard:
         win.erase()
         win.box()
         win.addstr(0, 2, " Plant Overview ", curses.color_pair(1) | curses.A_BOLD)
+        self._update_trends(state)
         height, width = win.getmaxyx()
         inner_height = height - 2
         inner_width = width - 2
@@ -335,6 +400,8 @@ class ReactorDashboard:
         right_win = win.derwin(inner_height, right_width, 1, 1 + left_width)
         for row in range(1, height - 1):
             win.addch(row, 1 + left_width, curses.ACS_VLINE)
+            if left_width + 1 < width - 1:
+                win.addch(row, 1 + left_width + 1, curses.ACS_VLINE)
 
         left_win.erase()
         right_win.erase()
@@ -345,15 +412,29 @@ class ReactorDashboard:
             left_y,
             "Core",
             [
-                ("Fuel Temp", f"{state.core.fuel_temperature:8.1f} K"),
+                (
-                ("Core Power", f"{state.core.power_output_mw:8.1f} MW"),
+                    "Fuel Temp",
+                    f"{state.core.fuel_temperature:6.1f} K (Max {constants.CORE_MELTDOWN_TEMPERATURE:4.0f})",
+                ),
+                (
+                    "Core Power",
+                    f"{state.core.power_output_mw:6.1f} MW (Nom {constants.NORMAL_CORE_POWER_MW:4.0f}/Max {constants.TEST_MAX_POWER_MW:4.0f})",
+                ),
                 ("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}"),
+                ("Boron", f"{state.boron_ppm:7.1f} ppm"),
+                ("P(meas)", f"{self._measured_power(state):6.1f} MW"),
             ],
         )
+        left_y = self._draw_section(
+            left_win,
+            left_y,
+            "Trends",
+            self._trend_lines(state),
+        )
         left_y = self._draw_section(left_win, left_y, "Key Poisons / Emitters", self._poison_lines(state))
         left_y = self._draw_section(
             left_win,
@@ -362,10 +443,16 @@ class ReactorDashboard:
             [
                 ("Pump1", self._pump_status(state.primary_pumps, 0)),
                 ("Pump2", self._pump_status(state.primary_pumps, 1)),
-                ("Flow", f"{state.primary_loop.mass_flow_rate:7.0f} kg/s"),
+                (
+                    "Flow",
+                    f"{state.primary_loop.mass_flow_rate:7.0f}/{self.reactor.primary_pump.nominal_flow * len(self.reactor.primary_pump_units):.0f} kg/s",
+                ),
+                ("Level", f"{state.primary_loop.level*100:6.1f}%"),
                 ("Inlet Temp", f"{state.primary_loop.temperature_in:7.1f} K"),
-                ("Outlet Temp", f"{state.primary_loop.temperature_out:7.1f} K"),
+                ("Outlet Temp", f"{state.primary_loop.temperature_out:7.1f} K (Target {constants.PRIMARY_OUTLET_TARGET_K:4.0f})"),
-                ("Pressure", f"{state.primary_loop.pressure:5.2f} MPa"),
+                ("Pressure", f"{state.primary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa"),
+                ("Pressurizer", f"{self.reactor.pressurizer_level*100:6.1f}% @ {constants.PRIMARY_PRESSURIZER_SETPOINT_MPA:4.1f} MPa"),
+                ("Relief", "OPEN" if self.reactor.primary_relief_open else "CLOSED"),
             ],
         )
         self._draw_section(
@@ -375,10 +462,20 @@ class ReactorDashboard:
             [
                 ("Pump1", self._pump_status(state.secondary_pumps, 0)),
                 ("Pump2", self._pump_status(state.secondary_pumps, 1)),
-                ("Flow", f"{state.secondary_loop.mass_flow_rate:7.0f} kg/s"),
+                (
-                ("Inlet Temp", f"{state.secondary_loop.temperature_in:7.1f} K"),
+                    "Flow",
-                ("Pressure", f"{state.secondary_loop.pressure:5.2f} MPa"),
+                    f"{state.secondary_loop.mass_flow_rate:7.0f}/{self.reactor.secondary_pump.nominal_flow * len(self.reactor.secondary_pump_units):.0f} kg/s",
-                ("Steam Quality", f"{state.secondary_loop.steam_quality:5.2f}"),
+                ),
+                ("Level", f"{state.secondary_loop.level*100:6.1f}% (Target {constants.SECONDARY_INVENTORY_TARGET*100:4.0f}%)"),
+                (
+                    "Feedwater",
+                    f"valve {self.reactor.feedwater_valve*100:5.1f}% steam {state.secondary_loop.mass_flow_rate * max(0.0, state.secondary_loop.steam_quality):6.0f} kg/s",
+                ),
+                ("Inlet Temp", f"{state.secondary_loop.temperature_in:7.1f} K (Target {constants.SECONDARY_OUTLET_TARGET_K:4.0f})"),
+                ("Outlet Temp", f"{state.secondary_loop.temperature_out:7.1f} K (Target {constants.SECONDARY_OUTLET_TARGET_K:4.0f})"),
+                ("Pressure", f"{state.secondary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa"),
+                ("Steam Quality", f"{state.secondary_loop.steam_quality:5.2f}/1.00"),
+                ("Relief", "OPEN" if self.reactor.secondary_relief_open else "CLOSED"),
             ],
         )
 
@@ -394,48 +491,85 @@ class ReactorDashboard:
             "Turbine / Grid",
             [
                 ("Turbines", " ".join(self._turbine_status_lines())),
-                ("Unit1 Elec", f"{state.turbines[0].electrical_output_mw:7.1f} MW" if state.turbines else "n/a"),
+                ("Rated Elec", f"{len(self.reactor.turbines)*self.reactor.turbines[0].rated_output_mw:7.1f} MW"),
                 (
-                    "Unit2 Elec",
+                    "Steam",
-                    f"{state.turbines[1].electrical_output_mw:7.1f} MW" if len(state.turbines) > 1 else "n/a",
+                    f"h={state.turbines[0].steam_enthalpy:5.0f} kJ/kg avail {self._steam_available_power(state):6.1f} MW "
+                    f"flow {state.secondary_loop.mass_flow_rate * max(0.0, state.secondary_loop.steam_quality):6.0f} kg/s"
+                    if state.turbines
+                    else "n/a",
                 ),
                 (
-                    "Unit3 Elec",
+                    "Units Elec",
-                    f"{state.turbines[2].electrical_output_mw:7.1f} MW" if len(state.turbines) > 2 else "n/a",
+                    " ".join([f"{t.electrical_output_mw:6.1f}MW" for t in state.turbines]) if state.turbines else "n/a",
                 ),
-                ("Electrical", f"{state.total_electrical_output():7.1f} MW"),
+                (
-                ("Load", f"{self._total_load_supplied(state):7.1f}/{self._total_load_demand(state):7.1f} MW"),
+                    "Governor",
-                ("Consumer", f"{consumer_status}"),
+                    (
-                ("Demand", f"{consumer_demand:7.1f} MW"),
+                        f"thr {self.reactor.turbines[0].throttle:4.2f}→{self._desired_throttle(state.turbines[0]):4.2f} "
+                        f"ΔP {(state.turbines[0].load_demand_mw - state.turbines[0].electrical_output_mw):6.1f} MW"
+                    )
+                    if state.turbines
+                    else "n/a",
+                ),
+                (
+                    "Condenser",
+                    (
+                        f"P={state.turbines[0].condenser_pressure:4.2f}/{constants.CONDENSER_MAX_PRESSURE_MPA:4.2f} MPa "
+                        f"T={state.turbines[0].condenser_temperature:6.1f}K Foul={state.turbines[0].fouling_penalty*100:4.1f}%"
+                    )
+                    if state.turbines
+                    else "n/a",
+                ),
+                ("Electrical", f"{state.total_electrical_output():7.1f} MW | Load {self._total_load_supplied(state):6.1f}/{self._total_load_demand(state):6.1f} MW"),
+                ("Consumer", f"{consumer_status} demand {consumer_demand:6.1f} MW"),
             ],
         )
         right_y = self._draw_section(right_win, right_y, "Generators", self._generator_lines(state))
+        right_y = self._draw_section(right_win, right_y, "Power Stats", self._power_lines(state))
+        right_y = self._draw_section(right_win, right_y, "Heat Exchanger", self._heat_exchanger_lines(state))
+        right_y = self._draw_section(right_win, right_y, "Protections / Warnings", self._protection_lines(state))
         right_y = self._draw_section(right_win, right_y, "Maintenance", self._maintenance_lines())
         self._draw_health_bars(right_win, right_y)
 
     def _draw_help_panel(self, win: "curses._CursesWindow") -> None:
+        def _add_safe(row: int, col: int, text: str, attr: int = 0) -> bool:
+            max_y, max_x = win.getmaxyx()
+            if row >= max_y - 1 or col >= max_x - 1:
+                return False
+            clipped = text[: max(0, max_x - col - 1)]
+            try:
+                win.addstr(row, col, clipped, attr)
+            except curses.error:
+                return False
+            return True
+
         win.erase()
         win.box()
-        win.addstr(0, 2, " Controls ", curses.color_pair(1) | curses.A_BOLD)
+        _add_safe(0, 2, " Controls ", curses.color_pair(1) | curses.A_BOLD)
         y = 2
         for title, entries in self.help_sections:
-            win.addstr(y, 2, title, curses.color_pair(1) | curses.A_BOLD)
+            if not _add_safe(y, 2, title, curses.color_pair(1) | curses.A_BOLD):
+                return
             y += 1
             for entry in entries:
-                win.addstr(y, 4, f"{entry.key:<8} {entry.description}")
+                if not _add_safe(y, 4, f"{entry.key:<8} {entry.description}"):
+                    return
                 y += 1
             y += 1
-        win.addstr(y, 2, "Tips:", curses.color_pair(2) | curses.A_BOLD)
+        if not _add_safe(y, 2, "Tips:", curses.color_pair(2) | curses.A_BOLD):
+            return
         tips = [
             "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/,/. for pump maintenance; B/V for generators.",
             "Press 'r' to reset/clear state if you want a cold start.",
-            "Watch component health to avoid automatic trips.",
+            "Watch component health, DNB margin, and subcooling to avoid automatic trips.",
         ]
         for idx, tip in enumerate(tips, start=y + 2):
-            win.addstr(idx, 4, f"- {tip}")
+            if not _add_safe(idx, 4, f"- {tip}"):
+                break
 
     def _draw_status_panel(self, win: "curses._CursesWindow", state: PlantState) -> None:
         win.erase()
@@ -445,7 +579,7 @@ class ReactorDashboard:
             f"Time {state.time_elapsed:7.1f}s | Rods {self.reactor.control.rod_fraction:.3f} | "
             f"Primary {'ON' if self.reactor.primary_pump_active else 'OFF'} | "
             f"Secondary {'ON' if self.reactor.secondary_pump_active else 'OFF'} | "
-            f"Turbines {turbine_text}"
+            f"Turbines {turbine_text} | Page {'Metrics' if self.page == 1 else 'Schematic'}"
         )
         win.addstr(1, 1, msg, curses.color_pair(3))
         if self.reactor.health_monitor.failure_log:
@@ -455,7 +589,9 @@ class ReactorDashboard:
                 f"Failures: {', '.join(self.reactor.health_monitor.failure_log)}",
                 curses.color_pair(4) | curses.A_BOLD,
             )
-        log_y = 3
+            log_y = 4
+        else:
+            log_y = 3
         for record in list(self.log_buffer):
             if log_y >= win.getmaxyx()[0] - 1:
                 break
@@ -468,7 +604,7 @@ class ReactorDashboard:
         win: "curses._CursesWindow",
         start_y: int,
         title: str,
-        lines: list[tuple[str, str] | str],
+        lines: list[tuple[str, str] | tuple[str, str, int] | str],
     ) -> int:
         height, width = win.getmaxyx()
         inner_width = width - 4
@@ -479,15 +615,43 @@ class ReactorDashboard:
         for line in lines:
             if row >= height - 1:
                 break
+            attr = 0
             if isinstance(line, tuple):
-                label, value = line
+                if len(line) == 3:
+                    label, value, attr = line
+                else:
+                    label, value = line
                 text = f"{label:<18}: {value}"
             else:
                 text = line
-            win.addstr(row, 4, text[:inner_width])
+            win.addstr(row, 4, text[:inner_width], attr)
             row += 1
         return row + 1
 
+    def _flow_arrow(self, flow: float) -> str:
+        if flow > 15000:
+            return "====>"
+        if flow > 5000:
+            return "===>"
+        if flow > 500:
+            return "->"
+        return "--"
+
+    def _pump_glyph(self, pump_state: PumpState | None) -> str:
+        if pump_state is None:
+            return "·"
+        status = getattr(pump_state, "status", "OFF")
+        if status == "RUN":
+            return "▶"
+        if status == "CAV":
+            return "!"
+        if status == "STARTING":
+            return ">"
+        if status == "STOPPING":
+            return "-"
+        return "·"
+
+
     def _turbine_status_lines(self) -> list[str]:
         if not self.reactor.turbine_unit_active:
             return ["n/a"]
@@ -511,15 +675,19 @@ class ReactorDashboard:
         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]:
+        def fmt(symbol: str, label: str, qty: float) -> 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("Xe", "Xe (xenon)", getattr(state.core, "xenon_inventory", 0.0)))
-        lines.append(fmt("Sm", "Sm (samarium)"))
+        lines.append(fmt("Sm", "Sm (samarium)", inventory.get("Sm", 0.0)))
-        lines.append(fmt("I", "I (iodine)"))
+        lines.append(fmt("I", "I (iodine)", getattr(state.core, "iodine_inventory", 0.0)))
+        try:
+            xe_penalty = -self.reactor.neutronics.xenon_penalty(state.core)
+            lines.append(("Xe Δρ", f"{xe_penalty:+.4f}"))
+        except Exception:
+            pass
         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}"))
@@ -548,10 +716,140 @@ class ReactorDashboard:
         for idx, gen in enumerate(state.generators):
             status = "RUN" if gen.running else "START" if gen.starting else "OFF"
             spool = f" spool {gen.spool_remaining:4.1f}s" if gen.starting else ""
-            lines.append((f"Gen{idx + 1}", f"{status} {gen.power_output_mw:6.1f} MW{spool}"))
+            lines.append((f"Gen{idx + 1}", f"{status} {gen.power_output_mw:6.1f}/{self.reactor.generators[idx].rated_output_mw:4.0f} MW{spool}"))
-            lines.append((f"  Battery", f"{gen.battery_charge*100:5.1f}%"))
+            lines.append((f"  Battery", f"{gen.battery_charge*100:5.1f}% out {gen.battery_output_mw:4.1f} MW"))
         return lines
 
+    def _power_lines(self, state: PlantState) -> list[tuple[str, str]]:
+        draws = getattr(state, "aux_draws", {}) or {}
+        primary_nom = constants.PUMP_POWER_MW * len(self.reactor.primary_pump_units)
+        secondary_nom = constants.PUMP_POWER_MW * len(self.reactor.secondary_pump_units)
+        lines = [
+            ("Base Aux", f"{draws.get('base', 0.0):6.1f}/{constants.BASE_AUX_LOAD_MW:4.1f} MW"),
+            ("Primary Aux", f"{draws.get('primary_pumps', 0.0):6.1f}/{primary_nom:4.1f} MW"),
+            ("Secondary Aux", f"{draws.get('secondary_pumps', 0.0):6.1f}/{secondary_nom:4.1f} MW"),
+            ("Aux Demand", f"{draws.get('total_demand', 0.0):6.1f} MW"),
+            ("Aux Supplied", f"{draws.get('supplied', 0.0):6.1f} MW"),
+            ("Gen Output", f"{draws.get('generator_output', 0.0):6.1f} MW"),
+            ("Turbine Elec", f"{draws.get('turbine_output', 0.0):6.1f} MW"),
+        ]
+        return lines
+
+    def _heat_exchanger_lines(self, state: PlantState) -> list[tuple[str, str]]:
+        delta_t = getattr(state, "primary_to_secondary_delta_t", 0.0)
+        eff = getattr(state, "heat_exchanger_efficiency", 0.0)
+        hx_fouling = getattr(state, "hx_fouling", 0.0)
+        return [
+            ("ΔT (pri-sec)", f"{delta_t:6.1f} K"),
+            ("HX Eff", f"{eff*100:6.1f}%"),
+            ("Chem/Foul", f"O2 {getattr(state, 'dissolved_oxygen_ppm', 0.0):5.1f} ppm Na {getattr(state, 'sodium_ppm', 0.0):5.1f} ppm Foul {hx_fouling*100:5.1f}%"),
+        ]
+
+    def _protection_lines(self, state: PlantState) -> list[tuple[str, str]]:
+        lines: list[tuple[str, str] | tuple[str, str, int]] = []
+        lines.append(("SCRAM", "ACTIVE" if self.reactor.shutdown else "CLEAR", curses.color_pair(4) if self.reactor.shutdown else 0))
+        if self.reactor.meltdown:
+            lines.append(("Meltdown", "IN PROGRESS", curses.color_pair(4) | curses.A_BOLD))
+        cooldown_status = "Normal"
+        if state.core.power_output_mw <= constants.RHR_CUTOFF_POWER_MW and (self.reactor.primary_pump_active or self.reactor.secondary_pump_active):
+            cooldown_status = "RHR/Passive"
+        lines.append(("Cooldown", cooldown_status))
+        sec_flow_low = state.secondary_loop.mass_flow_rate <= 1.0 or not self.reactor.secondary_pump_active
+        heat_sink_risk = sec_flow_low and state.core.power_output_mw > 50.0
+        if heat_sink_risk:
+            heat_text = "TRIPPED low secondary flow >50 MW"
+            heat_attr = curses.color_pair(4) | curses.A_BOLD
+        elif sec_flow_low:
+            heat_text = "ARMED (secondary off/low flow)"
+            heat_attr = curses.color_pair(2) | curses.A_BOLD
+        else:
+            heat_text = "OK"
+            heat_attr = curses.color_pair(3)
+        lines.append(("Heat sink", heat_text, heat_attr))
+
+        draws = getattr(state, "aux_draws", {}) or {}
+        demand = draws.get("total_demand", 0.0)
+        supplied = draws.get("supplied", 0.0)
+        if demand > 0.1 and supplied + 1e-6 < demand:
+            aux_text = f"DEFICIT {supplied:5.1f}/{demand:5.1f} MW"
+            aux_attr = curses.color_pair(2) | curses.A_BOLD
+        elif demand > 0.1:
+            aux_text = f"OK {supplied:5.1f}/{demand:5.1f} MW"
+            aux_attr = curses.color_pair(3)
+        else:
+            aux_text = "Idle"
+            aux_attr = 0
+        lines.append(("Aux power", aux_text, aux_attr))
+
+        reliefs = []
+        if self.reactor.primary_relief_open:
+            reliefs.append("Primary")
+        if self.reactor.secondary_relief_open:
+            reliefs.append("Secondary")
+        relief_attr = curses.color_pair(2) | curses.A_BOLD if reliefs else 0
+        lines.append(("Relief valves", ", ".join(reliefs) if reliefs else "Closed", relief_attr))
+        lines.append(("DNB margin", f"{state.core.dnb_margin:5.2f}" if state.core.dnb_margin is not None else "n/a"))
+        lines.append(("Subcooling", f"{state.core.subcooling_margin:5.1f} K" if state.core.subcooling_margin is not None else "n/a"))
+        lines.append(
+            (
+                "SG level",
+                f"{state.secondary_loop.level*100:5.1f}%",
+            )
+        )
+        lines.append(("SG pressure", f"{state.secondary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa"))
+        lines.append(
+            (
+                "SCRAM trips",
+                "DNB<0.5 | Subcool<2K | SG lvl<5/>98% | SG P>15.2MPa",
+            )
+        )
+        return lines
+
+    def _steam_available_power(self, state: PlantState) -> float:
+        mass_flow = state.secondary_loop.mass_flow_rate * max(0.0, state.secondary_loop.steam_quality)
+        if mass_flow <= 1.0:
+            return 0.0
+        if state.turbines:
+            enthalpy_kjkg = max(0.0, state.turbines[0].steam_enthalpy)
+        else:
+            enthalpy_kjkg = (constants.STEAM_LATENT_HEAT / 1_000.0)
+        return (enthalpy_kjkg * mass_flow) / 1_000.0
+
+    def _measured_power(self, state: PlantState) -> float:
+        ctl = getattr(self, "reactor", None)
+        if ctl and getattr(ctl, "control", None):
+            filtered = getattr(ctl.control, "_filtered_power_mw", 0.0)
+            if filtered > 0.0:
+                return filtered
+        return state.core.power_output_mw
+
+    def _desired_throttle(self, turbine_state) -> float:
+        if not self.reactor.turbines:
+            return 0.0
+        turbine = self.reactor.turbines[0]
+        demand = turbine_state.load_demand_mw
+        return 0.4 if demand <= 0 else min(1.0, 0.4 + demand / max(1e-6, turbine.rated_output_mw))
+
+    def _update_trends(self, state: PlantState) -> None:
+        self._trend_history.append((state.time_elapsed, state.core.fuel_temperature, state.core.power_output_mw))
+
+    def _trend_lines(self, state: PlantState) -> list[tuple[str, str]]:
+        if len(self._trend_history) < 2:
+            return [("Fuel Temp Δ", "n/a"), ("Core Power Δ", "n/a"), ("P(meas)", "n/a")]
+        start_t, start_temp, start_power = self._trend_history[0]
+        end_t, end_temp, end_power = self._trend_history[-1]
+        duration = max(1.0, end_t - start_t)
+        temp_delta = end_temp - start_temp
+        power_delta = end_power - start_power
+        temp_rate = temp_delta / duration
+        power_rate = power_delta / duration
+        measured = getattr(self.reactor.control, "_filtered_power_mw", 0.0) if hasattr(self.reactor, "control") else 0.0
+        return [
+            ("Fuel Temp Δ", f"{end_temp:7.1f} K (Δ{temp_delta:+6.1f} / {duration:4.0f}s, {temp_rate:+5.2f}/s)"),
+            ("Core Power Δ", f"{end_power:7.1f} MW (Δ{power_delta:+6.1f} / {duration:4.0f}s, {power_rate:+5.2f}/s)"),
+            ("P(meas)", f"{measured:7.1f} MW" if measured > 0 else "n/a"),
+        ]
+
     def _draw_health_bars(self, win: "curses._CursesWindow", start_y: int) -> int:
         height, width = win.getmaxyx()
         inner_width = width - 4
@@ -592,7 +890,9 @@ class ReactorDashboard:
 
     def _clamped_rod(self, delta: float) -> float:
         new_fraction = self.reactor.control.rod_fraction + delta
-        return max(0.0, min(0.95, new_fraction))
+        step = constants.ROD_MANUAL_STEP
+        quantized = round(new_fraction / step) * step
+        return max(0.0, min(0.95, quantized))
 
     def _install_log_capture(self) -> None:
         if self._log_handler:
@@ -620,7 +920,22 @@ class _DashboardLogHandler(logging.Handler):
     def __init__(self, buffer: deque[str]) -> None:
         super().__init__()
         self.buffer = buffer
+        self._last_msg: str | None = None
+        self._repeat_count: int = 0
 
     def emit(self, record: logging.LogRecord) -> None:
         msg = self.format(record)
+        if msg == self._last_msg:
+            self._repeat_count += 1
+            if self.buffer and self.buffer[-1].startswith(self._last_msg):
+                try:
+                    self.buffer[-1] = f"{self._last_msg} (x{self._repeat_count + 1})"
+                except Exception:
+                    self.buffer.append(f"{msg} (x{self._repeat_count + 1})")
+            else:
+                self.buffer.append(f"{msg} (x{self._repeat_count + 1})")
+            return
+        else:
+            self._last_msg = msg
+            self._repeat_count = 0
         self.buffer.append(msg)

src/reactor_sim/failures.py

@@ -50,7 +50,8 @@ class ComponentHealth:
 class HealthMonitor:
     """Tracks component wear and signals failures."""
 
-    def __init__(self) -> None:
+    def __init__(self, disable_degradation: bool = False) -> None:
+        self.disable_degradation = disable_degradation
         self.components: Dict[str, ComponentHealth] = {
             "core": ComponentHealth("core"),
             "primary_pump_1": ComponentHealth("primary_pump_1"),
@@ -77,6 +78,8 @@ class HealthMonitor:
         generator_states: Iterable,
         dt: float,
     ) -> List[str]:
+        if self.disable_degradation:
+            return []
         events: list[str] = []
         turbine_flags = list(turbine_active)
         core = self.component("core")
@@ -89,12 +92,18 @@ class HealthMonitor:
         sec_units = list(secondary_units)
         prim_states = state.primary_pumps or []
         sec_states = state.secondary_pumps or []
+        primary_temp = getattr(state.primary_loop, "temperature_out", 295.0)
+        secondary_temp = getattr(state.secondary_loop, "temperature_out", 295.0)
+        primary_heat_factor = max(0.0, (primary_temp - 600.0) / 400.0)  # harsher wear only when very hot
+        secondary_heat_factor = max(0.0, (secondary_temp - 520.0) / 300.0)
         for idx, active in enumerate(prim_units):
             comp = self.component(f"primary_pump_{idx + 1}")
             if idx < len(prim_states) and active:
                 flow = prim_states[idx].flow_rate
                 flow_ratio = 0.0 if flow <= 0 else min(1.0, flow / 9_000.0)
-                comp.degrade((0.0002 + (1 - flow_ratio) * 0.005) * dt)
+                low_flow_penalty = (1 - flow_ratio) * 0.001
+                rate = (0.0001 + low_flow_penalty) * (1 + primary_heat_factor)
+                comp.degrade(rate * dt)
             else:
                 comp.degrade(0.0)
         for idx, active in enumerate(sec_units):
@@ -102,7 +111,9 @@ class HealthMonitor:
             if idx < len(sec_states) and active:
                 flow = sec_states[idx].flow_rate
                 flow_ratio = 0.0 if flow <= 0 else min(1.0, flow / 8_000.0)
-                comp.degrade((0.0002 + (1 - flow_ratio) * 0.004) * dt)
+                low_flow_penalty = (1 - flow_ratio) * 0.0008
+                rate = (0.0001 + low_flow_penalty) * (1 + secondary_heat_factor)
+                comp.degrade(rate * dt)
             else:
                 comp.degrade(0.0)
 

src/reactor_sim/generator.py

@@ -18,6 +18,7 @@ class GeneratorState:
     power_output_mw: float
     battery_charge: float
     status: str = "OFF"
+    battery_output_mw: float = 0.0
 
 
 @dataclass
@@ -48,9 +49,14 @@ class DieselGenerator:
 
     def step(self, state: GeneratorState, load_demand_mw: float, dt: float) -> float:
         """Advance generator dynamics and return delivered power."""
+        state.battery_output_mw = 0.0
         if state.starting:
             state.spool_remaining = max(0.0, state.spool_remaining - dt)
-            state.power_output_mw = self.rated_output_mw * (1.0 - state.spool_remaining / max(self.spool_time, 1e-6))
+            progress = 1.0 - state.spool_remaining / max(self.spool_time, 1e-6)
+            available = self.rated_output_mw * progress
+            delivered = min(available, max(0.0, load_demand_mw))
+            state.power_output_mw = delivered
+            state.battery_output_mw = min(0.5, delivered) if delivered > 0 else 0.0
             if state.spool_remaining <= 0.0:
                 state.starting = False
                 state.running = True
@@ -58,7 +64,8 @@ class DieselGenerator:
                 LOGGER.info("Generator online at %.1f MW", self.rated_output_mw)
         elif state.running:
             available = self.rated_output_mw
-            state.power_output_mw = min(available, load_demand_mw)
+            delivered = min(available, max(0.0, load_demand_mw))
+            state.power_output_mw = delivered
             state.status = "RUN" if state.power_output_mw > 0 else "IDLE"
         else:
             state.power_output_mw = 0.0
@@ -67,8 +74,8 @@ class DieselGenerator:
         if state.running:
             state.battery_charge = min(1.0, state.battery_charge + 0.02 * dt)
         elif state.starting:
-            state.battery_charge = max(0.0, state.battery_charge - 0.01 * dt)
+            state.battery_charge = max(0.0, state.battery_charge - 0.003 * dt)
         else:
-            state.battery_charge = max(0.0, state.battery_charge - 0.001 * dt)
+            state.battery_charge = max(0.0, state.battery_charge - 0.00005 * dt)
 
         return state.power_output_mw

src/reactor_sim/neutronics.py

@@ -7,7 +7,7 @@ import logging
 
 from . import constants
 from .fuel import fuel_reactivity_penalty
-from .state import CoreState
+from .state import CoreState, clamp
 
 LOGGER = logging.getLogger(__name__)
 
@@ -25,43 +25,135 @@ def xenon_poisoning(flux: float) -> float:
 
 @dataclass
 class NeutronDynamics:
-    beta_effective: float = 0.0065
+    base_shutdown_bias: float = -0.014
-    delayed_neutron_fraction: float = 0.0008
-    external_source_coupling: float = 1e-6
     shutdown_bias: float = -0.014
+    beta_effective: float = 0.0065
+    delayed_decay_const: float = 0.08  # 1/s effective precursor decay
+    external_source_coupling: float = 1e-6
+    iodine_yield: float = 1e-6  # inventory units per MW*s
+    iodine_decay_const: float = 1.0 / 66000.0  # ~18h
+    xenon_decay_const: float = 1.0 / 33000.0  # ~9h
+    xenon_burnout_coeff: float = 1e-13  # per n/cm2
+    xenon_reactivity_coeff: float = 0.05
 
-    def reactivity(self, state: CoreState, control_fraction: float) -> float:
+    def reactivity(self, state: CoreState, control_fraction: float, rod_banks: list[float] | None = None) -> float:
+        if rod_banks:
+            weights = constants.CONTROL_ROD_BANK_WEIGHTS
+            worth = 0.0
+            total = sum(weights)
+            for w, pos in zip(weights, rod_banks):
+                worth += w * (1.0 - clamp(pos, 0.0, 0.95) / 0.95)
+            rod_term = constants.CONTROL_ROD_WORTH * worth / total
+        else:
+            rod_term = constants.CONTROL_ROD_WORTH * (1.0 - control_fraction)
         rho = (
             self.shutdown_bias +
-            constants.CONTROL_ROD_WORTH * (1.0 - control_fraction)
+            rod_term
             + temperature_feedback(state.fuel_temperature)
             - fuel_reactivity_penalty(state.burnup)
-            - xenon_poisoning(state.neutron_flux)
+            - self.xenon_penalty(state)
         )
         return rho
 
     def flux_derivative(
-        self, state: CoreState, rho: float, external_source_rate: float = 0.0, baseline_source: float = 1e5
+        self,
+        state: CoreState,
+        rho: float,
+        delayed_source: float,
+        external_source_rate: float = 0.0,
+        baseline_source: float = 1e5,
     ) -> float:
         generation_time = constants.NEUTRON_LIFETIME
         beta = self.beta_effective
         source_term = self.external_source_coupling * external_source_rate
-        return ((rho - beta) / generation_time) * state.neutron_flux + baseline_source + source_term
+        prompt = ((rho - beta) / generation_time) * state.neutron_flux
+        return prompt + delayed_source + baseline_source + source_term
 
-    def step(self, state: CoreState, control_fraction: float, dt: float, external_source_rate: float = 0.0) -> None:
+    def step(
-        rho = self.reactivity(state, control_fraction)
+        self,
+        state: CoreState,
+        control_fraction: float,
+        dt: float,
+        external_source_rate: float = 0.0,
+        rod_banks: list[float] | None = None,
+    ) -> None:
+        rho = self.reactivity(state, control_fraction, rod_banks)
         rho = min(rho, 0.02)
         shutdown = control_fraction >= 0.95
         if shutdown:
             rho = min(rho, -0.04)
         baseline = 0.0 if shutdown else 1e5
         source = 0.0 if shutdown else external_source_rate
-        d_flux = self.flux_derivative(state, rho, source, baseline_source=baseline)
+        rod_positions = rod_banks if rod_banks else [control_fraction] * len(constants.CONTROL_ROD_BANK_WEIGHTS)
+        self._ensure_precursors(state, len(rod_positions))
+        bank_factors = self._bank_factors(rod_positions)
+        bank_betas = self._bank_betas(len(bank_factors))
+        delayed_source = self._delayed_source(state, bank_factors)
+
+        d_flux = self.flux_derivative(
+            state,
+            rho,
+            delayed_source,
+            external_source_rate=source,
+            baseline_source=baseline,
+        )
         state.neutron_flux = max(0.0, state.neutron_flux + d_flux * dt)
         state.reactivity_margin = rho
+        self._update_precursors(state, bank_factors, bank_betas, dt)
         LOGGER.debug(
             "Neutronics: rho=%.5f, flux=%.2e n/cm2/s, d_flux=%.2e",
             rho,
             state.neutron_flux,
             d_flux,
         )
+
+    def update_poisons(self, state: CoreState, dt: float) -> None:
+        prod_I = max(0.0, state.power_output_mw) * self.iodine_yield
+        decay_I = state.iodine_inventory * self.iodine_decay_const
+        state.iodine_inventory = max(0.0, state.iodine_inventory + (prod_I - decay_I) * dt)
+        prod_Xe = decay_I
+        burn_Xe = state.neutron_flux * self.xenon_burnout_coeff
+        decay_Xe = state.xenon_inventory * self.xenon_decay_const
+        state.xenon_inventory = max(0.0, state.xenon_inventory + (prod_Xe - decay_Xe - burn_Xe) * dt)
+
+    def xenon_penalty(self, state: CoreState) -> float:
+        """Return delta-rho penalty from xenon inventory (positive magnitude)."""
+        return self._xenon_penalty(state)
+
+    def _xenon_penalty(self, state: CoreState) -> float:
+        return min(0.05, state.xenon_inventory * self.xenon_reactivity_coeff)
+
+    def _bank_betas(self, bank_count: int) -> list[float]:
+        weights = list(constants.CONTROL_ROD_BANK_WEIGHTS)
+        if bank_count != len(weights):
+            weights = [1.0 for _ in range(bank_count)]
+        total = sum(weights) if weights else 1.0
+        return [self.beta_effective * (w / total) for w in weights]
+
+    def _bank_factors(self, positions: list[float]) -> list[float]:
+        factors: list[float] = []
+        for pos in positions:
+            insertion = clamp(pos, 0.0, 0.95)
+            factors.append(max(0.0, 1.0 - insertion / 0.95))
+        return factors
+
+    def _ensure_precursors(self, state: CoreState, bank_count: int) -> None:
+        if not state.delayed_precursors or len(state.delayed_precursors) != bank_count:
+            state.delayed_precursors = [0.0 for _ in range(bank_count)]
+
+    def _delayed_source(self, state: CoreState, bank_factors: list[float]) -> float:
+        decay = self.delayed_decay_const
+        return sum(decay * precursor * factor for precursor, factor in zip(state.delayed_precursors, bank_factors))
+
+    def _update_precursors(
+        self, state: CoreState, bank_factors: list[float], bank_betas: list[float], dt: float
+    ) -> None:
+        generation_time = constants.NEUTRON_LIFETIME
+        decay = self.delayed_decay_const
+        new_pools: list[float] = []
+        for precursor, factor, beta in zip(state.delayed_precursors, bank_factors, bank_betas):
+            production = (beta / generation_time) * state.neutron_flux * factor
+            loss = decay * precursor
+            updated = max(0.0, precursor + (production - loss) * dt)
+            new_pools.append(updated)
+        state.delayed_precursors = new_pools

src/reactor_sim/reactor.py

@@ -16,7 +16,7 @@ from .fuel import FuelAssembly, decay_heat_fraction
 from .generator import DieselGenerator, GeneratorState
 from .neutronics import NeutronDynamics
 from .state import CoolantLoopState, CoreState, PlantState, PumpState, TurbineState
-from .thermal import ThermalSolver, heat_transfer
+from .thermal import ThermalSolver, heat_transfer, saturation_pressure, saturation_temperature, temperature_rise
 from .turbine import SteamGenerator, Turbine
 
 LOGGER = logging.getLogger(__name__)
@@ -36,6 +36,9 @@ class Reactor:
     atomic_model: AtomicPhysics
     consumer: ElectricalConsumer | None = None
     health_monitor: HealthMonitor = field(default_factory=HealthMonitor)
+    pressurizer_level: float = 0.6
+    allow_external_aux: bool = False
+    relaxed_npsh: bool = False
     primary_pump_active: bool = True
     secondary_pump_active: bool = True
     primary_pump_units: list[bool] = field(default_factory=lambda: [True, True])
@@ -44,7 +47,9 @@ class Reactor:
     turbine_unit_active: list[bool] = field(default_factory=lambda: [True, True, True])
     shutdown: bool = False
     meltdown: bool = False
-    generator_auto: bool = True
+    generator_auto: bool = False
+    primary_relief_open: bool = False
+    secondary_relief_open: bool = False
     poison_alerts: set[str] = field(default_factory=set)
     maintenance_active: set[str] = field(default_factory=set)
 
@@ -56,6 +61,11 @@ class Reactor:
         self.turbine_active = any(self.turbine_unit_active)
         if not self.generators:
             self.generators = [DieselGenerator() for _ in range(2)]
+        # Balance-of-plant controls
+        self.feedwater_valve = 0.5
+        self._last_steam_out_kg_s = 0.0
+        # Slow chemistry/boron trim control
+        self._boron_trim_active = True
         if not self.primary_pump_units or len(self.primary_pump_units) != 2:
             self.primary_pump_units = [True, True]
         if not self.secondary_pump_units or len(self.secondary_pump_units) != 2:
@@ -68,8 +78,10 @@ class Reactor:
             fuel=FuelAssembly(enrichment=0.045, mass_kg=80_000.0, atomic_physics=atomic_model),
             neutronics=NeutronDynamics(),
             control=ControlSystem(),
-            primary_pump=Pump(nominal_flow=18_000.0),
+            primary_pump=Pump(nominal_flow=18_000.0, shutoff_head_mpa=constants.PRIMARY_PUMP_SHUTOFF_HEAD_MPA),
-            secondary_pump=Pump(nominal_flow=16_000.0, efficiency=0.85),
+            secondary_pump=Pump(
+                nominal_flow=16_000.0, efficiency=0.85, shutoff_head_mpa=constants.SECONDARY_PUMP_SHUTOFF_HEAD_MPA
+            ),
             thermal=ThermalSolver(),
             steam_generator=SteamGenerator(),
             turbines=[Turbine() for _ in range(3)],
@@ -81,20 +93,29 @@ class Reactor:
 
     def initial_state(self) -> PlantState:
         ambient = constants.ENVIRONMENT_TEMPERATURE
+        primary_nominal_mass = constants.PRIMARY_LOOP_VOLUME_M3 * constants.COOLANT_DENSITY
+        secondary_nominal_mass = constants.SECONDARY_LOOP_VOLUME_M3 * constants.COOLANT_DENSITY
+        self.pressurizer_level = 0.6
         core = CoreState(
             fuel_temperature=ambient,
             neutron_flux=1e5,
             reactivity_margin=-0.02,
             power_output_mw=0.1,
             burnup=0.0,
+            delayed_precursors=[0.0 for _ in constants.CONTROL_ROD_BANK_WEIGHTS],
             fission_product_inventory={},
             emitted_particles={},
         )
         # Default to a cold, safe configuration: rods fully inserted, manual control, pumps/turbines off.
         self.control.manual_control = True
         self.control.rod_fraction = 0.95
+        self.control.rod_banks = [0.95 for _ in self.control.rod_banks]
+        self.control.rod_target = 0.95
         self.shutdown = True
         self.meltdown = False
+        self.generator_auto = False
+        self.primary_relief_open = False
+        self.secondary_relief_open = False
         self.primary_pump_active = False
         self.secondary_pump_active = False
         self.primary_pump_units = [False] * len(self.primary_pump_units)
@@ -110,6 +131,9 @@ class Reactor:
             pressure=0.5,
             mass_flow_rate=0.0,
             steam_quality=0.0,
+            inventory_kg=primary_nominal_mass * constants.PRIMARY_INVENTORY_TARGET,
+            level=constants.PRIMARY_INVENTORY_TARGET,
+            energy_j=primary_nominal_mass * constants.PRIMARY_INVENTORY_TARGET * constants.COOLANT_HEAT_CAPACITY * ambient,
         )
         secondary = CoolantLoopState(
             temperature_in=ambient,
@@ -117,6 +141,9 @@ class Reactor:
             pressure=0.5,
             mass_flow_rate=0.0,
             steam_quality=0.0,
+            inventory_kg=secondary_nominal_mass * constants.SECONDARY_INVENTORY_TARGET,
+            level=constants.SECONDARY_INVENTORY_TARGET,
+            energy_j=secondary_nominal_mass * constants.SECONDARY_INVENTORY_TARGET * constants.COOLANT_HEAT_CAPACITY * ambient,
         )
         primary_pumps = [
             PumpState(active=self.primary_pump_active and self.primary_pump_units[idx], flow_rate=0.0, pressure=0.5)
@@ -155,7 +182,7 @@ class Reactor:
 
     def step(self, state: PlantState, dt: float, command: ReactorCommand | None = None) -> None:
         if self.shutdown:
-            rod_fraction = self.control.rod_fraction
+            rod_fraction = self.control.update_rods(state.core, dt)
         else:
             rod_fraction = self.control.update_rods(state.core, dt)
 
@@ -165,18 +192,24 @@ class Reactor:
         overrides = {}
         if command:
             overrides = self._apply_command(command, state)
-            rod_fraction = overrides.get("rod_fraction", rod_fraction)
+            if not self.shutdown and not self.control.manual_control:
+                rod_fraction = self.control.update_rods(state.core, 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)
+        self.neutronics.update_poisons(state.core, dt)
+        # Apply soluble boron reactivity bias (slow trim).
+        boron_delta = state.boron_ppm - constants.CHEM_BORON_DEFAULT_PPM
+        self.neutronics.shutdown_bias = self.neutronics.base_shutdown_bias - boron_delta * constants.BORON_WORTH_PER_PPM
+        self.neutronics.step(state.core, rod_fraction, dt, external_source_rate=decay_neutron_source, rod_banks=self.control.rod_banks)
 
         prompt_power, fission_rate, fission_event = self.fuel.prompt_energy_rate(
             state.core.neutron_flux, rod_fraction
         )
         decay_heat = decay_heat_fraction(state.core.burnup) * state.core.power_output_mw
         total_power = prompt_power + decay_heat + decay_power
+        total_power = min(total_power, constants.TEST_MAX_POWER_MW * 0.98, self.control.setpoint_mw * 1.1)
         state.core.power_output_mw = total_power
         state.core.update_burnup(dt)
         # Track fission products and emitted particles for diagnostics.
@@ -190,7 +223,18 @@ class Reactor:
         state.core.add_emitted_particles(particles)
         self._check_poison_alerts(state)
 
-        pump_demand = overrides.get("coolant_demand", self.control.coolant_demand(state.primary_loop))
+        self._update_loop_inventory(
+            state.primary_loop, constants.PRIMARY_LOOP_VOLUME_M3, constants.PRIMARY_INVENTORY_TARGET, dt
+        )
+
+        pump_demand = overrides.get(
+            "coolant_demand",
+            self.control.coolant_demand(
+                state.primary_loop,
+                state.core.power_output_mw,
+                state.total_electrical_output(),
+            ),
+        )
         self.primary_pump_active = self.primary_pump_active and any(self.primary_pump_units)
         self.secondary_pump_active = self.secondary_pump_active and any(self.secondary_pump_units)
         primary_units_active = [
@@ -201,42 +245,77 @@ class Reactor:
             self.secondary_pump_active and idx < len(self.secondary_pump_units) and self.secondary_pump_units[idx]
             for idx in range(2)
         ]
-        aux_demand = constants.BASE_AUX_LOAD_MW + constants.PUMP_POWER_MW * (
+        any_units = any(primary_units_active) or any(secondary_units_active)
-            sum(primary_units_active) + sum(secondary_units_active)
+        load_present = any_units or (self.consumer and self.consumer.online) or state.total_electrical_output() > 0.1
+        idle_core = state.core.power_output_mw < 1.0 and self.control.rod_fraction >= 0.9
+        baseline_off = (
+            (self.shutdown or idle_core)
+            and not load_present
+            and state.total_electrical_output() <= 0.01
+            and sum(primary_units_active) == 0
+            and sum(secondary_units_active) == 0
         )
+        aux_base = 0.0 if baseline_off else constants.BASE_AUX_LOAD_MW
+        aux_pump_primary = constants.PUMP_POWER_MW * sum(primary_units_active)
+        aux_pump_secondary = constants.PUMP_POWER_MW * sum(secondary_units_active)
+        aux_demand = aux_base + aux_pump_primary + aux_pump_secondary
         turbine_electrical = state.total_electrical_output()
         generator_power = self._step_generators(state, aux_demand, turbine_electrical, dt)
         aux_available = turbine_electrical + generator_power
-        power_ratio = 1.0 if aux_demand <= 0 else min(1.0, aux_available / aux_demand)
+        if self.allow_external_aux:
-        if aux_demand > 0 and aux_available < 0.5 * aux_demand:
+            aux_available = max(aux_available, aux_demand)
+        supplied = aux_available if aux_demand <= 0 else min(aux_available, aux_demand)
+        power_ratio = 1.0 if aux_demand <= 0 else min(1.0, supplied / max(1e-6, aux_demand))
+        if aux_demand > 0 and aux_available < 0.99 * aux_demand:
             LOGGER.warning("Aux power deficit: available %.1f/%.1f MW", aux_available, aux_demand)
+        state.aux_draws = {
+            "base": aux_base * power_ratio,
+            "primary_pumps": aux_pump_primary * power_ratio,
+            "secondary_pumps": aux_pump_secondary * power_ratio,
+            "total_demand": aux_demand,
+            "supplied": supplied,
+            "generator_output": generator_power,
+            "turbine_output": turbine_electrical,
+        }
 
         if self.primary_pump_active:
             total_flow = 0.0
-            target_pressure = (12.0 * pump_demand + 2.0) * power_ratio
+            base_flow, base_head = self.primary_pump.performance(pump_demand)
-            loop_pressure = 0.5
+            target_flow = base_flow * power_ratio
-            target_flow = self.primary_pump.flow_rate(pump_demand) * power_ratio
+            loop_pressure = max(
+                state.primary_loop.pressure, saturation_pressure(state.primary_loop.temperature_out), 0.1
+            )
+            target_pressure = max(0.5, base_head * power_ratio)
+            if self.primary_relief_open:
+                target_pressure = min(target_pressure, 1.0)
+            primary_flow_scale = min(
+                self._inventory_flow_scale(state.primary_loop), self._npsh_factor(state.primary_loop)
+            )
             for idx, pump_state in enumerate(state.primary_pumps):
                 unit_enabled = (
                     self.primary_pump_active and idx < len(self.primary_pump_units) and self.primary_pump_units[idx]
                 )
-                desired_flow = target_flow if unit_enabled else 0.0
+                powered = power_ratio > 0.1
+                desired_flow = target_flow * primary_flow_scale if unit_enabled else 0.0
                 desired_pressure = target_pressure if unit_enabled else 0.5
+                if not powered:
+                    desired_flow = 0.0
+                    desired_pressure = 0.5
                 pump_state.flow_rate = self._ramp_value(
                     pump_state.flow_rate, desired_flow, dt, self.primary_pump.spool_time
                 )
                 pump_state.pressure = self._ramp_value(
                     pump_state.pressure, desired_pressure, dt, self.primary_pump.spool_time
                 )
-                pump_state.active = (unit_enabled and power_ratio > 0.05) or pump_state.flow_rate > 1.0
+                pump_state.active = unit_enabled and powered and pump_state.flow_rate > 1.0
-                if unit_enabled and pump_state.flow_rate < max(1.0, desired_flow * 0.8):
+                if not powered or not unit_enabled:
+                    pump_state.status = "STOPPING" if pump_state.flow_rate > 1.0 else "OFF"
+                elif primary_flow_scale < 0.99:
+                    pump_state.status = "CAV"
+                elif pump_state.flow_rate < max(1.0, desired_flow * 0.8):
                     pump_state.status = "STARTING"
-                elif not unit_enabled and pump_state.flow_rate > 1.0:
-                    pump_state.status = "STOPPING"
-                elif pump_state.active:
-                    pump_state.status = "RUN"
                 else:
-                    pump_state.status = "OFF"
+                    pump_state.status = "RUN"
                 total_flow += pump_state.flow_rate
                 loop_pressure = max(loop_pressure, pump_state.pressure)
             state.primary_loop.mass_flow_rate = total_flow
@@ -247,44 +326,56 @@ class Reactor:
             state.primary_loop.mass_flow_rate = self._ramp_value(
                 state.primary_loop.mass_flow_rate, 0.0, dt, self.primary_pump.spool_time
             )
+            target_pressure = max(0.5, saturation_pressure(state.primary_loop.temperature_out))
             state.primary_loop.pressure = self._ramp_value(
-                state.primary_loop.pressure, 0.5, dt, self.primary_pump.spool_time
+                state.primary_loop.pressure, target_pressure, dt, self.primary_pump.spool_time
             )
             for pump_state in state.primary_pumps:
                 pump_state.active = False
                 pump_state.flow_rate = self._ramp_value(
                     pump_state.flow_rate, 0.0, dt, self.primary_pump.spool_time
                 )
-                pump_state.pressure = self._ramp_value(
+                pump_state.pressure = state.primary_loop.pressure
-                    pump_state.pressure, state.primary_loop.pressure, dt, self.primary_pump.spool_time
+                pump_state.status = "STOPPING" if pump_state.flow_rate > 0.1 else "OFF"
-                )
-                pump_state.status = "STOPPING" if pump_state.flow_rate > 1.0 else "OFF"
         if self.secondary_pump_active:
             total_flow = 0.0
-            target_pressure = 12.0 * 0.75 + 2.0
+            demand = 0.75
-            loop_pressure = 0.5
+            base_flow, base_head = self.secondary_pump.performance(demand)
-            target_flow = self.secondary_pump.flow_rate(0.75)
+            target_pressure = max(0.5, base_head * power_ratio)
+            if self.secondary_relief_open:
+                target_pressure = min(target_pressure, 1.0)
+            loop_pressure = max(
+                state.secondary_loop.pressure, saturation_pressure(state.secondary_loop.temperature_out), 0.1
+            )
+            target_flow = base_flow * power_ratio
+            secondary_flow_scale = min(
+                self._inventory_flow_scale(state.secondary_loop), self._npsh_factor(state.secondary_loop)
+            )
             for idx, pump_state in enumerate(state.secondary_pumps):
                 unit_enabled = (
                     self.secondary_pump_active and idx < len(self.secondary_pump_units) and self.secondary_pump_units[idx]
                 )
-                desired_flow = target_flow if unit_enabled else 0.0
+                powered = power_ratio > 0.1
+                desired_flow = target_flow * secondary_flow_scale if unit_enabled else 0.0
                 desired_pressure = target_pressure if unit_enabled else 0.5
+                if not powered:
+                    desired_flow = 0.0
+                    desired_pressure = 0.5
                 pump_state.flow_rate = self._ramp_value(
                     pump_state.flow_rate, desired_flow, dt, self.secondary_pump.spool_time
                 )
                 pump_state.pressure = self._ramp_value(
                     pump_state.pressure, desired_pressure, dt, self.secondary_pump.spool_time
                 )
-                pump_state.active = unit_enabled or pump_state.flow_rate > 1.0
+                pump_state.active = unit_enabled and powered and pump_state.flow_rate > 1.0
-                if unit_enabled and pump_state.flow_rate < max(1.0, desired_flow * 0.8):
+                if not powered or not unit_enabled:
+                    pump_state.status = "STOPPING" if pump_state.flow_rate > 1.0 else "OFF"
+                elif secondary_flow_scale < 0.99:
+                    pump_state.status = "CAV"
+                elif pump_state.flow_rate < max(1.0, desired_flow * 0.8):
                     pump_state.status = "STARTING"
-                elif not unit_enabled and pump_state.flow_rate > 1.0:
-                    pump_state.status = "STOPPING"
-                elif pump_state.active:
-                    pump_state.status = "RUN"
                 else:
-                    pump_state.status = "OFF"
+                    pump_state.status = "RUN"
                 total_flow += pump_state.flow_rate
                 loop_pressure = max(loop_pressure, pump_state.pressure)
             state.secondary_loop.mass_flow_rate = total_flow
@@ -296,30 +387,85 @@ class Reactor:
                 state.secondary_loop.mass_flow_rate, 0.0, dt, self.secondary_pump.spool_time
             )
             state.secondary_loop.pressure = self._ramp_value(
-                state.secondary_loop.pressure, 0.5, dt, self.secondary_pump.spool_time
+                state.secondary_loop.pressure, max(0.1, saturation_pressure(state.secondary_loop.temperature_out)), dt, self.secondary_pump.spool_time
             )
             for pump_state in state.secondary_pumps:
                 pump_state.active = False
                 pump_state.flow_rate = self._ramp_value(
                     pump_state.flow_rate, 0.0, dt, self.secondary_pump.spool_time
                 )
-                pump_state.pressure = self._ramp_value(
+                pump_state.pressure = state.secondary_loop.pressure
-                    pump_state.pressure, state.secondary_loop.pressure, dt, self.secondary_pump.spool_time
+                pump_state.status = "STOPPING" if pump_state.flow_rate > 0.1 else "OFF"
-                )
+
-                pump_state.status = "STOPPING" if pump_state.flow_rate > 1.0 else "OFF"
+        self._apply_pressurizer(state.primary_loop, dt)
 
-        self.thermal.step_core(state.core, state.primary_loop, total_power, dt)
         if not self.secondary_pump_active or state.secondary_loop.mass_flow_rate <= 1.0:
             transferred = 0.0
         else:
-            transferred = heat_transfer(state.primary_loop, state.secondary_loop, total_power)
+            transferred = heat_transfer(
-        self.thermal.step_secondary(state.secondary_loop, transferred)
+                state.primary_loop,
+                state.secondary_loop,
+                total_power,
+                fouling_factor=getattr(state, "hx_fouling", 0.0),
+            )
+        residual = max(0.0, total_power - transferred)
+        self.thermal.step_core(state.core, state.primary_loop, total_power, dt, residual_power_mw=residual)
+        self.thermal.step_secondary(state.secondary_loop, transferred, dt)
+        if self.primary_relief_open:
+            self._vent_relief(
+                state.primary_loop,
+                target_pressure=1.0,
+                vent_rate_max=0.02,
+                ramp_time=12.0,
+                dt=dt,
+            )
+            for pump_state in state.primary_pumps:
+                pump_state.pressure = state.primary_loop.pressure
+        if self.secondary_relief_open:
+            self._vent_relief(
+                state.secondary_loop,
+                target_pressure=1.0,
+                vent_rate_max=0.05,
+                ramp_time=10.0,
+                dt=dt,
+            )
+            for pump_state in state.secondary_pumps:
+                pump_state.pressure = state.secondary_loop.pressure
+        if not self.control.manual_control and not self.shutdown:
+            self.control.safety_backoff(state.core.subcooling_margin, state.core.dnb_margin, dt)
+        self._apply_secondary_boiloff(state, dt)
+        self._update_secondary_level(state, dt)
+        self._update_chemistry(state, dt)
+        self._apply_boron_trim(state, dt)
 
-        self._step_turbine_bank(state, transferred, dt)
+        steam_draw = self._step_turbine_bank(state, transferred, dt)
+        if steam_draw > 0.0:
+            self.thermal.remove_steam_energy(state.secondary_loop, steam_draw, dt)
         self._maintenance_tick(state, dt)
 
         if (not self.secondary_pump_active or state.secondary_loop.mass_flow_rate <= 1.0) and total_power > 50.0:
             self._handle_heat_sink_loss(state)
+        # SCRAM matrix: DNB, subcooling, steam generator level/pressure
+        if state.core.dnb_margin is not None and state.core.dnb_margin < 0.45:
+            LOGGER.critical("DNB margin low: %.2f, initiating SCRAM", state.core.dnb_margin)
+            self.shutdown = True
+            self.control.scram()
+        elif state.core.dnb_margin is not None and state.core.dnb_margin < 0.6:
+            LOGGER.warning("DNB margin low: %.2f", state.core.dnb_margin)
+        if state.core.subcooling_margin is not None and state.core.subcooling_margin < 2.0:
+            LOGGER.critical("Subcooling margin lost: %.1fK, initiating SCRAM", state.core.subcooling_margin)
+            self.shutdown = True
+            self.control.scram()
+        elif state.core.subcooling_margin is not None and state.core.subcooling_margin < 5.0:
+            LOGGER.warning("Subcooling margin low: %.1fK", state.core.subcooling_margin)
+        if state.secondary_loop.level < 0.05 or state.secondary_loop.level > 0.98:
+            LOGGER.critical("Secondary level out of bounds (%.1f%%), initiating SCRAM", state.secondary_loop.level * 100)
+            self.shutdown = True
+            self.control.scram()
+        if state.secondary_loop.pressure > 0.95 * constants.MAX_PRESSURE:
+            LOGGER.critical("Secondary pressure high (%.2f MPa), initiating SCRAM", state.secondary_loop.pressure)
+            self.shutdown = True
+            self.control.scram()
 
         failures = self.health_monitor.evaluate(
             state,
@@ -334,6 +480,52 @@ class Reactor:
 
         state.time_elapsed += dt
 
+        # Update inlet temperatures based on heat removed/added for next iteration.
+        env = constants.ENVIRONMENT_TEMPERATURE
+        primary_cooling = temperature_rise(transferred, state.primary_loop.mass_flow_rate)
+        if transferred <= 0.0 or state.secondary_loop.mass_flow_rate <= 1.0:
+            passive = constants.PASSIVE_COOL_RATE_PRIMARY * max(0.0, state.primary_loop.temperature_out - env) * dt
+            primary_cooling = max(primary_cooling, passive)
+        state.primary_loop.temperature_in = max(env, state.primary_loop.temperature_out - primary_cooling)
+        if state.core.power_output_mw <= constants.RHR_CUTOFF_POWER_MW and not self.turbine_active:
+            bleed = constants.RHR_COOL_RATE * dt
+            state.primary_loop.temperature_out = max(env, state.primary_loop.temperature_out - bleed)
+            state.primary_loop.temperature_in = max(env, state.primary_loop.temperature_out - bleed)
+
+        if state.secondary_loop.mass_flow_rate <= 1.0:
+            # Passive cooldown toward ambient when pumps off/low steam.
+            rhr = 0.0
+            if constants.RHR_ACTIVE and state.core.power_output_mw <= constants.RHR_CUTOFF_POWER_MW:
+                rhr = constants.RHR_COOL_RATE * dt
+            target_temp = max(env, state.secondary_loop.temperature_out - rhr)
+            state.secondary_loop.temperature_out = self._ramp_value(
+                state.secondary_loop.temperature_out, target_temp, dt, self.secondary_pump.spool_time
+            )
+            state.secondary_loop.temperature_in = state.secondary_loop.temperature_out
+        else:
+            # Allow the secondary to retain more heat so it can approach saturation and form steam.
+            excess = max(0.0, state.secondary_loop.temperature_out - env)
+            cooling_drop = min(40.0, max(10.0, 0.2 * excess))
+            state.secondary_loop.temperature_in = max(env, state.secondary_loop.temperature_out - cooling_drop)
+            if state.core.power_output_mw <= constants.RHR_CUTOFF_POWER_MW and not self.turbine_active:
+                bleed = constants.RHR_COOL_RATE * dt
+                state.secondary_loop.temperature_out = max(env, state.secondary_loop.temperature_out - bleed)
+                state.secondary_loop.temperature_in = max(env, state.secondary_loop.temperature_out - bleed)
+
+        # Keep stored energies consistent with updated temperatures/quality.
+        cp = constants.COOLANT_HEAT_CAPACITY
+        primary_avg = 0.5 * (state.primary_loop.temperature_in + state.primary_loop.temperature_out)
+        state.primary_loop.energy_j = max(0.0, state.primary_loop.inventory_kg * cp * primary_avg)
+        sat_temp_sec = saturation_temperature(max(0.05, state.secondary_loop.pressure))
+        sec_liquid_energy = state.secondary_loop.inventory_kg * cp * min(state.secondary_loop.temperature_out, sat_temp_sec)
+        sec_latent = state.secondary_loop.inventory_kg * state.secondary_loop.steam_quality * constants.STEAM_LATENT_HEAT
+        superheat = max(0.0, state.secondary_loop.temperature_out - sat_temp_sec)
+        sec_superheat = state.secondary_loop.inventory_kg * cp * superheat if state.secondary_loop.steam_quality >= 1.0 else 0.0
+        state.secondary_loop.energy_j = max(0.0, sec_liquid_energy + sec_latent + sec_superheat)
+
+        state.primary_to_secondary_delta_t = max(0.0, state.primary_loop.temperature_out - state.secondary_loop.temperature_in)
+        state.heat_exchanger_efficiency = 0.0 if total_power <= 0 else min(1.0, max(0.0, transferred / max(1e-6, total_power)))
+
         LOGGER.info(
             (
                 "t=%5.1fs rods=%.2f core_power=%.1fMW prompt=%.1fMW :: "
@@ -350,9 +542,10 @@ class Reactor:
             sum(t.load_demand_mw for t in state.turbines),
         )
 
-    def _step_turbine_bank(self, state: PlantState, steam_power_mw: float, dt: float) -> None:
+    def _step_turbine_bank(self, state: PlantState, steam_power_mw: float, dt: float) -> float:
         if not state.turbines:
-            return
+            return 0.0
+        steam_draw_mw = 0.0
         active_indices = [
             idx for idx, active in enumerate(self.turbine_unit_active) if active and idx < len(state.turbines)
         ]
@@ -362,17 +555,32 @@ class Reactor:
                 break
             turbine_state = state.turbines[idx]
             if idx in active_indices:
+                # Simple throttle map: reduce throttle when electrical demand is low, open as demand rises.
+                demand = turbine_state.load_demand_mw
+                desired = 0.4 if demand <= 0 else min(1.0, 0.4 + demand / max(1e-6, turbine.rated_output_mw))
+                # Governor: nudge throttle toward desired based on electrical error.
+                error = (demand - turbine_state.electrical_output_mw) / max(1.0, turbine.rated_output_mw)
+                turbine.throttle = max(0.3, min(1.0, turbine.throttle + (desired - turbine.throttle) * 0.5 + 0.2 * error * dt))
                 turbine.step(state.secondary_loop, turbine_state, steam_power_mw=power_per_unit, dt=dt)
+                if turbine_state.electrical_output_mw > 1.05 * turbine.rated_output_mw:
+                    LOGGER.critical("Turbine %d overspeed/overload, tripping unit", idx + 1)
+                    self._spin_down_turbine(turbine_state, dt, turbine.spool_time)
+                    turbine_state.status = "TRIPPED"
+                    self.turbine_unit_active[idx] = False
+                    continue
                 if power_per_unit <= 0.0 and turbine_state.electrical_output_mw < 0.1:
                     turbine_state.status = "OFF"
-                elif turbine_state.electrical_output_mw < max(0.5, power_per_unit * 0.5):
+                elif turbine_state.electrical_output_mw < max(0.1 * turbine.rated_output_mw, 1.0):
                     turbine_state.status = "STARTING"
                 else:
                     turbine_state.status = "RUN"
+                total_eff = max(1e-6, turbine.generator_efficiency * turbine.mechanical_efficiency)
+                steam_draw_mw += turbine_state.electrical_output_mw / total_eff
             else:
                 self._spin_down_turbine(turbine_state, dt, turbine.spool_time)
-                turbine_state.status = "STOPPING" if turbine_state.electrical_output_mw > 0 else "OFF"
+                turbine_state.status = "STOPPING" if turbine_state.electrical_output_mw > 0.1 else "OFF"
         self._dispatch_consumer_load(state, active_indices)
+        return steam_draw_mw
 
     def _reset_turbine_state(self, turbine_state: TurbineState) -> None:
         turbine_state.shaft_power_mw = 0.0
@@ -397,6 +605,8 @@ class Reactor:
         turbine_state.load_supplied_mw = self._ramp_value(
             turbine_state.load_supplied_mw, 0.0, dt, time_constant
         )
+        if turbine_state.electrical_output_mw < 0.1:
+            turbine_state.electrical_output_mw = 0.0
 
     def _dispatch_consumer_load(self, state: PlantState, active_indices: list[int]) -> None:
         total_electrical = sum(state.turbines[idx].electrical_output_mw for idx in active_indices)
@@ -418,6 +628,146 @@ class Reactor:
             turbine_state.load_demand_mw = demand_per_unit
             turbine_state.load_supplied_mw = share if demand_per_unit <= 0 else min(share, demand_per_unit)
 
+    def _nominal_inventory(self, volume_m3: float) -> float:
+        return volume_m3 * constants.COOLANT_DENSITY
+
+    def _update_loop_inventory(
+        self, loop: CoolantLoopState, volume_m3: float, target_level: float, dt: float
+    ) -> None:
+        nominal_mass = self._nominal_inventory(volume_m3)
+        if nominal_mass <= 0.0:
+            loop.level = 0.0
+            return
+        if loop.inventory_kg <= 0.0:
+            loop.inventory_kg = nominal_mass * target_level
+        current_level = loop.inventory_kg / nominal_mass
+        correction = (target_level - current_level) * constants.LOOP_INVENTORY_CORRECTION_RATE
+        loop.inventory_kg = max(0.0, loop.inventory_kg + correction * nominal_mass * dt)
+        loop.level = min(1.2, max(0.0, loop.inventory_kg / nominal_mass))
+
+    def _update_secondary_level(self, state: PlantState, dt: float) -> None:
+        """Steam drum level controller with shrink/swell and feedwater valve."""
+        loop = state.secondary_loop
+        nominal_mass = self._nominal_inventory(constants.SECONDARY_LOOP_VOLUME_M3)
+        if nominal_mass <= 0.0:
+            loop.level = 0.0
+            return
+        if loop.inventory_kg <= 0.0:
+            loop.inventory_kg = nominal_mass * constants.SECONDARY_INVENTORY_TARGET
+        current_level = loop.inventory_kg / nominal_mass
+        steam_out = loop.mass_flow_rate * max(0.0, loop.steam_quality)
+        # Shrink/swell: apparent level drops when steam draw surges.
+        swell = -0.02 * (steam_out - self._last_steam_out_kg_s) / max(1.0, nominal_mass)
+        sensed_level = current_level + swell
+        # PI-ish valve adjustment toward target level.
+        error = constants.SECONDARY_INVENTORY_TARGET - sensed_level
+        valve_delta = 0.3 * error * dt
+        self.feedwater_valve = max(0.0, min(1.0, self.feedwater_valve + valve_delta))
+        # Feedwater adds mass and energy at inlet temperature.
+        steam_factor = min(1.0, max(0.1, steam_out / max(1.0, nominal_mass * 0.1)))
+        feed_rate = (
+            self.feedwater_valve
+            * nominal_mass
+            * 0.002
+            * steam_factor
+        )  # up to ~0.2% of nominal mass per second, scaled by steam draw
+        added_mass = feed_rate * dt
+        loop.inventory_kg = max(0.0, loop.inventory_kg + added_mass)
+        cp = constants.COOLANT_HEAT_CAPACITY
+        loop.energy_j += added_mass * cp * loop.temperature_in
+        loop.level = min(1.2, max(0.0, loop.inventory_kg / nominal_mass))
+        self._last_steam_out_kg_s = steam_out
+
+    def _apply_boron_trim(self, state: PlantState, dt: float) -> None:
+        """Slow soluble boron trim to hold power near setpoint; acts only near target."""
+        if not self._boron_trim_active or self.control.manual_control or self.shutdown:
+            return
+        if state.time_elapsed < 300.0:
+            return
+        if self.control.setpoint_mw <= 0.0:
+            return
+        error = (state.core.power_output_mw - self.control.setpoint_mw) / self.control.setpoint_mw
+        if abs(error) < 0.005:
+            return
+        delta = constants.BORON_TRIM_RATE_PPM_PER_S * error * dt
+        state.boron_ppm = min(constants.CHEM_MAX_PPM, max(0.0, state.boron_ppm + delta))
+
+    def _update_chemistry(self, state: PlantState, dt: float) -> None:
+        """Track dissolved species and fouling impacts on HX and condenser."""
+        env = constants.ENVIRONMENT_TEMPERATURE
+        steam_out = state.secondary_loop.mass_flow_rate * max(0.0, state.secondary_loop.steam_quality)
+        temp = state.secondary_loop.temperature_out
+        temp_factor = max(0.0, (temp - env) / 300.0)
+        impurity_load = max(0.0, state.dissolved_oxygen_ppm + 0.5 * state.sodium_ppm)
+        fouling_rate = constants.HX_FOULING_RATE * temp_factor * impurity_load
+        heal = constants.HX_FOULING_HEAL_RATE * (1.0 if steam_out < 200.0 or temp_factor < 0.2 else 0.0)
+        state.hx_fouling = max(
+            0.0,
+            min(constants.HX_FOULING_MAX_PENALTY, state.hx_fouling + (fouling_rate - heal) * dt),
+        )
+        # Degas oxygen with steam production; small impurity ingress over time (worse when venting).
+        degas = 0.0005 * steam_out * dt / max(1.0, constants.SECONDARY_LOOP_VOLUME_M3)
+        state.dissolved_oxygen_ppm = max(0.0, state.dissolved_oxygen_ppm - degas)
+        ingress = (0.01 if self.secondary_relief_open else 0.002) * dt
+        state.sodium_ppm = min(constants.CHEM_MAX_PPM, state.sodium_ppm + ingress)
+        state.boron_ppm = max(0.0, state.boron_ppm - 0.001 * dt)
+        chem_penalty = constants.CONDENSER_CHEM_FOULING_RATE * impurity_load / 1_000.0
+        for turb_state in state.turbines:
+            turb_state.fouling_penalty = min(
+                constants.CONDENSER_FOULING_MAX_PENALTY,
+                max(0.0, turb_state.fouling_penalty + chem_penalty * dt),
+            )
+            backpressure = constants.CONDENSER_CHEM_BACKPRESSURE_FACTOR * impurity_load * dt
+            turb_state.condenser_pressure = min(
+                constants.CONDENSER_MAX_PRESSURE_MPA, turb_state.condenser_pressure + backpressure
+            )
+
+    def _inventory_flow_scale(self, loop: CoolantLoopState) -> float:
+        if loop.level <= constants.LOW_LEVEL_FLOW_FLOOR:
+            return 0.0
+        if loop.level <= 0.25:
+            return max(0.0, (loop.level - constants.LOW_LEVEL_FLOW_FLOOR) / (0.25 - constants.LOW_LEVEL_FLOW_FLOOR))
+        return 1.0
+
+    def _npsh_factor(self, loop: CoolantLoopState) -> float:
+        if self.relaxed_npsh:
+            return 1.0
+        vapor_pressure = saturation_pressure(loop.temperature_in)
+        available = max(0.0, loop.pressure - vapor_pressure)
+        if available <= 0.0:
+            return 0.001
+        return max(0.001, min(1.0, available / constants.NPSH_REQUIRED_MPA))
+
+    def _apply_pressurizer(self, primary: CoolantLoopState, dt: float) -> None:
+        if self.shutdown and primary.mass_flow_rate <= 100.0:
+            return
+        target = constants.PRIMARY_PRESSURIZER_SETPOINT_MPA
+        band = constants.PRIMARY_PRESSURIZER_DEADBAND_MPA
+        heat_rate = constants.PRIMARY_PRESSURIZER_HEAT_RATE_MPA_PER_S
+        spray_rate = constants.PRIMARY_PRESSURIZER_SPRAY_RATE_MPA_PER_S
+        if primary.pressure < target - band and self.pressurizer_level > 0.05:
+            primary.pressure = min(target, primary.pressure + heat_rate * dt)
+            self.pressurizer_level = max(0.0, self.pressurizer_level - constants.PRIMARY_PRESSURIZER_LEVEL_DRAW_PER_S * dt)
+        elif primary.pressure > target + band:
+            primary.pressure = max(target - band, primary.pressure - spray_rate * dt)
+            self.pressurizer_level = min(1.0, self.pressurizer_level + constants.PRIMARY_PRESSURIZER_LEVEL_FILL_PER_S * dt)
+        primary.pressure = min(constants.MAX_PRESSURE, max(saturation_pressure(primary.temperature_out), primary.pressure))
+
+    def _apply_secondary_boiloff(self, state: PlantState, dt: float) -> None:
+        loop = state.secondary_loop
+        if loop.mass_flow_rate <= 0.0 or loop.steam_quality <= 0.0:
+            return
+        steam_mass = loop.mass_flow_rate * loop.steam_quality * constants.SECONDARY_STEAM_LOSS_FRACTION * dt
+        if steam_mass <= 0.0:
+            return
+        prev_mass = max(1e-6, loop.inventory_kg)
+        loop.inventory_kg = max(0.0, loop.inventory_kg - steam_mass)
+        # Scale stored energy with the remaining mass to keep specific enthalpy consistent.
+        ratio = max(0.0, loop.inventory_kg) / prev_mass
+        loop.energy_j *= ratio
+        nominal = self._nominal_inventory(constants.SECONDARY_LOOP_VOLUME_M3)
+        loop.level = min(1.2, max(0.0, loop.inventory_kg / nominal)) if nominal > 0 else 0.0
+
     def _handle_failure(self, component: str) -> None:
         if component == "core":
             LOGGER.critical("Core failure detected. Initiating SCRAM.")
@@ -444,16 +794,22 @@ class Reactor:
             self._set_turbine_state(False)
         if command.generator_auto is not None:
             self.generator_auto = command.generator_auto
+        if command.primary_relief is not None:
+            self.primary_relief_open = command.primary_relief
+        if command.secondary_relief is not None:
+            self.secondary_relief_open = command.secondary_relief
         if command.power_setpoint is not None:
             self.control.set_power_setpoint(command.power_setpoint)
         if command.rod_manual is not None:
             self.control.set_manual_mode(command.rod_manual)
+            if command.rod_manual is False and not self.meltdown:
+                self.shutdown = False
         if command.rod_position is not None:
             self.control.set_manual_mode(True)
-            overrides["rod_fraction"] = self.control.set_rods(command.rod_position)
+            self.control.set_rods(command.rod_position)
             self.shutdown = self.shutdown or command.rod_position >= 0.95
         elif command.rod_step is not None:
-            overrides["rod_fraction"] = self.control.increment_rods(command.rod_step)
+            self.control.increment_rods(command.rod_step)
         if command.primary_pumps:
             for idx, flag in command.primary_pumps.items():
                 self._toggle_primary_pump_unit(idx - 1, flag)
@@ -477,6 +833,10 @@ class Reactor:
             overrides["coolant_demand"] = max(0.0, min(1.0, command.coolant_demand))
         for component in command.maintenance_components:
             self._toggle_maintenance(component)
+        if not self.meltdown and not command.scram:
+            rod_target = overrides.get("rod_fraction", self.control.rod_fraction)
+            if rod_target < 0.95:
+                self.shutdown = False
         return overrides
 
     def _set_primary_pump(self, active: bool) -> None:
@@ -547,6 +907,11 @@ class Reactor:
                     GeneratorState(running=False, starting=False, spool_remaining=0.0, power_output_mw=0.0, battery_charge=1.0)
                 )
         deficit = max(0.0, aux_demand - turbine_electric)
+        if self.generator_auto and aux_demand <= 0.0:
+            for idx, gen_state in enumerate(state.generators):
+                if gen_state.running or gen_state.starting:
+                    self.generators[idx].stop(gen_state)
+            return 0.0
         if self.generator_auto:
             if deficit > 0.0:
                 for idx, gen_state in enumerate(state.generators):
@@ -571,6 +936,47 @@ class Reactor:
             remaining = max(0.0, remaining - delivered)
         return total_power
 
+    def _vent_relief(
+        self,
+        loop: CoolantLoopState,
+        target_pressure: float,
+        vent_rate_max: float,
+        ramp_time: float,
+        dt: float,
+    ) -> None:
+        """Model relief valve venting: gradual depressurization with mass/enthalpy loss."""
+        if loop.inventory_kg <= 0.0:
+            loop.pressure = max(target_pressure, loop.pressure)
+            return
+        # Vent rate scales with overpressure; capped to keep a multi-second depressurization.
+        overfrac = max(0.0, (loop.pressure - target_pressure) / max(1e-6, loop.pressure))
+        vent_rate = min(vent_rate_max, 0.01 + vent_rate_max * overfrac)  # fraction of mass per second
+        vent_mass = min(loop.inventory_kg, loop.inventory_kg * vent_rate * dt)
+        if vent_mass > 0.0:
+            specific_enthalpy = (
+                loop.steam_quality * constants.STEAM_LATENT_HEAT
+                + constants.COOLANT_HEAT_CAPACITY * max(loop.temperature_out, constants.ENVIRONMENT_TEMPERATURE)
+            )
+            loop.inventory_kg = max(0.0, loop.inventory_kg - vent_mass)
+            loop.energy_j = max(0.0, loop.energy_j - vent_mass * specific_enthalpy)
+        # Pressure ramps toward target with the requested time constant.
+        ramp = min(1.0, dt / max(1e-6, ramp_time))
+        loop.pressure = max(target_pressure, loop.pressure - (loop.pressure - target_pressure) * ramp)
+        # Cool toward saturation at the new pressure to avoid re-pressurizing from superheat.
+        sat_target = saturation_temperature(target_pressure)
+        if loop.temperature_out > sat_target:
+            temp_drop = (loop.temperature_out - sat_target) * ramp
+            loop.temperature_out -= temp_drop
+            loop.temperature_in = min(loop.temperature_in, loop.temperature_out)
+            loop.steam_quality = 0.0
+            cp = constants.COOLANT_HEAT_CAPACITY
+            loop.energy_j = max(0.0, loop.inventory_kg * cp * loop.average_temperature())
+        # Re-resolve temperature/quality/pressure to reflect the vented state.
+        try:
+            self.thermal._resolve_secondary_state(loop)  # type: ignore[attr-defined]
+        except AttributeError:
+            pass
+
     def _set_turbine_state(self, active: bool, index: int | None = None) -> None:
         if index is None:
             for idx in range(len(self.turbine_unit_active)):
@@ -679,6 +1085,9 @@ class Reactor:
             "shutdown": self.shutdown,
             "meltdown": self.meltdown,
             "generator_auto": self.generator_auto,
+            "primary_relief_open": self.primary_relief_open,
+            "secondary_relief_open": self.secondary_relief_open,
+            "pressurizer_level": self.pressurizer_level,
             "maintenance_active": list(self.maintenance_active),
             "generators": [
                 {
@@ -712,6 +1121,9 @@ class Reactor:
         self.shutdown = metadata.get("shutdown", self.shutdown)
         self.meltdown = metadata.get("meltdown", self.meltdown)
         self.generator_auto = metadata.get("generator_auto", self.generator_auto)
+        self.primary_relief_open = metadata.get("primary_relief_open", self.primary_relief_open)
+        self.secondary_relief_open = metadata.get("secondary_relief_open", self.secondary_relief_open)
+        self.pressurizer_level = metadata.get("pressurizer_level", self.pressurizer_level)
         maint = metadata.get("maintenance_active")
         if maint is not None:
             self.maintenance_active = set(maint)

src/reactor_sim/rich_dashboard.py

@@ -0,0 +1,133 @@
+"""Alternate dashboard using rich Live rendering (non-interactive)."""
+
+from __future__ import annotations
+
+import logging
+from typing import Optional
+
+from rich.console import Group
+from rich.live import Live
+from rich.panel import Panel
+from rich.table import Table
+from rich.layout import Layout
+from rich.text import Text
+
+from .reactor import Reactor
+from .simulation import ReactorSimulation
+from .state import PlantState
+from . import constants
+
+LOGGER = logging.getLogger(__name__)
+
+
+def _table(title: str) -> Table:
+    tbl = Table.grid(expand=True, padding=(0, 1))
+    tbl.title = title
+    tbl.title_style = "bold cyan"
+    return tbl
+
+
+class RichDashboard:
+    """Read-only dashboard that refreshes plant metrics using rich."""
+
+    def __init__(
+        self,
+        reactor: Reactor,
+        start_state: Optional[PlantState],
+        timestep: float = 1.0,
+        save_path: Optional[str] = None,
+    ) -> None:
+        self.reactor = reactor
+        self.start_state = start_state
+        self.timestep = timestep
+        self.save_path = save_path
+        self.sim: Optional[ReactorSimulation] = None
+
+    def _render(self, state: PlantState) -> Layout:
+        layout = Layout()
+        layout.split_column(Layout(name="upper"), Layout(name="lower"))
+        layout["upper"].split_row(Layout(name="core"), Layout(name="loops"))
+        layout["lower"].split_row(Layout(name="turbine"), Layout(name="misc"))
+
+        core = _table("Core")
+        core.add_row(f"Power {state.core.power_output_mw:6.1f} MW", f"Fuel {state.core.fuel_temperature:6.1f} K")
+        core.add_row(f"Rods {self.reactor.control.rod_fraction:.3f}", f"Mode {'AUTO' if not self.reactor.control.manual_control else 'MAN'}")
+        core.add_row(f"Setpoint {self.reactor.control.setpoint_mw:5.0f} MW", f"Reactivity {state.core.reactivity_margin:+.4f}")
+        layout["upper"]["core"].update(Panel(core, padding=0))
+
+        loops = _table("Loops")
+        loops.add_row(
+            f"P pri {state.primary_loop.pressure:4.1f}/{constants.MAX_PRESSURE:4.1f} MPa",
+            f"Tin {state.primary_loop.temperature_in:5.1f}K",
+            f"Tout {state.primary_loop.temperature_out:5.1f}K",
+            f"Flow {state.primary_loop.mass_flow_rate:6.0f} kg/s",
+        )
+        loops.add_row(
+            f"P sec {state.secondary_loop.pressure:4.1f}/{constants.MAX_PRESSURE:4.1f} MPa",
+            f"Tin {state.secondary_loop.temperature_in:5.1f}K",
+            f"Tout {state.secondary_loop.temperature_out:5.1f}K",
+            f"q {state.secondary_loop.steam_quality:4.2f}",
+        )
+        loops.add_row(
+            f"HX ΔT {state.primary_to_secondary_delta_t:4.0f}K",
+            f"Eff {state.heat_exchanger_efficiency*100:5.1f}%",
+            f"Relief pri {'OPEN' if self.reactor.primary_relief_open else 'CLOSED'}",
+            f"Relief sec {'OPEN' if self.reactor.secondary_relief_open else 'CLOSED'}",
+        )
+        layout["upper"]["loops"].update(Panel(loops, padding=0))
+
+        turbine = _table("Turbines / Grid")
+        avail = getattr(self, "_steam_available_power", lambda s: 0.0)(state)  # type: ignore[arg-type]
+        steam_h = state.turbines[0].steam_enthalpy if state.turbines else 0.0
+        turbine.add_row(
+            f"Steam avail {avail:5.1f} MW",
+            f"h {steam_h:5.0f} kJ/kg",
+            f"Cond P {state.turbines[0].condenser_pressure:4.2f} MPa" if state.turbines else "Cond P n/a",
+        )
+        turbine.add_row(
+            f"Unit1 {state.turbines[0].electrical_output_mw:5.1f} MW" if state.turbines else "Unit1 n/a",
+            f"Unit2 {state.turbines[1].electrical_output_mw:5.1f} MW" if len(state.turbines) > 1 else "Unit2 n/a",
+            f"Unit3 {state.turbines[2].electrical_output_mw:5.1f} MW" if len(state.turbines) > 2 else "Unit3 n/a",
+        )
+        turbine.add_row(
+            f"Electrical {state.total_electrical_output():5.1f} MW",
+            f"Demand {self._total_load_demand(state):5.1f} MW",
+            f"Supplied {self._total_load_supplied(state):5.1f} MW",
+        )
+        layout["lower"]["turbine"].update(Panel(turbine, padding=0))
+
+        misc_group = []
+        misc_group.append(Text(f"Time {state.time_elapsed:6.1f}s", style="cyan"))
+        misc_group.append(Text(f"Primary pumps: {[p.status for p in state.primary_pumps] if state.primary_pumps else []}"))
+        misc_group.append(Text(f"Secondary pumps: {[p.status for p in state.secondary_pumps] if state.secondary_pumps else []}"))
+        misc_group.append(Text(f"Pressurizer level {self.reactor.pressurizer_level*100:5.1f}% @ {constants.PRIMARY_PRESSURIZER_SETPOINT_MPA:4.1f} MPa"))
+        misc_group.append(Text(f"Reliefs: pri={'OPEN' if self.reactor.primary_relief_open else 'CLOSED'} sec={'OPEN' if self.reactor.secondary_relief_open else 'CLOSED'}"))
+        if self.reactor.health_monitor.failure_log:
+            misc_group.append(Text(f"Failures: {', '.join(self.reactor.health_monitor.failure_log)}", style="bold red"))
+        layout["lower"]["misc"].update(Panel(Group(*misc_group), padding=0))
+        return layout
+
+    def _total_load_supplied(self, state: PlantState) -> float:
+        return sum(t.load_supplied_mw for t in state.turbines)
+
+    def _total_load_demand(self, state: PlantState) -> float:
+        return sum(t.load_demand_mw for t in state.turbines)
+
+    def run(self) -> None:
+        self.sim = ReactorSimulation(
+            self.reactor,
+            timestep=self.timestep,
+            duration=None,
+            realtime=True,
+        )
+        self.sim.start_state = self.start_state
+        try:
+            with Live(console=None, refresh_per_second=4) as live:
+                for state in self.sim.run():
+                    live.update(self._render(state))
+        except KeyboardInterrupt:
+            if self.sim:
+                self.sim.stop()
+        finally:
+            if self.save_path and self.sim and self.sim.last_state:
+                self.reactor.save_state(self.save_path, self.sim.last_state)

src/reactor_sim/simulation.py

@@ -95,6 +95,10 @@ def main() -> None:
     log_file = os.getenv("FISSION_LOG_FILE")
     configure_logging(log_level, log_file)
     realtime = os.getenv("FISSION_REALTIME", "0") == "1"
+    alternate_dashboard = os.getenv("ALTERNATE_DASHBOARD", "0") == "1"
+    snapshot_at_env = os.getenv("FISSION_SNAPSHOT_AT")
+    snapshot_at = float(snapshot_at_env) if snapshot_at_env else None
+    snapshot_path = os.getenv("FISSION_SNAPSHOT_PATH", "artifacts/snapshot.txt")
     duration_env = os.getenv("FISSION_SIM_DURATION")
     if duration_env:
         duration = None if duration_env.lower() in {"none", "infinite"} else float(duration_env)
@@ -116,19 +120,43 @@ def main() -> None:
     save_path = os.getenv("FISSION_SAVE_STATE") or str(state_path)
     if load_path:
         sim.start_state = reactor.load_state(load_path)
-    if dashboard_mode:
+    if dashboard_mode and snapshot_at is None:
-        from .dashboard import ReactorDashboard
+        if alternate_dashboard:
+            try:
+                from .textual_dashboard import run_textual_dashboard
+            except ImportError as exc:  # pragma: no cover - optional dependency
+                LOGGER.error("Textual dashboard requested but dependency missing: %s", exc)
+                return
+            run_textual_dashboard(
+                reactor,
+                start_state=sim.start_state,
+                timestep=sim.timestep,
+                save_path=save_path,
+            )
+            return
+        else:
+            from .dashboard import ReactorDashboard
 
-        dashboard = ReactorDashboard(
+            dashboard = ReactorDashboard(
-            reactor,
+                reactor,
-            start_state=sim.start_state,
+                start_state=sim.start_state,
-            timestep=sim.timestep,
+                timestep=sim.timestep,
-            save_path=save_path,
+                save_path=save_path,
-        )
+            )
-        dashboard.run()
+            dashboard.run()
-        return
+            return
     try:
-        if realtime:
+        if snapshot_at is not None:
+            sim.duration = snapshot_at
+            LOGGER.info("Running headless to t=%.1fs for snapshot...", snapshot_at)
+            for _ in sim.run():
+                pass
+            if sim.last_state:
+                from .snapshot import write_snapshot
+
+                write_snapshot(snapshot_path, reactor, sim.last_state)
+                LOGGER.info("Snapshot written to %s", snapshot_path)
+        elif realtime:
             LOGGER.info("Running in real-time mode (Ctrl+C to stop)...")
             for _ in sim.run():
                 pass

src/reactor_sim/snapshot.py

@@ -0,0 +1,53 @@
+"""Snapshot formatting helpers shared across dashboards."""
+
+from __future__ import annotations
+
+from pathlib import Path
+from typing import Iterable
+
+from . import constants
+from .reactor import Reactor
+from .state import PlantState
+
+
+def snapshot_lines(reactor: Reactor, state: PlantState) -> list[str]:
+    core = state.core
+    prim = state.primary_loop
+    sec = state.secondary_loop
+    t0 = state.turbines[0] if state.turbines else None
+    lines: list[str] = [
+        f"Time {state.time_elapsed:6.1f}s",
+        f"Core: {core.power_output_mw:6.1f}MW fuel {core.fuel_temperature:6.1f}K rods {reactor.control.rod_fraction:.3f} ({'AUTO' if not reactor.control.manual_control else 'MAN'})",
+        f"Primary: P={prim.pressure:4.1f}/{constants.MAX_PRESSURE:4.1f}MPa Tin={prim.temperature_in:6.1f}K Tout={prim.temperature_out:6.1f}K Flow={prim.mass_flow_rate:6.0f}kg/s",
+        f"Secondary: P={sec.pressure:4.1f}/{constants.MAX_PRESSURE:4.1f}MPa Tin={sec.temperature_in:6.1f}K Tout={sec.temperature_out:6.1f}K q={sec.steam_quality:4.2f} Flow={sec.mass_flow_rate:6.0f}kg/s",
+        f"HX ΔT={state.primary_to_secondary_delta_t:4.0f}K Eff={state.heat_exchanger_efficiency*100:5.1f}%",
+    ]
+    if t0:
+        lines.append(
+            f"Turbines: h={t0.steam_enthalpy:5.0f}kJ/kg avail={_steam_available_power(state):5.1f}MW "
+            f"CondP={t0.condenser_pressure:4.2f}/{constants.CONDENSER_MAX_PRESSURE_MPA:4.2f}MPa "
+            f"CondT={t0.condenser_temperature:6.1f}K"
+        )
+        lines.append("Outputs: " + " ".join([f"T{idx+1}:{t.electrical_output_mw:5.1f}MW" for idx, t in enumerate(state.turbines)]))
+    failures = ", ".join(reactor.health_monitor.failure_log) if reactor.health_monitor.failure_log else "None"
+    lines.append(
+        f"Status: pumps pri {[p.status for p in state.primary_pumps]} sec {[p.status for p in state.secondary_pumps]} "
+        f"relief pri={'OPEN' if reactor.primary_relief_open else 'CLOSED'} sec={'OPEN' if reactor.secondary_relief_open else 'CLOSED'} "
+        f"failures={failures}"
+    )
+    return lines
+
+
+def write_snapshot(path: Path | str, reactor: Reactor, state: PlantState) -> None:
+    p = Path(path)
+    p.parent.mkdir(parents=True, exist_ok=True)
+    p.write_text("\n".join(snapshot_lines(reactor, state)))
+
+
+def _steam_available_power(state: PlantState) -> float:
+    mass_flow = state.secondary_loop.mass_flow_rate * max(0.0, state.secondary_loop.steam_quality)
+    if mass_flow <= 1.0:
+        return 0.0
+    enthalpy = state.turbines[0].steam_enthalpy if state.turbines else (constants.STEAM_LATENT_HEAT / 1_000.0)
+    return (enthalpy * mass_flow) / 1_000.0
+

src/reactor_sim/state.py

@@ -4,6 +4,8 @@ from __future__ import annotations
 
 from dataclasses import dataclass, field, asdict
 
+from . import constants
+
 from .generator import GeneratorState
 
 
@@ -18,9 +20,27 @@ class CoreState:
     reactivity_margin: float  # delta rho
     power_output_mw: float  # MW thermal
     burnup: float  # fraction of fuel consumed
+    clad_temperature: float | None = None  # Kelvin
+    pellet_center_temperature: float | None = None  # Kelvin, peak centerline surrogate
+    gap_conductance: float = 1.0  # surrogate factor 0-1
+    dnb_margin: float | None = None  # ratio to critical heat flux surrogate
+    subcooling_margin: float | None = None  # K until boiling
+    xenon_inventory: float = 0.0
+    iodine_inventory: float = 0.0
+    delayed_precursors: list[float] = field(default_factory=list)
     fission_product_inventory: dict[str, float] = field(default_factory=dict)
     emitted_particles: dict[str, float] = field(default_factory=dict)
 
+    def __post_init__(self) -> None:
+        if self.clad_temperature is None:
+            self.clad_temperature = self.fuel_temperature
+        if self.pellet_center_temperature is None:
+            self.pellet_center_temperature = self.fuel_temperature
+        if self.dnb_margin is None:
+            self.dnb_margin = 1.0
+        if self.subcooling_margin is None:
+            self.subcooling_margin = 0.0
+
     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)
@@ -41,6 +61,9 @@ class CoolantLoopState:
     pressure: float  # MPa
     mass_flow_rate: float  # kg/s
     steam_quality: float  # fraction of vapor
+    inventory_kg: float = 0.0  # bulk mass of coolant
+    level: float = 1.0  # fraction full relative to nominal volume
+    energy_j: float = 0.0  # stored thermal/latent energy for the loop
 
     def average_temperature(self) -> float:
         return 0.5 * (self.temperature_in + self.temperature_out)
@@ -52,6 +75,8 @@ class TurbineState:
     shaft_power_mw: float
     electrical_output_mw: float
     condenser_temperature: float
+    condenser_pressure: float = constants.CONDENSER_BASE_PRESSURE_MPA
+    fouling_penalty: float = 0.0
     load_demand_mw: float = 0.0
     load_supplied_mw: float = 0.0
     status: str = "OFF"
@@ -74,6 +99,13 @@ class PlantState:
     primary_pumps: list[PumpState] = field(default_factory=list)
     secondary_pumps: list[PumpState] = field(default_factory=list)
     generators: list[GeneratorState] = field(default_factory=list)
+    aux_draws: dict[str, float] = field(default_factory=dict)
+    heat_exchanger_efficiency: float = 0.0
+    primary_to_secondary_delta_t: float = 0.0
+    dissolved_oxygen_ppm: float = constants.CHEM_OXYGEN_DEFAULT_PPM
+    boron_ppm: float = constants.CHEM_BORON_DEFAULT_PPM
+    sodium_ppm: float = constants.CHEM_SODIUM_DEFAULT_PPM
+    hx_fouling: float = 0.0
     time_elapsed: float = field(default=0.0)
 
     def snapshot(self) -> dict[str, float]:
@@ -103,6 +135,12 @@ class PlantState:
         core_blob = dict(data["core"])
         inventory = core_blob.pop("fission_product_inventory", {})
         particles = core_blob.pop("emitted_particles", {})
+        # Backwards/forwards compatibility for optional core fields.
+        core_blob.pop("dnb_margin", None)
+        core_blob.pop("subcooling_margin", None)
+        core_blob.setdefault("clad_temperature", core_blob.get("fuel_temperature", 295.0))
+        core_blob.setdefault("pellet_center_temperature", core_blob.get("fuel_temperature", 295.0))
+        core_blob.setdefault("gap_conductance", 1.0)
         turbines_blob = data.get("turbines")
         if turbines_blob is None:
             # Compatibility with previous single-turbine snapshots.
@@ -113,13 +151,38 @@ class PlantState:
         sec_pumps_blob = data.get("secondary_pumps", [])
         generators_blob = data.get("generators", [])
         generators = [GeneratorState(**g) for g in generators_blob]
+        aux_draws = data.get("aux_draws", {})
+        hx_eff = data.get("heat_exchanger_efficiency", 0.0)
+        delta_t = data.get("primary_to_secondary_delta_t", 0.0)
+        dissolved_oxygen = data.get("dissolved_oxygen_ppm", constants.CHEM_OXYGEN_DEFAULT_PPM)
+        boron_ppm = data.get("boron_ppm", constants.CHEM_BORON_DEFAULT_PPM)
+        sodium_ppm = data.get("sodium_ppm", constants.CHEM_SODIUM_DEFAULT_PPM)
+        hx_fouling = data.get("hx_fouling", 0.0)
         return cls(
             core=CoreState(**core_blob, fission_product_inventory=inventory, emitted_particles=particles),
-            primary_loop=CoolantLoopState(**data["primary_loop"]),
+            primary_loop=CoolantLoopState(**_with_energy(data["primary_loop"])),
-            secondary_loop=CoolantLoopState(**data["secondary_loop"]),
+            secondary_loop=CoolantLoopState(**_with_energy(data["secondary_loop"])),
             turbines=turbines,
             primary_pumps=[PumpState(**p) for p in prim_pumps_blob],
             secondary_pumps=[PumpState(**p) for p in sec_pumps_blob],
             generators=generators,
+            aux_draws=aux_draws,
+            heat_exchanger_efficiency=hx_eff,
+            primary_to_secondary_delta_t=delta_t,
+            dissolved_oxygen_ppm=dissolved_oxygen,
+            boron_ppm=boron_ppm,
+            sodium_ppm=sodium_ppm,
+            hx_fouling=hx_fouling,
             time_elapsed=data.get("time_elapsed", 0.0),
         )
+
+
+def _with_energy(loop_blob: dict) -> dict:
+    """Backwards compatibility: derive energy if missing."""
+    if "energy_j" in loop_blob:
+        return loop_blob
+    energy = 0.5 * (loop_blob.get("temperature_in", 295.0) + loop_blob.get("temperature_out", 295.0))
+    energy *= loop_blob.get("inventory_kg", 0.0) * constants.COOLANT_HEAT_CAPACITY
+    out = dict(loop_blob)
+    out["energy_j"] = energy
+    return out

src/reactor_sim/textual_dashboard.py

@@ -0,0 +1,603 @@
+"""Interactive Textual-based dashboard mirroring the curses layout."""
+
+from __future__ import annotations
+
+import logging
+import os
+from collections import deque
+from pathlib import Path
+from typing import Optional
+
+from textual.app import App, ComposeResult
+from textual.containers import Grid, Horizontal, HorizontalScroll, Vertical, VerticalScroll
+from textual.widgets import Button, Footer, Header, Static
+from textual.timer import Timer
+
+from . import constants
+from .reactor import Reactor
+from .simulation import ReactorSimulation
+from .state import PlantState
+from .commands import ReactorCommand
+from .snapshot import snapshot_lines
+
+LOGGER = logging.getLogger(__name__)
+
+
+def _bar(label: str, value: float, width: int = 24) -> str:
+    filled = int(max(0.0, min(1.0, value)) * width)
+    return f"{label:<14} [{'#'*filled}{'-'*(width-filled)}] {value*100:5.1f}%"
+
+
+class TextualDashboard(App):
+    """Textual dashboard with controls and sections similar to the curses view."""
+
+    CSS = """
+    Screen { layout: vertical; }
+    #controls { padding: 0 0; height: auto; }
+    #controls .row { height: auto; }
+    #controls Button { min-width: 8; padding: 0 1; }
+    .columns { height: 1fr; }
+    .col { width: 1fr; height: 1fr; }
+    .panel { padding: 0 1; border: round $primary; }
+    """
+
+    BINDINGS = [
+        ("q", "quit", "Quit"),
+        ("space", "scram", "SCRAM"),
+        ("g", "toggle_primary_p1", "P1"),
+        ("h", "toggle_primary_p2", "P2"),
+        ("j", "toggle_secondary_p1", "S1"),
+        ("k", "toggle_secondary_p2", "S2"),
+        ("b", "toggle_generator1", "Gen1"),
+        ("v", "toggle_generator2", "Gen2"),
+        ("t", "toggle_turbine_bank", "Turbines"),
+        ("1", "toggle_turbine_unit('1')", "T1"),
+        ("2", "toggle_turbine_unit('2')", "T2"),
+        ("3", "toggle_turbine_unit('3')", "T3"),
+        ("+", "rods_insert", "+ rods"),
+        ("-", "rods_withdraw", "- rods"),
+        ("[", "demand_down", "Demand -"),
+        ("]", "demand_up", "Demand +"),
+        ("s", "setpoint_down", "SP -"),
+        ("d", "setpoint_up", "SP +"),
+        ("l", "toggle_primary_relief", "Relief pri"),
+        (";", "toggle_secondary_relief", "Relief sec"),
+        ("c", "toggle_consumer", "Consumer"),
+        ("a", "toggle_auto_rods", "Auto rods"),
+        ("f12", "snapshot", "Snapshot"),
+        # Maintenance is mouse-only here; keyboard remains in curses.
+    ]
+
+    timestep: float = 1.0
+
+    def __init__(
+        self,
+        reactor: Reactor,
+        start_state: Optional[PlantState],
+        timestep: float = 1.0,
+        save_path: Optional[str] = None,
+    ) -> None:
+        super().__init__()
+        self.reactor = reactor
+        self.state = start_state or self.reactor.initial_state()
+        self.timestep = timestep
+        self.save_path = save_path
+        self._pending: deque[ReactorCommand] = deque()
+        self._timer: Optional[Timer] = None
+        self._trend_history: deque[tuple[float, float, float]] = deque(maxlen=120)
+        self.log_buffer: deque[str] = deque(maxlen=8)
+        self._log_handler: Optional[logging.Handler] = None
+        self._previous_handlers: list[logging.Handler] = []
+        snap_at_env = os.getenv("FISSION_SNAPSHOT_AT")
+        self.snapshot_at = float(snap_at_env) if snap_at_env else None
+        self.snapshot_done = False
+        self.snapshot_path = Path(os.getenv("FISSION_SNAPSHOT_PATH", "artifacts/textual_snapshot.txt"))
+
+        # Panels
+        self.core_panel = Static(classes="panel")
+        self.trend_panel = Static(classes="panel")
+        self.poison_panel = Static(classes="panel")
+        self.primary_panel = Static(classes="panel")
+        self.secondary_panel = Static(classes="panel")
+        self.turbine_panel = Static(classes="panel")
+        self.generator_panel = Static(classes="panel")
+        self.power_panel = Static(classes="panel")
+        self.hx_panel = Static(classes="panel")
+        self.protection_panel = Static(classes="panel")
+        self.maintenance_panel = Static(classes="panel")
+        self.health_panel = Static(classes="panel")
+        self.help_panel = Static(classes="panel")
+        self.log_panel = Static(classes="panel")
+        self.status_panel = Static(classes="panel")
+        self.controls_panel: Grid | None = None
+
+    def compose(self) -> ComposeResult:
+        yield Header()
+        controls = self._build_controls()
+        yield controls
+        with Horizontal(classes="columns"):
+            with VerticalScroll(classes="col"):
+                yield self.core_panel
+                yield self.trend_panel
+                yield self.poison_panel
+                yield self.primary_panel
+                yield self.secondary_panel
+                yield self.power_panel
+                yield self.generator_panel
+                yield self.status_panel
+            with VerticalScroll(classes="col"):
+                yield self.turbine_panel
+                yield self.hx_panel
+                yield self.protection_panel
+                yield self.maintenance_panel
+                yield self.health_panel
+                yield self.log_panel
+                yield self.help_panel
+        yield Footer()
+
+    def on_mount(self) -> None:
+        self._install_log_capture()
+        self._timer = self.set_interval(self.timestep, self._tick, pause=False)
+        self._render_panels()
+
+    def action_quit(self) -> None:
+        if self._timer:
+            self._timer.pause()
+        if self.save_path and self.state:
+            self.reactor.save_state(self.save_path, self.state)
+        self._restore_logging()
+        self.exit()
+
+    # Command helpers
+    def _enqueue(self, cmd: ReactorCommand) -> None:
+        self._pending.append(cmd)
+
+    def action_scram(self) -> None:
+        self._enqueue(ReactorCommand.scram_all())
+
+    def action_toggle_primary_p1(self) -> None:
+        self._enqueue(ReactorCommand(primary_pumps={1: not self.reactor.primary_pump_units[0]}))
+
+    def action_toggle_primary_p2(self) -> None:
+        self._enqueue(ReactorCommand(primary_pumps={2: not self.reactor.primary_pump_units[1]}))
+
+    def action_toggle_secondary_p1(self) -> None:
+        self._enqueue(ReactorCommand(secondary_pumps={1: not self.reactor.secondary_pump_units[0]}))
+
+    def action_toggle_secondary_p2(self) -> None:
+        self._enqueue(ReactorCommand(secondary_pumps={2: not self.reactor.secondary_pump_units[1]}))
+
+    def action_toggle_generator1(self) -> None:
+        self._enqueue(ReactorCommand(generator_units={1: True}))
+
+    def action_toggle_generator2(self) -> None:
+        self._enqueue(ReactorCommand(generator_units={2: True}))
+
+    def action_toggle_turbine_bank(self) -> None:
+        self._enqueue(ReactorCommand(turbine_on=not self.reactor.turbine_active))
+
+    def action_toggle_turbine_unit(self, unit: str) -> None:
+        idx = int(unit)
+        current = self.reactor.turbine_unit_active[idx - 1] if idx - 1 < len(self.reactor.turbine_unit_active) else False
+        self._enqueue(ReactorCommand(turbine_units={idx: not current}))
+
+    def action_rods_insert(self) -> None:
+        self._enqueue(ReactorCommand(rod_position=min(0.95, self.reactor.control.rod_fraction + constants.ROD_MANUAL_STEP)))
+
+    def action_rods_withdraw(self) -> None:
+        self._enqueue(ReactorCommand(rod_position=max(0.0, self.reactor.control.rod_fraction - constants.ROD_MANUAL_STEP)))
+
+    def action_demand_down(self) -> None:
+        if self.reactor.consumer:
+            self._enqueue(ReactorCommand(consumer_demand=max(0.0, self.reactor.consumer.demand_mw - 50.0)))
+
+    def action_demand_up(self) -> None:
+        if self.reactor.consumer:
+            self._enqueue(ReactorCommand(consumer_demand=self.reactor.consumer.demand_mw + 50.0))
+
+    def action_setpoint_down(self) -> None:
+        self._enqueue(ReactorCommand(power_setpoint=self.reactor.control.setpoint_mw - 250.0))
+
+    def action_setpoint_up(self) -> None:
+        self._enqueue(ReactorCommand(power_setpoint=self.reactor.control.setpoint_mw + 250.0))
+
+    def action_toggle_primary_relief(self) -> None:
+        self._enqueue(ReactorCommand(primary_relief=not self.reactor.primary_relief_open))
+
+    def action_toggle_secondary_relief(self) -> None:
+        self._enqueue(ReactorCommand(secondary_relief=not self.reactor.secondary_relief_open))
+
+    def action_toggle_consumer(self) -> None:
+        if self.reactor.consumer:
+            self._enqueue(ReactorCommand(consumer_online=not self.reactor.consumer.online))
+
+    def action_toggle_auto_rods(self) -> None:
+        self._enqueue(ReactorCommand(rod_manual=not self.reactor.control.manual_control))
+
+    # Maintenance (mouse-driven)
+    def action_maintain_core(self) -> None:
+        self._enqueue(ReactorCommand.maintain("core"))
+
+    def action_maintain_primary_p1(self) -> None:
+        self._enqueue(ReactorCommand.maintain("primary_pump_1"))
+
+    def action_maintain_primary_p2(self) -> None:
+        self._enqueue(ReactorCommand.maintain("primary_pump_2"))
+
+    def action_maintain_secondary_p1(self) -> None:
+        self._enqueue(ReactorCommand.maintain("secondary_pump_1"))
+
+    def action_maintain_secondary_p2(self) -> None:
+        self._enqueue(ReactorCommand.maintain("secondary_pump_2"))
+
+    def action_maintain_turbine_1(self) -> None:
+        self._enqueue(ReactorCommand.maintain("turbine_1"))
+
+    def action_maintain_turbine_2(self) -> None:
+        self._enqueue(ReactorCommand.maintain("turbine_2"))
+
+    def action_maintain_turbine_3(self) -> None:
+        self._enqueue(ReactorCommand.maintain("turbine_3"))
+
+    def action_maintain_generator_1(self) -> None:
+        self._enqueue(ReactorCommand.maintain("generator_1"))
+
+    def action_maintain_generator_2(self) -> None:
+        self._enqueue(ReactorCommand.maintain("generator_2"))
+
+    def action_snapshot(self) -> None:
+        self._save_snapshot(auto=False)
+
+    def _merge_commands(self) -> Optional[ReactorCommand]:
+        if not self._pending:
+            return None
+        cmd = ReactorCommand()
+        while self._pending:
+            nxt = self._pending.popleft()
+            for field in nxt.__dataclass_fields__:  # type: ignore[attr-defined]
+                val = getattr(nxt, field)
+                if val is None or val is False:
+                    continue
+                setattr(cmd, field, val)
+        return cmd
+
+    def _tick(self) -> None:
+        cmd = self._merge_commands()
+        self.reactor.step(self.state, self.timestep, cmd)
+        self._render_panels()
+        if self.snapshot_at is not None and not self.snapshot_done and self.state.time_elapsed >= self.snapshot_at:
+            self._save_snapshot(auto=True)
+
+    def _build_controls(self) -> Vertical:
+        row1 = Horizontal(
+            Button("SCRAM", id="scram"),
+            Button("P1", id="p1"),
+            Button("P2", id="p2"),
+            Button("S1", id="s1"),
+            Button("S2", id="s2"),
+            Button("RelP", id="relief_pri"),
+            Button("RelS", id="relief_sec"),
+            Button("Turb", id="turbines"),
+            Button("T1", id="t1"),
+            Button("T2", id="t2"),
+            Button("T3", id="t3"),
+            Button("G1", id="g1"),
+            Button("G2", id="g2"),
+            Button("Grid", id="consumer"),
+            Button("AutoR", id="auto_rods"),
+            Button("+Rod", id="rods_plus"),
+            Button("-Rod", id="rods_minus"),
+            classes="row",
+        )
+        row2 = Horizontal(
+            Button("D+50", id="demand_up"),
+            Button("D-50", id="demand_down"),
+            Button("SP+250", id="sp_up"),
+            Button("SP-250", id="sp_down"),
+            Button("Snap", id="snapshot"),
+            classes="row",
+        )
+        row3 = Horizontal(
+            Button("MCore", id="m_core"),
+            Button("MP1", id="m_p1"),
+            Button("MP2", id="m_p2"),
+            Button("MS1", id="m_s1"),
+            Button("MS2", id="m_s2"),
+            Button("MT1", id="m_t1"),
+            Button("MT2", id="m_t2"),
+            Button("MT3", id="m_t3"),
+            Button("MG1", id="m_g1"),
+            Button("MG2", id="m_g2"),
+            classes="row",
+        )
+        return Vertical(
+            HorizontalScroll(row1, classes="row"),
+            HorizontalScroll(row2, classes="row"),
+            HorizontalScroll(row3, classes="row"),
+            id="controls",
+        )
+
+    def _render_panels(self) -> None:
+        self._trend_history.append((self.state.time_elapsed, self.state.core.fuel_temperature, self.state.core.power_output_mw))
+
+        self.core_panel.update(
+            "[bold cyan]Core[/bold cyan]\n"
+            f"Power {self.state.core.power_output_mw:6.1f} MW (Nom {constants.NORMAL_CORE_POWER_MW:4.0f}/Max {constants.TEST_MAX_POWER_MW:4.0f})\n"
+            f"Fuel {self.state.core.fuel_temperature:6.1f} K  Rods {self.reactor.control.rod_fraction:.3f} ({'AUTO' if not self.reactor.control.manual_control else 'MAN'})\n"
+            f"Setpoint {self.reactor.control.setpoint_mw:5.0f} MW  Reactivity {self.state.core.reactivity_margin:+.4f}\n"
+            f"DNB {self.state.core.dnb_margin:4.2f} Subcool {self.state.core.subcooling_margin:4.1f}K"
+        )
+        self.trend_panel.update(self._trend_text())
+        self.poison_panel.update(self._poison_text())
+        self.primary_panel.update(
+            "[bold cyan]Primary Loop[/bold cyan]\n"
+            f"Flow {self.state.primary_loop.mass_flow_rate:7.0f}/{self.reactor.primary_pump.nominal_flow * len(self.reactor.primary_pump_units):.0f} kg/s\n"
+            f"Level {self.state.primary_loop.level*100:6.1f}%\n"
+            f"Tin {self.state.primary_loop.temperature_in:7.1f} K  Tout {self.state.primary_loop.temperature_out:7.1f} K (Target {constants.PRIMARY_OUTLET_TARGET_K:4.0f})\n"
+            f"P {self.state.primary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa  Pressurizer {self.reactor.pressurizer_level*100:6.1f}% @ {constants.PRIMARY_PRESSURIZER_SETPOINT_MPA:4.1f} MPa\n"
+            f"Relief {'OPEN' if self.reactor.primary_relief_open else 'CLOSED'} Pumps {[p.status for p in self.state.primary_pumps]}"
+        )
+        self.secondary_panel.update(
+            "[bold cyan]Secondary Loop[/bold cyan]\n"
+            f"Flow {self.state.secondary_loop.mass_flow_rate:7.0f}/{self.reactor.secondary_pump.nominal_flow * len(self.reactor.secondary_pump_units):.0f} kg/s\n"
+            f"Level {self.state.secondary_loop.level*100:6.1f}%\n"
+            f"Tin {self.state.secondary_loop.temperature_in:7.1f} K  Tout {self.state.secondary_loop.temperature_out:7.1f} K (Target {constants.SECONDARY_OUTLET_TARGET_K:4.0f})\n"
+            f"P {self.state.secondary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa  q {self.state.secondary_loop.steam_quality:5.2f}/1.00\n"
+            f"Relief {'OPEN' if self.reactor.secondary_relief_open else 'CLOSED'} Pumps {[p.status for p in self.state.secondary_pumps]}"
+        )
+        self.turbine_panel.update(self._turbine_text())
+        self.generator_panel.update(self._generator_text())
+        self.power_panel.update(self._power_text())
+        self.hx_panel.update(
+            "[bold cyan]Heat Exchanger[/bold cyan]\n"
+            f"ΔT (pri-sec) {self.state.primary_to_secondary_delta_t:4.0f} K\n"
+            f"Efficiency {self.state.heat_exchanger_efficiency*100:5.1f}%"
+        )
+        self.protection_panel.update(self._protection_text())
+        self.maintenance_panel.update(self._maintenance_text())
+        self.health_panel.update(self._health_text())
+        self.log_panel.update(self._log_text())
+        self.help_panel.update(self._help_text())
+        failures = ", ".join(self.reactor.health_monitor.failure_log) if self.reactor.health_monitor.failure_log else "None"
+        self.status_panel.update(
+            "[bold cyan]Status[/bold cyan]\n"
+            f"Time {self.state.time_elapsed:6.1f}s\n"
+            f"Consumer {'ON' if (self.reactor.consumer and self.reactor.consumer.online) else 'OFF'} Demand {self.reactor.consumer.demand_mw if self.reactor.consumer else 0.0:5.1f} MW\n"
+            f"Failures: {failures}"
+        )
+
+    def _trend_text(self) -> str:
+        if len(self._trend_history) < 2:
+            return "[bold cyan]Trends[/bold cyan]\nFuel Temp Δ n/a\nCore Power Δ n/a"
+        start_t, start_temp, start_power = self._trend_history[0]
+        end_t, end_temp, end_power = self._trend_history[-1]
+        duration = max(1.0, end_t - start_t)
+        temp_delta = end_temp - start_temp
+        power_delta = end_power - start_power
+        temp_rate = temp_delta / duration
+        power_rate = power_delta / duration
+        return (
+            "[bold cyan]Trends[/bold cyan]\n"
+            f"Fuel Temp Δ {end_temp:7.1f} K (Δ{temp_delta:+6.1f} / {duration:4.0f}s, {temp_rate:+5.2f}/s)\n"
+            f"Core Power Δ {end_power:7.1f} MW (Δ{power_delta:+6.1f} / {duration:4.0f}s, {power_rate:+5.2f}/s)"
+        )
+
+    def _poison_text(self) -> str:
+        inventory = self.state.core.fission_product_inventory or {}
+        particles = self.state.core.emitted_particles or {}
+        xe = getattr(self.state.core, "xenon_inventory", 0.0)
+        sm = inventory.get("Sm", 0.0)
+        iodine = inventory.get("I", 0.0)
+        xe_drho = getattr(self.state.core, "reactivity_margin", 0.0)
+        return (
+            "[bold cyan]Key Poisons / Emitters[/bold cyan]\n"
+            f"Xe (xenon): {xe:9.2e}\n"
+            f"Sm (samarium): {sm:9.2e}\n"
+            f"I (iodine): {iodine:9.2e}\n"
+            f"Xe Δρ: {xe_drho:+.4f}\n"
+            f"Neutrons (src): {particles.get('neutrons', 0.0):.2e}"
+        )
+
+    def _turbine_text(self) -> str:
+        steam_avail = self._steam_available_power(self.state)
+        enthalpy = self.state.turbines[0].steam_enthalpy if self.state.turbines else 0.0
+        cond = ""
+        if self.state.turbines:
+            cond = (
+                f"Cond P {self.state.turbines[0].condenser_pressure:4.2f}/{constants.CONDENSER_MAX_PRESSURE_MPA:4.2f} MPa "
+                f"T {self.state.turbines[0].condenser_temperature:6.1f}K"
+            )
+        lines = [
+            "[bold cyan]Turbine / Grid[/bold cyan]",
+            f"Turbines {' '.join(self._turbine_status_lines()) if self._turbine_status_lines() else 'n/a'}",
+            f"Rated Elec {len(self.reactor.turbines)*self.reactor.turbines[0].rated_output_mw:7.1f} MW" if self.reactor.turbines else "Rated Elec n/a",
+            f"Steam Avail {steam_avail:5.1f} MW h={enthalpy:5.0f} kJ/kg {cond}",
+        ]
+        if self.state.turbines:
+            lines.append("Unit Elec " + " ".join([f"{t.electrical_output_mw:6.1f}MW" for t in self.state.turbines]))
+        lines.append(f"Throttle {self.reactor.turbines[0].throttle:5.2f}" if self.reactor.turbines else "Throttle n/a")
+        lines.append(f"Electrical {self.state.total_electrical_output():7.1f} MW Load {self._total_load_supplied(self.state):7.1f}/{self._total_load_demand(self.state):7.1f} MW")
+        if self.reactor.consumer:
+            lines.append(f"Consumer {'ONLINE' if self.reactor.consumer.online else 'OFF'} Demand {self.reactor.consumer.demand_mw:7.1f} MW")
+        return "\n".join(lines)
+
+    def _generator_text(self) -> str:
+        lines = ["[bold cyan]Generators[/bold cyan]"]
+        for idx, g in enumerate(self.state.generators):
+            status = "RUN" if g.running else "OFF"
+            if g.starting:
+                status = "START"
+            lines.append(f"Gen{idx+1}: {status:5} {g.power_output_mw:5.1f} MW batt {g.battery_charge*100:5.1f}%")
+        return "\n".join(lines)
+
+    def _power_text(self) -> str:
+        draws = getattr(self.state, "aux_draws", {}) or {}
+        base = draws.get("base", 0.0)
+        prim = draws.get("primary_pumps", 0.0)
+        sec = draws.get("secondary_pumps", 0.0)
+        demand = draws.get("total_demand", 0.0)
+        supplied = draws.get("supplied", 0.0)
+        gen_out = draws.get("generator_output", 0.0)
+        turb_out = draws.get("turbine_output", 0.0)
+        return (
+            "[bold cyan]Power Stats[/bold cyan]\n"
+            f"Base Aux {base:5.1f} MW  Prim Aux {prim:5.1f} MW  Sec Aux {sec:5.1f} MW\n"
+            f"Aux demand {demand:5.1f} MW supplied {supplied:5.1f} MW\n"
+            f"Gen out {gen_out:5.1f} MW Turbine out {turb_out:5.1f} MW"
+        )
+
+    def _protection_text(self) -> str:
+        lines = ["[bold cyan]Protections / Warnings[/bold cyan]"]
+        lines.append(f"SCRAM {'ACTIVE' if self.reactor.shutdown else 'CLEAR'}")
+        if self.reactor.meltdown:
+            lines.append("Meltdown IN PROGRESS")
+        sec_flow_low = self.state.secondary_loop.mass_flow_rate <= 1.0 or not self.reactor.secondary_pump_active
+        heat_sink_risk = sec_flow_low and self.state.core.power_output_mw > 50.0
+        if heat_sink_risk:
+            heat_text = "TRIPPED low secondary flow >50 MW"
+        elif sec_flow_low:
+            heat_text = "ARMED (secondary off/low flow)"
+        else:
+            heat_text = "OK"
+        lines.append(f"Heat sink {heat_text}")
+        lines.append(f"DNB margin {self.state.core.dnb_margin:4.2f}")
+        lines.append(f"Subcooling {self.state.core.subcooling_margin:5.1f} K")
+        lines.append(f"Reliefs pri={'OPEN' if self.reactor.primary_relief_open else 'CLOSED'} sec={'OPEN' if self.reactor.secondary_relief_open else 'CLOSED'}")
+        return "\n".join(lines)
+
+    def _maintenance_text(self) -> str:
+        active = list(self.reactor.maintenance_active)
+        return "[bold cyan]Maintenance[/bold cyan]\nActive: " + (", ".join(active) if active else "None")
+
+    def _health_text(self) -> str:
+        lines = ["[bold cyan]Component Health[/bold cyan]"]
+        for name, comp in self.reactor.health_monitor.components.items():
+            lines.append(_bar(name, comp.integrity))
+        return "\n".join(lines)
+
+    def _help_text(self) -> str:
+        tips = [
+            "Start pumps before withdrawing rods.",
+            "Bring turbine/consumer online after thermal stabilization.",
+            "Toggle turbine units 1/2/3 individually.",
+            "Use m/n/,/. in curses; mapped to j/k etc here.",
+            "F12 saves a snapshot; set FISSION_SNAPSHOT_AT for auto.",
+        ]
+        return "[bold cyan]Controls & Tips[/bold cyan]\n" + "\n".join(f"- {t}" for t in tips)
+
+    def _log_text(self) -> str:
+        lines = ["[bold cyan]Logs[/bold cyan]"]
+        if not self.log_buffer:
+            lines.append("No recent logs.")
+        else:
+            lines.extend(list(self.log_buffer))
+        return "\n".join(lines)
+
+    def _turbine_status_lines(self) -> list[str]:
+        if not self.reactor.turbine_unit_active:
+            return []
+        lines: list[str] = []
+        for idx, active in enumerate(self.reactor.turbine_unit_active):
+            label = f"{idx + 1}:"
+            status = "ON" if active else "OFF"
+            if idx < len(getattr(self.state, "turbines", [])):
+                t_state = self.state.turbines[idx]
+                status = getattr(t_state, "status", status)
+            lines.append(f"{label}{status}")
+        return lines
+
+    def _total_load_supplied(self, state: PlantState) -> float:
+        return sum(t.load_supplied_mw for t in state.turbines)
+
+    def _total_load_demand(self, state: PlantState) -> float:
+        return sum(t.load_demand_mw for t in state.turbines)
+
+    def _steam_available_power(self, state: PlantState) -> float:
+        mass_flow = state.secondary_loop.mass_flow_rate * max(0.0, state.secondary_loop.steam_quality)
+        if mass_flow <= 1.0:
+            return 0.0
+        enthalpy = state.turbines[0].steam_enthalpy if state.turbines else (constants.STEAM_LATENT_HEAT / 1_000.0)
+        return (enthalpy * mass_flow) / 1_000.0
+
+    def _snapshot_lines(self) -> list[str]:
+        return snapshot_lines(self.reactor, self.state)
+
+    def _save_snapshot(self, auto: bool = False) -> None:
+        try:
+            self.snapshot_path.parent.mkdir(parents=True, exist_ok=True)
+            self.snapshot_path.write_text("\n".join(self._snapshot_lines()))
+            self.snapshot_done = True
+            LOGGER.info("Saved dashboard snapshot to %s%s", self.snapshot_path, " (auto)" if auto else "")
+        except Exception as exc:  # pragma: no cover
+            LOGGER.error("Failed to save snapshot: %s", exc)
+
+    def _install_log_capture(self) -> None:
+        # Silence existing root handlers to avoid spewing logs over the UI.
+        root = logging.getLogger()
+        root.handlers = []
+        # Capture reactor_sim logs into the on-screen log buffer.
+        logger = logging.getLogger("reactor_sim")
+        logger.propagate = False
+        self._previous_handlers = list(logger.handlers)
+        handler = logging.StreamHandler()
+        handler.setLevel(logging.INFO)
+
+        def emit(record: logging.LogRecord) -> None:
+            msg = handler.format(record)
+            self.log_buffer.append(msg)
+
+        handler.emit = emit  # type: ignore[assignment]
+        logger.handlers = [handler]
+        logger.setLevel(logging.INFO)
+        self._log_handler = handler
+
+    def _restore_logging(self) -> None:
+        logger = logging.getLogger("reactor_sim")
+        if self._previous_handlers:
+            logger.handlers = self._previous_handlers
+        if self._log_handler and self._log_handler in logger.handlers:
+            logger.removeHandler(self._log_handler)
+
+    def on_button_pressed(self, event: Button.Pressed) -> None:  # type: ignore[override]
+        mapping = {
+            "scram": self.action_scram,
+            "p1": self.action_toggle_primary_p1,
+            "p2": self.action_toggle_primary_p2,
+            "s1": self.action_toggle_secondary_p1,
+            "s2": self.action_toggle_secondary_p2,
+            "relief_pri": self.action_toggle_primary_relief,
+            "relief_sec": self.action_toggle_secondary_relief,
+            "turbines": self.action_toggle_turbine_bank,
+            "t1": lambda: self.action_toggle_turbine_unit("1"),
+            "t2": lambda: self.action_toggle_turbine_unit("2"),
+            "t3": lambda: self.action_toggle_turbine_unit("3"),
+            "g1": self.action_toggle_generator1,
+            "g2": self.action_toggle_generator2,
+            "consumer": self.action_toggle_consumer,
+            "auto_rods": self.action_toggle_auto_rods,
+            "rods_plus": self.action_rods_insert,
+            "rods_minus": self.action_rods_withdraw,
+            "demand_up": self.action_demand_up,
+            "demand_down": self.action_demand_down,
+            "sp_up": self.action_setpoint_up,
+            "sp_down": self.action_setpoint_down,
+            "snapshot": lambda: self._save_snapshot(auto=False),
+            "m_core": self.action_maintain_core,
+            "m_p1": self.action_maintain_primary_p1,
+            "m_p2": self.action_maintain_primary_p2,
+            "m_s1": self.action_maintain_secondary_p1,
+            "m_s2": self.action_maintain_secondary_p2,
+            "m_t1": self.action_maintain_turbine_1,
+            "m_t2": self.action_maintain_turbine_2,
+            "m_t3": self.action_maintain_turbine_3,
+            "m_g1": self.action_maintain_generator_1,
+            "m_g2": self.action_maintain_generator_2,
+        }
+        handler = mapping.get(event.button.id or "")
+        if handler:
+            handler()
+
+
+def run_textual_dashboard(reactor: Reactor, start_state: Optional[PlantState], timestep: float, save_path: Optional[str]) -> None:
+    app = TextualDashboard(reactor, start_state, timestep=timestep, save_path=save_path)
+    app.run()

src/reactor_sim/thermal.py

@@ -13,15 +13,37 @@ from .state import CoolantLoopState, CoreState
 LOGGER = logging.getLogger(__name__)
 
 
-def heat_transfer(primary: CoolantLoopState, secondary: CoolantLoopState, core_power_mw: float) -> float:
+def heat_transfer(
+    primary: CoolantLoopState, secondary: CoolantLoopState, core_power_mw: float, fouling_factor: float = 0.0
+) -> float:
     """Return MW transferred to the secondary loop."""
-    if secondary.mass_flow_rate <= 0.0:
+    if primary.mass_flow_rate <= 0.0 or secondary.mass_flow_rate <= 0.0:
         return 0.0
-    delta_t = max(0.0, primary.temperature_out - secondary.temperature_in)
+    delta_t1 = max(1e-3, primary.temperature_out - secondary.temperature_in)
-    conductance = 0.15  # steam generator effectiveness
+    delta_t2 = max(1e-3, primary.temperature_in - secondary.temperature_out)
-    efficiency = 1.0 - math.exp(-conductance * delta_t)
+    if delta_t1 <= 0.0 or delta_t2 <= 0.0:
-    transferred = min(core_power_mw, core_power_mw * efficiency)
+        return 0.0
-    LOGGER.debug("Heat transfer %.2f MW with ΔT=%.1fK", transferred, delta_t)
+    if abs(delta_t1 - delta_t2) < 1e-6:
+        lmtd = delta_t1
+    else:
+        lmtd = (delta_t1 - delta_t2) / math.log(delta_t1 / delta_t2)
+    fouling = max(0.0, min(constants.HX_FOULING_MAX_PENALTY, fouling_factor))
+    ua = constants.STEAM_GENERATOR_UA_MW_PER_K * (1.0 - fouling)
+    ua_limited = ua * lmtd
+
+    # Prevent the heat exchanger from over-transferring and inverting the outlet temperatures.
+    primary_capacity = primary.mass_flow_rate * constants.COOLANT_HEAT_CAPACITY
+    secondary_capacity = secondary.mass_flow_rate * constants.COOLANT_HEAT_CAPACITY
+    approach = 2.0  # K minimum approach between loop outlets
+    pinch_limited = 0.0
+    if primary_capacity > 0.0 and secondary_capacity > 0.0:
+        temp_gap = primary.temperature_out - secondary.temperature_in - approach
+        if temp_gap > 0.0:
+            pinch_watts = temp_gap / ((1.0 / primary_capacity) + (1.0 / secondary_capacity))
+            pinch_limited = max(0.0, pinch_watts / constants.MEGAWATT)
+
+    transferred = max(0.0, min(core_power_mw, ua_limited, pinch_limited))
+    LOGGER.debug("Heat transfer %.2f MW with LMTD=%.1fK (ΔT1=%.1f ΔT2=%.1f)", transferred, lmtd, delta_t1, delta_t2)
     return transferred
 
 
@@ -31,19 +53,108 @@ def temperature_rise(power_mw: float, mass_flow: float) -> float:
     return (power_mw * constants.MEGAWATT) / (mass_flow * constants.COOLANT_HEAT_CAPACITY)
 
 
+def saturation_pressure(temp_k: float) -> float:
+    """Return approximate saturation pressure (MPa) for water at temp_k."""
+    temp_c = max(0.0, temp_k - 273.15)
+    # Antoine equation parameters for water in 1-374C range, pressure in mmHg.
+    a, b, c = 8.14019, 1810.94, 244.485
+    psat_mmHg = 10 ** (a - (b / (c + temp_c)))
+    psat_mpa = psat_mmHg * 133.322e-6
+    return min(constants.MAX_PRESSURE, max(0.01, psat_mpa))
+
+
+def saturation_temperature(pressure_mpa: float) -> float:
+    """Approximate saturation temperature (K) for water at the given pressure."""
+    target = max(0.01, min(constants.MAX_PRESSURE, pressure_mpa))
+    low, high = 273.15, 900.0
+    for _ in range(40):
+        mid = 0.5 * (low + high)
+        if saturation_pressure(mid) < target:
+            low = mid
+        else:
+            high = mid
+    return high
+
+
 @dataclass
 class ThermalSolver:
     primary_volume_m3: float = 300.0
 
-    def step_core(self, core: CoreState, primary: CoolantLoopState, power_mw: float, dt: float) -> None:
+    def _resolve_secondary_state(self, secondary: CoolantLoopState) -> None:
+        """Project stored energy onto temperature, quality, and pressure."""
+        cp = constants.COOLANT_HEAT_CAPACITY
+        mass = max(1e-6, secondary.inventory_kg)
+        secondary.energy_j = max(0.0, secondary.energy_j)
+        sat_temp = saturation_temperature(max(0.05, secondary.pressure))
+        liquid_energy = mass * cp * sat_temp
+        available = secondary.energy_j
+
+        if available <= liquid_energy:
+            temp = available / (mass * cp)
+            secondary.temperature_out = max(constants.ENVIRONMENT_TEMPERATURE, temp)
+            secondary.steam_quality = 0.0
+        else:
+            latent_energy = min(available - liquid_energy, mass * constants.STEAM_LATENT_HEAT)
+            quality = latent_energy / (mass * constants.STEAM_LATENT_HEAT)
+            superheat_energy = max(0.0, available - liquid_energy - latent_energy)
+            superheat_temp = superheat_energy / (mass * cp) if quality >= 1.0 else 0.0
+            secondary.temperature_out = sat_temp + superheat_temp
+            secondary.steam_quality = max(0.0, min(1.0, quality))
+            secondary.energy_j = liquid_energy + latent_energy + superheat_energy
+
+        secondary.pressure = min(
+            constants.MAX_PRESSURE, max(secondary.pressure, saturation_pressure(secondary.temperature_out))
+        )
+
+    def step_core(
+        self,
+        core: CoreState,
+        primary: CoolantLoopState,
+        power_mw: float,
+        dt: float,
+        residual_power_mw: float | None = None,
+    ) -> None:
+        def _lag(prev: float, new: float, tau: float) -> float:
+            if tau <= 0.0:
+                return new
+            alpha = min(1.0, max(0.0, dt / max(1e-6, tau)))
+            return prev + alpha * (new - prev)
+
+        if residual_power_mw is None:
+            residual_power_mw = power_mw
+        prev_fuel = core.fuel_temperature
+        prev_clad = core.clad_temperature or primary.temperature_out
         temp_rise = temperature_rise(power_mw, primary.mass_flow_rate)
         primary.temperature_out = primary.temperature_in + temp_rise
-        # Fuel heats from any power not immediately convected away, and cools toward the primary outlet.
+        # Fuel heats from total fission power (even when most is convected) plus any residual left in the coolant.
-        heating = 0.005 * max(0.0, power_mw - temp_rise) * dt
+        heating = (0.003 * max(0.0, power_mw) + 0.01 * max(0.0, residual_power_mw)) * dt
-        cooling = 0.05 * max(0.0, core.fuel_temperature - primary.temperature_out) * dt
+        cooling = 0.025 * 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.
+        # Simple radial split: fuel -> clad conduction with burnup-driven conductivity drop and gap penalty.
-        core.fuel_temperature = min(max(primary.temperature_out, core.fuel_temperature), constants.MAX_CORE_TEMPERATURE)
+        gap_penalty = 1.0 + 2.0 * core.burnup
+        conductivity = 0.03 / gap_penalty
+        conduction = conductivity * max(0.0, core.fuel_temperature - (core.clad_temperature or primary.temperature_out)) * dt
+        clad = core.clad_temperature or primary.temperature_out
+        clad_cooling = 0.06 * max(0.0, clad - primary.temperature_out) * dt
+        clad = max(primary.temperature_out, clad + conduction - clad_cooling)
+        core.fuel_temperature = max(primary.temperature_out, core.fuel_temperature - conduction)
+        # Peak pellet centerline surrogate based on power and burnup conductivity loss.
+        peak_delta = 20.0 * (core.power_output_mw / max(1.0, constants.NORMAL_CORE_POWER_MW)) * (1.0 + core.burnup)
+        core.pellet_center_temperature = min(constants.MAX_CORE_TEMPERATURE, core.fuel_temperature + peak_delta)
+        # Keep temperatures bounded and never below coolant outlet.
+        core.fuel_temperature = min(core.fuel_temperature, constants.MAX_CORE_TEMPERATURE)
+        core.clad_temperature = min(clad, constants.MAX_CORE_TEMPERATURE)
+        # Apply mild lags so heat moves from fuel to clad to coolant over a short time constant.
+        core.fuel_temperature = _lag(prev_fuel, core.fuel_temperature, constants.FUEL_TO_CLAD_TIME_CONSTANT)
+        core.clad_temperature = _lag(prev_clad, core.clad_temperature or prev_clad, constants.CLAD_TO_COOLANT_TIME_CONSTANT)
+        core.subcooling_margin = max(0.0, saturation_temperature(primary.pressure) - primary.temperature_out)
+        chf = self._critical_heat_flux(primary)
+        heat_flux = (power_mw * constants.MEGAWATT) / max(1.0, self._core_surface_area())
+        core.dnb_margin = max(0.0, chf / max(1e-6, heat_flux))
+        avg_temp = 0.5 * (primary.temperature_in + primary.temperature_out)
+        primary.energy_j = max(
+            0.0, primary.inventory_kg * constants.COOLANT_HEAT_CAPACITY * avg_temp
+        )
         LOGGER.debug(
             "Primary loop: flow=%.0f kg/s temp_out=%.1fK core_temp=%.1fK",
             primary.mass_flow_rate,
@@ -51,14 +162,50 @@ class ThermalSolver:
             core.fuel_temperature,
         )
 
-    def step_secondary(self, secondary: CoolantLoopState, transferred_mw: float) -> None:
+    def step_secondary(self, secondary: CoolantLoopState, transferred_mw: float, dt: float = 1.0) -> None:
-        delta_t = temperature_rise(transferred_mw, secondary.mass_flow_rate)
+        """Update secondary loop using a stored-energy steam-drum balance."""
-        secondary.temperature_out = secondary.temperature_in + delta_t
+        cp = constants.COOLANT_HEAT_CAPACITY
-        secondary.steam_quality = min(1.0, max(0.0, delta_t / 100.0))
+        mass = max(1e-6, secondary.inventory_kg)
-        secondary.pressure = min(constants.MAX_PRESSURE, 6.0 + delta_t * 0.01)
+        if secondary.energy_j <= 0.0:
+            secondary.energy_j = mass * cp * secondary.average_temperature()
+
+        # Add transferred heat; if no heat, bleed toward ambient.
+        if transferred_mw > 0.0:
+            secondary.energy_j += transferred_mw * constants.MEGAWATT * dt
+        else:
+            bleed = 0.01 * (secondary.temperature_out - constants.ENVIRONMENT_TEMPERATURE)
+            secondary.energy_j = max(
+                mass * cp * constants.ENVIRONMENT_TEMPERATURE, secondary.energy_j - max(0.0, bleed) * mass * cp * dt
+            )
+
+        self._resolve_secondary_state(secondary)
         LOGGER.debug(
-            "Secondary loop: transferred=%.1fMW temp_out=%.1fK quality=%.2f",
+            "Secondary loop: transferred=%.1fMW temp_out=%.1fK quality=%.2f energy=%.1eJ",
             transferred_mw,
             secondary.temperature_out,
             secondary.steam_quality,
+            secondary.energy_j,
         )
+
+    def remove_steam_energy(self, secondary: CoolantLoopState, steam_power_mw: float, dt: float) -> None:
+        """Remove steam enthalpy consumed by turbines and rebalance the drum."""
+        if steam_power_mw <= 0.0:
+            return
+        secondary.energy_j = max(0.0, secondary.energy_j - steam_power_mw * constants.MEGAWATT * dt)
+        self._resolve_secondary_state(secondary)
+
+    def _critical_heat_flux(self, primary: CoolantLoopState) -> float:
+        """CHF surrogate with pressure, mass flux, and subcooling influence (Chen-like)."""
+        area = self._core_surface_area()
+        mass_flux = max(1.0, primary.mass_flow_rate / max(1.0, area))  # kg/s per m2 surrogate
+        pressure_factor = 0.6 + 0.4 * min(1.0, primary.pressure / max(0.1, constants.MAX_PRESSURE))
+        subcool = max(0.0, saturation_temperature(primary.pressure) - primary.temperature_out)
+        subcool_factor = max(0.2, min(1.0, subcool / 30.0))
+        flux_factor = max(0.4, min(3.0, (mass_flux / 500.0) ** 0.5))
+        base_chf = 1.2e7  # W/m2 surrogate baseline
+        return base_chf * flux_factor * pressure_factor * subcool_factor
+
+    def _core_surface_area(self) -> float:
+        # Simple surrogate: area scaling with volume^(2/3)
+        volume = self.primary_volume_m3
+        return max(1.0, (volume ** (2.0 / 3.0)) * 10.0)

src/reactor_sim/turbine.py

@@ -6,6 +6,7 @@ from dataclasses import dataclass
 import logging
 
 from . import constants
+from .thermal import saturation_temperature, saturation_pressure
 from .state import CoolantLoopState, TurbineState
 
 LOGGER = logging.getLogger(__name__)
@@ -27,6 +28,7 @@ class Turbine:
     mechanical_efficiency: float = constants.STEAM_TURBINE_EFFICIENCY
     rated_output_mw: float = 400.0  # cap per unit electrical output
     spool_time: float = constants.TURBINE_SPOOL_TIME
+    throttle: float = 1.0  # 0-1 valve position
 
     def step(
         self,
@@ -35,15 +37,47 @@ class Turbine:
         steam_power_mw: float = 0.0,
         dt: float = 1.0,
     ) -> None:
-        enthalpy = 2_700.0 + loop.steam_quality * 600.0
+        effective_mass_flow = loop.mass_flow_rate * max(0.0, loop.steam_quality)
-        mass_flow = loop.mass_flow_rate * 0.6
+        if steam_power_mw <= 0.0 and (loop.steam_quality <= 0.01 or effective_mass_flow <= 10.0):
-        available_power = max(steam_power_mw, (enthalpy * mass_flow / 1_000.0) / 1_000.0)
+            # No steam available; turbine should idle.
-        shaft_power_mw = available_power * self.mechanical_efficiency
+            state.shaft_power_mw = 0.0
+            state.electrical_output_mw = 0.0
+            state.load_demand_mw = 0.0
+            state.load_supplied_mw = 0.0
+            state.steam_enthalpy = 0.0
+            state.condenser_temperature = max(constants.CONDENSER_COOLING_WATER_TEMP_K, loop.temperature_in - 20.0)
+            state.condenser_pressure = max(constants.CONDENSER_BASE_PRESSURE_MPA, state.condenser_pressure - 0.01 * dt)
+            state.fouling_penalty = max(0.0, state.fouling_penalty - 0.0001 * dt)
+            return
+
+        throttle = min(constants.TURBINE_THROTTLE_MAX, max(constants.TURBINE_THROTTLE_MIN, self.throttle))
+        throttle_eff = 1.0 - constants.TURBINE_THROTTLE_EFFICIENCY_DROP * (constants.TURBINE_THROTTLE_MAX - throttle)
+
+        sat_temp = saturation_temperature(max(0.05, loop.pressure))
+        superheat = max(0.0, loop.temperature_out - sat_temp)
+        enthalpy = (constants.STEAM_LATENT_HEAT / 1_000.0) + (superheat * constants.COOLANT_HEAT_CAPACITY / 1_000.0)
+        mass_flow = effective_mass_flow * 0.6 * throttle
+        computed_power = (enthalpy * mass_flow) / 1_000.0  # MW from enthalpy flow
+        available_power = steam_power_mw if steam_power_mw > 0 else computed_power
+        available_power = min(available_power, computed_power)
+        backpressure_loss = 1.0 - _backpressure_penalty(state)
+        shaft_power_mw = available_power * self.mechanical_efficiency * throttle_eff * backpressure_loss
         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)
+        condenser_temp = max(constants.CONDENSER_COOLING_WATER_TEMP_K, loop.temperature_in - 20.0)
+        # Vacuum pump tends toward base pressure; fouling raises it slowly when hot.
+        target_pressure = constants.CONDENSER_BASE_PRESSURE_MPA
+        if condenser_temp > constants.CONDENSER_COOLING_WATER_TEMP_K + 20.0:
+            state.fouling_penalty = min(
+                constants.CONDENSER_FOULING_MAX_PENALTY,
+                state.fouling_penalty + constants.CONDENSER_FOULING_RATE * dt,
+            )
+        state.condenser_pressure = max(
+            target_pressure,
+            min(constants.CONDENSER_MAX_PRESSURE_MPA, state.condenser_pressure - constants.CONDENSER_VACUUM_PUMP_RATE * dt),
+        )
         state.steam_enthalpy = enthalpy
         state.shaft_power_mw = _ramp(state.shaft_power_mw, shaft_power_mw, dt, self.spool_time)
         state.electrical_output_mw = _ramp(state.electrical_output_mw, electrical, dt, self.spool_time)
@@ -59,5 +93,16 @@ class Turbine:
 def _ramp(current: float, target: float, dt: float, time_constant: float) -> float:
     if time_constant <= 0.0:
         return target
-    alpha = min(1.0, max(0.0, dt / time_constant))
+    alpha = min(1.0, max(0.0, dt / max(1e-6, time_constant)))
     return current + (target - current) * alpha
+
+
+def _backpressure_penalty(state: TurbineState) -> float:
+    base = constants.CONDENSER_BASE_PRESSURE_MPA
+    max_p = constants.CONDENSER_MAX_PRESSURE_MPA
+    pressure = max(base, min(max_p, state.condenser_pressure))
+    if pressure <= base:
+        return min(constants.CONDENSER_BACKPRESSURE_PENALTY, state.fouling_penalty)
+    frac = (pressure - base) / max(1e-6, max_p - base)
+    penalty = frac * constants.CONDENSER_BACKPRESSURE_PENALTY
+    return min(constants.CONDENSER_BACKPRESSURE_PENALTY + state.fouling_penalty, penalty + state.fouling_penalty)

tests/test_aux_power.py

@@ -0,0 +1,18 @@
+from reactor_sim.reactor import Reactor
+from reactor_sim.commands import ReactorCommand
+
+
+def test_pump_stays_off_without_aux_power():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = False
+    reactor.primary_pump_active = True
+    reactor.primary_pump_units = [True, False]
+    reactor.secondary_pump_active = False
+    reactor.generator_auto = False
+    reactor.turbine_unit_active = [False, False, False]
+
+    reactor.step(state, dt=1.0, command=ReactorCommand(primary_pumps={1: True}))
+
+    assert state.primary_loop.mass_flow_rate == 0.0
+    assert state.primary_pumps[0].status in {"OFF", "STOPPING"}

tests/test_control.py

@@ -0,0 +1,40 @@
+import pytest
+
+from reactor_sim.control import ControlSystem
+from reactor_sim import constants
+from reactor_sim.state import CoreState
+
+
+def _core_state() -> CoreState:
+    return CoreState(
+        fuel_temperature=300.0,
+        neutron_flux=1e5,
+        reactivity_margin=0.0,
+        power_output_mw=0.0,
+        burnup=0.0,
+    )
+
+
+def test_manual_rods_quantized_to_step():
+    control = ControlSystem()
+    control.manual_control = True
+    core = _core_state()
+
+    control.set_rods(0.333)
+    assert control.rod_target == 0.325
+    control.update_rods(core, dt=100.0)
+    assert control.rod_fraction == pytest.approx(0.325, rel=1e-6)
+
+    control.increment_rods(0.014)
+    assert control.rod_target == pytest.approx(0.35)
+    control.update_rods(core, dt=100.0)
+    assert control.rod_fraction == pytest.approx(0.35, rel=1e-6)
+
+    # Clamp upper bound
+    control.set_rods(1.0)
+    control.update_rods(core, dt=100.0)
+    assert control.rod_fraction == 0.95
+
+
+def test_dashboard_step_constant_exposed():
+    assert constants.ROD_MANUAL_STEP == 0.025

tests/test_neutronics.py

@@ -1,3 +1,5 @@
+import pytest
+
 from reactor_sim.neutronics import NeutronDynamics
 from reactor_sim.state import CoreState
 
@@ -23,3 +25,33 @@ def test_reactivity_increases_with_rod_withdrawal():
     rho_full_out = dynamics.reactivity(state, control_fraction=0.0)
     rho_half = dynamics.reactivity(state, control_fraction=0.5)
     assert rho_full_out > rho_half
+
+
+def test_poisons_accumulate_under_power():
+    dynamics = NeutronDynamics()
+    state = _core_state(power=800.0, flux=1e6)
+    dynamics.update_poisons(state, dt=100.0)
+    dynamics.update_poisons(state, dt=100.0)
+    assert state.iodine_inventory > 0.0
+    assert state.xenon_inventory > 0.0
+
+
+def test_xenon_penalty_caps():
+    dynamics = NeutronDynamics()
+    state = _core_state()
+    state.xenon_inventory = 50.0
+    assert dynamics.xenon_penalty(state) == 0.05
+    state.xenon_inventory = 0.2
+    assert dynamics.xenon_penalty(state) == pytest.approx(0.01)
+
+
+def test_delayed_precursors_follow_rod_banks():
+    dynamics = NeutronDynamics()
+    state = _core_state(power=600.0, flux=5e6)
+    dynamics.step(state, control_fraction=0.2, dt=1.0, rod_banks=[0.2, 0.2, 0.2])
+    initial_sum = sum(state.delayed_precursors)
+    assert len(state.delayed_precursors) == 3
+    assert initial_sum > 0.0
+
+    dynamics.step(state, control_fraction=0.95, dt=2.0, rod_banks=[0.95, 0.95, 0.95])
+    assert sum(state.delayed_precursors) < initial_sum

tests/test_pressurizer.py

@@ -0,0 +1,36 @@
+from reactor_sim.reactor import Reactor
+from reactor_sim.commands import ReactorCommand
+
+
+def test_pressurizer_raises_pressure_with_level():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    state.primary_loop.pressure = 0.5
+    reactor.pressurizer_level = 0.8
+    reactor.primary_pump_active = False
+    reactor.secondary_pump_active = False
+    reactor.shutdown = False
+
+    reactor.step(state, dt=1.0, command=ReactorCommand())
+
+    assert state.primary_loop.pressure > 0.5
+    assert reactor.pressurizer_level < 0.8
+
+
+def test_low_npsh_limits_primary_flow():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = False
+    reactor.control.manual_control = True
+    reactor.control.rod_fraction = 0.0
+    reactor.primary_pump_active = True
+    reactor.primary_pump_units = [True, True]
+    reactor.secondary_pump_active = False
+    state.primary_loop.pressure = 0.05  # near-vacuum to force cavitation
+    state.primary_loop.temperature_in = 400.0
+    state.primary_loop.temperature_out = 600.0
+
+    reactor.step(state, dt=1.0, command=ReactorCommand(generator_units={1: True}))
+
+    assert state.primary_pumps[0].status == "CAV"
+    assert state.primary_loop.mass_flow_rate < 100.0

tests/test_simulation.py

@@ -127,7 +127,10 @@ def test_secondary_pump_unit_toggle_can_restart_pump():
     reactor.secondary_pump_active = False
     reactor.secondary_pump_units = [False, False]
 
-    reactor.step(state, dt=1.0, command=ReactorCommand(secondary_pumps={1: True}))
+    cmd = ReactorCommand(secondary_pumps={1: True}, generator_units={1: True})
+    for _ in range(5):
+        reactor.step(state, dt=1.0, command=cmd)
+        cmd = ReactorCommand(coolant_demand=0.75)
 
     assert reactor.secondary_pump_units == [True, False]
     assert reactor.secondary_pump_active is True
@@ -140,11 +143,9 @@ def test_primary_pumps_spool_up_over_seconds():
     reactor.secondary_pump_units = [False, False]
     # Enable both pumps and command full flow; spool should take multiple steps.
     target_flow = reactor.primary_pump.flow_rate(1.0) * len(reactor.primary_pump_units)
-    reactor.step(
+    cmd = ReactorCommand(primary_pumps={1: True, 2: True}, generator_units={1: True}, coolant_demand=1.0)
-        state,
+    reactor.step(state, dt=1.0, command=cmd)
-        dt=1.0,
+    reactor.step(state, dt=1.0, command=ReactorCommand(coolant_demand=1.0))
-        command=ReactorCommand(primary_pumps={1: True, 2: True}, generator_units={1: True}, coolant_demand=1.0),
-    )
     first_flow = state.primary_loop.mass_flow_rate
     assert 0.0 < first_flow < target_flow
 
@@ -164,6 +165,8 @@ def test_full_rod_withdrawal_reaches_gigawatt_power():
     reactor.secondary_pump_active = True
     reactor.primary_pump_units = [True, True]
     reactor.secondary_pump_units = [True, True]
+    reactor.generator_auto = True
+    reactor.step(state, dt=1.0, command=ReactorCommand(generator_units={1: True}))
 
     early_power = 0.0
     for step in range(60):
@@ -178,18 +181,19 @@ def test_partially_inserted_rods_hold_near_three_gw():
     reactor = Reactor.default()
     state = reactor.initial_state()
     reactor.shutdown = False
-    reactor.control.manual_control = True
+    reactor.control.set_manual_mode(False)
-    reactor.control.rod_fraction = 0.4
     reactor.primary_pump_active = True
     reactor.secondary_pump_active = True
     reactor.primary_pump_units = [True, True]
     reactor.secondary_pump_units = [True, True]
+    reactor.generator_auto = True
+    reactor.step(state, dt=1.0, command=ReactorCommand(generator_units={1: True}))
 
-    for _ in range(120):
+    for _ in range(180):
         reactor.step(state, dt=1.0)
 
-    assert 2_000.0 < state.core.power_output_mw < 4_000.0
+    assert 2_500.0 < state.core.power_output_mw < 3_500.0
-    assert 500.0 < state.core.fuel_temperature < 800.0
+    assert 0.05 < reactor.control.rod_fraction < 0.9
 
 
 def test_generator_spools_and_powers_pumps():
@@ -248,3 +252,127 @@ def test_meltdown_triggers_shutdown():
 
     assert reactor.shutdown is True
     assert reactor.meltdown is True
+
+
+def test_auto_control_resets_shutdown_and_moves_rods():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = True
+    reactor.control.manual_control = True
+    reactor.control.rod_fraction = 0.95
+
+    reactor.step(state, dt=1.0, command=ReactorCommand(rod_manual=False))
+
+    assert reactor.shutdown is False
+    assert reactor.control.manual_control is False
+    assert reactor.control.rod_fraction < 0.95
+
+
+def test_chemistry_builds_fouling_and_backpressure():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    # Push impurities high to accelerate fouling dynamics.
+    state.dissolved_oxygen_ppm = 200.0
+    state.sodium_ppm = 100.0
+    state.secondary_loop.mass_flow_rate = 20_000.0
+    state.secondary_loop.steam_quality = 0.3
+    state.secondary_loop.temperature_out = 600.0
+    state.secondary_loop.temperature_in = 560.0
+    base_hx = state.hx_fouling
+    base_foul = state.turbines[0].fouling_penalty if state.turbines else 0.0
+    base_pressure = state.turbines[0].condenser_pressure if state.turbines else constants.CONDENSER_BASE_PRESSURE_MPA
+
+    reactor._update_chemistry(state, dt=20.0)
+
+    assert state.hx_fouling > base_hx
+    if state.turbines:
+        assert state.turbines[0].fouling_penalty > base_foul
+        assert state.turbines[0].condenser_pressure >= base_pressure
+
+
+def test_full_power_reaches_steam_and_turbine_output():
+    """Integration: ramp to full power with staged rod control and verify sustained steam/electric output."""
+    reactor = Reactor.default()
+    reactor.health_monitor.disable_degradation = True
+    reactor.allow_external_aux = True
+    reactor.relaxed_npsh = True
+    reactor.control.set_power_setpoint(3_000.0)
+    state = reactor.initial_state()
+    reactor.step(
+        state,
+        dt=1.0,
+        command=ReactorCommand(
+            generator_units={1: True, 2: True},
+            primary_pumps={1: True, 2: True},
+            secondary_pumps={1: True, 2: True},
+            rod_manual=True,
+            rod_position=0.55,
+        ),
+    )
+    checkpoints = {300, 600, 900, 1800, 2700, 3600}
+    results = {}
+    turbines_started = False
+    for i in range(3600):
+        cmd = None
+        if state.core.power_output_mw >= 2_500.0 and reactor.control.manual_control:
+            cmd = ReactorCommand(rod_manual=False)
+        if (
+            not turbines_started
+            and state.secondary_loop.steam_quality > 0.02
+            and state.secondary_loop.pressure > 1.0
+        ):
+            cmd = ReactorCommand(turbine_on=True, turbine_units={1: True, 2: True, 3: True})
+            turbines_started = True
+        if i == 600 and not turbines_started:
+            cmd = ReactorCommand(turbine_on=True, turbine_units={1: True, 2: True, 3: True})
+            turbines_started = True
+        reactor.step(state, dt=1.0, command=cmd)
+        if state.time_elapsed in checkpoints:
+            results[state.time_elapsed] = {
+                "quality": state.secondary_loop.steam_quality,
+                "electric": state.total_electrical_output(),
+                "core_temp": state.core.fuel_temperature,
+            }
+
+    # At or after 10 minutes of operation, ensure we have meaningful steam and electrical output.
+    assert results[600]["quality"] > 0.05
+    assert results[600]["electric"] > 50.0
+    assert results[3600]["quality"] > 0.1
+    assert results[3600]["electric"] > 150.0
+    # No runaway core temperatures.
+    assert results[3600]["core_temp"] < constants.CORE_MELTDOWN_TEMPERATURE
+
+
+def test_cooldown_reaches_ambient_without_runaway():
+    """Shutdown with pumps running should cool loops toward ambient, no runaway."""
+    reactor = Reactor.default()
+    reactor.health_monitor.disable_degradation = True
+    reactor.allow_external_aux = True
+    reactor.relaxed_npsh = True
+    state = reactor.initial_state()
+    # Start hot.
+    reactor.step(
+        state,
+        dt=1.0,
+        command=ReactorCommand(
+            generator_units={1: True, 2: True},
+            primary_pumps={1: True, 2: True},
+            secondary_pumps={1: True, 2: True},
+            rod_manual=True,
+            rod_position=0.55,
+        ),
+    )
+    turbines_started = False
+    for i in range(1800):
+        cmd = None
+        if not turbines_started and state.secondary_loop.steam_quality > 0.02 and state.secondary_loop.pressure > 1.0:
+            cmd = ReactorCommand(turbine_on=True, turbine_units={1: True, 2: True, 3: True})
+            turbines_started = True
+        if i == 900:
+            cmd = ReactorCommand(rod_position=0.95, turbine_on=False, turbine_units={1: False, 2: False, 3: False})
+        reactor.step(state, dt=1.0, command=cmd)
+
+    assert not reactor.meltdown
+    assert state.core.power_output_mw < 1.0
+    assert state.primary_loop.temperature_out < 320.0
+    assert state.secondary_loop.temperature_out < 320.0

tests/test_thermal.py

@@ -0,0 +1,50 @@
+import pytest
+
+from reactor_sim import constants
+from reactor_sim.state import CoolantLoopState
+from reactor_sim.thermal import ThermalSolver, heat_transfer, saturation_temperature
+
+
+def _secondary_loop(temp_in: float = 350.0, pressure: float = 0.5, flow: float = 200.0) -> CoolantLoopState:
+    nominal_mass = constants.SECONDARY_LOOP_VOLUME_M3 * constants.COOLANT_DENSITY
+    return CoolantLoopState(
+        temperature_in=temp_in,
+        temperature_out=temp_in,
+        pressure=pressure,
+        mass_flow_rate=flow,
+        steam_quality=0.0,
+        inventory_kg=nominal_mass,
+        level=constants.SECONDARY_INVENTORY_TARGET,
+    )
+
+
+def test_secondary_heats_to_saturation_before_boiling():
+    solver = ThermalSolver()
+    loop = _secondary_loop(temp_in=330.0, flow=180.0, pressure=0.5)
+    sat_temp = saturation_temperature(loop.pressure)
+    solver.step_secondary(loop, transferred_mw=50.0, dt=1.0)
+    assert loop.steam_quality == pytest.approx(0.0)
+    assert loop.temperature_out < sat_temp
+
+
+def test_secondary_generates_steam_when_energy_exceeds_sensible_heat():
+    solver = ThermalSolver()
+    loop = _secondary_loop(temp_in=330.0, flow=180.0, pressure=0.5)
+    loop.inventory_kg *= 0.1  # reduce mass to let boil-up happen quickly
+    sat_temp = saturation_temperature(loop.pressure)
+    solver.step_secondary(loop, transferred_mw=120.0, dt=100.0)
+    assert loop.temperature_out == pytest.approx(sat_temp, rel=0.05)
+    assert loop.steam_quality > 0.0
+    assert loop.steam_quality < 1.0
+
+
+def test_heat_transfer_reduced_by_fouling():
+    primary = CoolantLoopState(
+        temperature_in=360.0, temperature_out=380.0, pressure=15.0, mass_flow_rate=50_000.0, steam_quality=0.0
+    )
+    secondary = CoolantLoopState(
+        temperature_in=320.0, temperature_out=330.0, pressure=6.5, mass_flow_rate=50_000.0, steam_quality=0.1
+    )
+    clean = heat_transfer(primary, secondary, core_power_mw=7_000.0, fouling_factor=0.0)
+    fouled = heat_transfer(primary, secondary, core_power_mw=7_000.0, fouling_factor=0.25)
+    assert fouled < clean

tests/test_turbine.py

@@ -6,10 +6,11 @@ from reactor_sim.turbine import Turbine
 
 def test_turbine_spools_toward_target_output():
     turbine = Turbine()
+    turbine.throttle = 1.0
     loop = CoolantLoopState(
         temperature_in=600.0,
         temperature_out=650.0,
-        pressure=6.0,
+        pressure=0.02,
         mass_flow_rate=20_000.0,
         steam_quality=0.9,
     )
@@ -20,14 +21,37 @@ def test_turbine_spools_toward_target_output():
         condenser_temperature=300.0,
     )
     target_electric = min(
-        turbine.rated_output_mw, 300.0 * turbine.mechanical_efficiency * turbine.generator_efficiency
+        turbine.rated_output_mw,
+        300.0 * turbine.mechanical_efficiency * turbine.generator_efficiency,
     )
 
-    dt = 5.0
+    dt = 1.0
     turbine.step(loop, state, steam_power_mw=300.0, dt=dt)
     assert 0.0 < state.electrical_output_mw < target_electric
 
-    for _ in range(5):
+    for _ in range(60):
         turbine.step(loop, state, steam_power_mw=300.0, dt=dt)
 
     assert state.electrical_output_mw == pytest.approx(target_electric, rel=0.05)
+
+
+def test_turbine_produces_no_power_without_steam():
+    turbine = Turbine()
+    loop = CoolantLoopState(
+        temperature_in=295.0,
+        temperature_out=295.0,
+        pressure=6.0,
+        mass_flow_rate=20_000.0,
+        steam_quality=0.0,
+    )
+    state = TurbineState(
+        steam_enthalpy=0.0,
+        shaft_power_mw=0.0,
+        electrical_output_mw=0.0,
+        condenser_temperature=300.0,
+    )
+
+    turbine.step(loop, state, steam_power_mw=0.0, dt=1.0)
+
+    assert state.electrical_output_mw == 0.0
+    assert state.shaft_power_mw == 0.0