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

CONTEXT_NOTES.md

@@ -1,21 +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.
-- Turbines: produce zero output unless steam quality is present and effective steam flow is >10 kg/s. Steam pressure shown on dashboard only when quality ≥0.05 and flow ≥100 kg/s; otherwise 0 MPa. Steam supply displayed in Turbine panel.
+- 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. Steam supply pressure shown in turbine panel. 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.
+- 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) clamp loop pressure to saturation when open; status displayed per loop.
+- 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: require meaningful steam flow/quality; otherwise zero output. Steam supply pressure in turbine panel reads 0 when no steam.
+- 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

@@ -1,10 +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; manual mode with staged bank motion and SCRAM; state persists across runs.
+- **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.
+- **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, condenser back-pressure penalty, load dispatch to consumer, steam quality gating for output, generator states with batteries/spool.
-- **Protections & failures**: health monitor degrading components under stress, automatic SCRAM on core or heat-sink loss, relief valves per loop, maintenance actions to restore integrity.
+- **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

@@ -4,11 +4,21 @@
 - [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.
-- [ ] Introduce CHF/DNB margin, clad/fuel split temps, and SCRAM matrix for subcooling loss or SG level/pressure trips.
-- [ ] Flesh out condenser behavior: vacuum pump limits, cooling water temperature coupling, and dynamic back-pressure with fouling.
-- [ ] Dashboard polish: compact turbine/generator rows, color critical warnings (SCRAM/heat-sink), and reduce repeated log noise.
+- [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:
-  - 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.
-  - 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.
-  - 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.
-  - 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] 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/constants.py

@@ -7,11 +7,11 @@ 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
@@ -26,8 +26,8 @@ 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 = 8.0  # approximate shutoff head for primary pumps
-SECONDARY_PUMP_SHUTOFF_HEAD_MPA = 3.0
+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
@@ -37,6 +37,12 @@ TURBINE_THROTTLE_EFFICIENCY_DROP = 0.15  # efficiency loss when at minimum throt
 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
@@ -45,14 +51,14 @@ 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 = 7.0  # MPa typical RBMK channel header pressure
-SECONDARY_NOMINAL_PRESSURE = 7.0  # MPa steam drum/steam line pressure surrogate
-STEAM_GENERATOR_UA_MW_PER_K = 25.0  # overall UA for steam generator (MW/K)
+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 = 7.0
-PRIMARY_PRESSURIZER_DEADBAND_MPA = 0.15
+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
@@ -63,6 +69,26 @@ 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

@@ -24,6 +24,11 @@ class ControlSystem:
     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):
@@ -34,11 +39,68 @@ class ControlSystem:
             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.
-        adjustment = error * 0.2
-        adjustment = clamp(adjustment, -constants.CONTROL_ROD_SPEED * dt, constants.CONTROL_ROD_SPEED * dt)
+
+        # On transition from manual -> auto, start a gentle setpoint ramp from current power.
+        if self._last_manual:
+            self._ramp_start_mw = max(0.0, state.power_output_mw)
+            self._ramp_progress_mw = 0.0
+        self._last_manual = False
+
+        # 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)
@@ -70,6 +132,22 @@ class ControlSystem:
             self.manual_control = manual
         LOGGER.info("Rod control %s", "manual" if manual else "automatic")
 
+    def safety_backoff(self, subcooling_margin: float | None, dnb_margin: float | None, dt: float) -> None:
+        """Insert rods proactively when thermal margins are thin."""
+        if self.manual_control:
+            return
+        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,
@@ -84,12 +162,13 @@ class ControlSystem:
         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.15 + 0.2 * power_fraction
-        # Allow warmer operation when electrical load is already being served (turbines online),
-        # but keep a higher floor when idling so test scenarios still converge near 3 GW.
-        if electrical_output_mw is not None and electrical_output_mw > 10.0:
-            power_floor *= 0.6
+            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",

src/reactor_sim/coolant.py

@@ -26,8 +26,9 @@ class Pump:
         """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 = min(1.2, flow / max(1e-3, self.nominal_flow))
-        head = max(0.0, self.shutoff_head_mpa * max(0.0, 1.0 - flow_frac**2))
+        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:

src/reactor_sim/dashboard.py

@@ -13,7 +13,7 @@ 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__)
 
@@ -76,6 +76,7 @@ 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._trend_history: deque[tuple[float, float, float]] = deque(maxlen=120)
         self.log_buffer: deque[str] = deque(maxlen=8)
@@ -90,6 +91,7 @@ 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"),
@@ -213,15 +215,19 @@ class ReactorDashboard:
                 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"):
                 target = (ch - ord("0")) / 10.0
                 self._queue_command(ReactorCommand(rod_position=target, rod_manual=True))
-            elif ch in (ord("!"), ord("@"), ord("#")):
-                idx = ch - ord("!")
-                self._toggle_turbine_unit(idx)
             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:
@@ -419,13 +425,15 @@ class ReactorDashboard:
                 ("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(),
+            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(
@@ -458,7 +466,11 @@ class ReactorDashboard:
                     "Flow",
                     f"{state.secondary_loop.mass_flow_rate:7.0f}/{self.reactor.secondary_pump.nominal_flow * len(self.reactor.secondary_pump_units):.0f} kg/s",
                 ),
-                ("Level", f"{state.secondary_loop.level*100:6.1f}%"),
+                ("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"),
@@ -480,25 +492,37 @@ class ReactorDashboard:
             [
                 ("Turbines", " ".join(self._turbine_status_lines())),
                 ("Rated Elec", f"{len(self.reactor.turbines)*self.reactor.turbines[0].rated_output_mw:7.1f} MW"),
-                ("Steam P", f"{self._steam_pressure(state):5.2f} MPa"),
-                ("Unit1 Elec", f"{state.turbines[0].electrical_output_mw:7.1f} MW" if state.turbines else "n/a"),
-                (
-                    "Unit2 Elec",
-                    f"{state.turbines[1].electrical_output_mw:7.1f} MW" if len(state.turbines) > 1 else "n/a",
-                ),
-                (
-                    "Unit3 Elec",
-                    f"{state.turbines[2].electrical_output_mw:7.1f} MW" if len(state.turbines) > 2 else "n/a",
-                ),
-                ("Throttle", f"{self.reactor.turbines[0].throttle:5.2f}" if self.reactor.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"),
-                ("Consumer", f"{consumer_status}"),
-                ("Demand", f"{consumer_demand:7.1f} MW"),
                 (
                     "Steam",
-                    f"P={state.secondary_loop.pressure:4.2f} MPa q={state.secondary_loop.steam_quality:4.2f} mdot={state.secondary_loop.mass_flow_rate:6.0f} kg/s",
+                    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",
                 ),
+                (
+                    "Units Elec",
+                    " ".join([f"{t.electrical_output_mw:6.1f}MW" for t in state.turbines]) if state.turbines else "n/a",
+                ),
+                (
+                    "Governor",
+                    (
+                        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))
@@ -555,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:
@@ -580,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
@@ -591,15 +615,43 @@ class ReactorDashboard:
         for line in lines:
             if row >= height - 1:
                 break
+            attr = 0
             if isinstance(line, tuple):
+                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"]
@@ -686,57 +738,104 @@ class ReactorDashboard:
     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]] = []
-        lines.append(("SCRAM", "ACTIVE" if self.reactor.shutdown else "CLEAR"))
+        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"))
+            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"
-        lines.append(("Heat sink", heat_text))
+            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"
-        lines.append(("Aux power", aux_text))
+            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")
-        lines.append(("Relief valves", ", ".join(reliefs) if reliefs else "Closed"))
+        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_pressure(self, state: PlantState) -> float:
-        # Only report steam pressure if quality/flow indicate steam is present.
-        if state.secondary_loop.steam_quality < 0.05 or state.secondary_loop.mass_flow_rate < 100.0:
+    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
-        return state.secondary_loop.pressure
+        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) -> list[tuple[str, str]]:
+    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")]
+            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)
@@ -744,9 +843,11 @@ class ReactorDashboard:
         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:
@@ -826,7 +927,13 @@ class _DashboardLogHandler(logging.Handler):
         msg = self.format(record)
         if msg == self._last_msg:
             self._repeat_count += 1
-            if self._repeat_count > 3:
+            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

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")

src/reactor_sim/neutronics.py

@@ -25,10 +25,11 @@ def xenon_poisoning(flux: float) -> float:
 
 @dataclass
 class NeutronDynamics:
+    base_shutdown_bias: float = -0.014
+    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
-    shutdown_bias: float = -0.014
     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

src/reactor_sim/reactor.py

@@ -37,6 +37,8 @@ class Reactor:
     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])
@@ -59,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:
@@ -192,6 +199,9 @@ class Reactor:
             state.core
         )
         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(
@@ -199,7 +209,7 @@ class Reactor:
         )
         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)
+        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.
@@ -216,9 +226,6 @@ class Reactor:
         self._update_loop_inventory(
             state.primary_loop, constants.PRIMARY_LOOP_VOLUME_M3, constants.PRIMARY_INVENTORY_TARGET, dt
         )
-        self._update_loop_inventory(
-            state.secondary_loop, constants.SECONDARY_LOOP_VOLUME_M3, constants.SECONDARY_INVENTORY_TARGET, dt
-        )
 
         pump_demand = overrides.get(
             "coolant_demand",
@@ -255,6 +262,8 @@ class Reactor:
         turbine_electrical = state.total_electrical_output()
         generator_power = self._step_generators(state, aux_demand, turbine_electrical, dt)
         aux_available = turbine_electrical + generator_power
+        if self.allow_external_aux:
+            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:
@@ -273,8 +282,12 @@ class Reactor:
             total_flow = 0.0
             base_flow, base_head = self.primary_pump.performance(pump_demand)
             target_flow = base_flow * power_ratio
-            loop_pressure = max(0.1, saturation_pressure(state.primary_loop.temperature_out))
+            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)
             )
@@ -329,7 +342,11 @@ class Reactor:
             demand = 0.75
             base_flow, base_head = self.secondary_pump.performance(demand)
             target_pressure = max(0.5, base_head * power_ratio)
-            loop_pressure = max(0.1, saturation_pressure(state.secondary_loop.temperature_out))
+            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)
@@ -381,33 +398,74 @@ class Reactor:
                 pump_state.status = "STOPPING" if pump_state.flow_rate > 0.1 else "OFF"
 
         self._apply_pressurizer(state.primary_loop, dt)
-        if self.primary_relief_open:
-            state.primary_loop.pressure = max(0.1, saturation_pressure(state.primary_loop.temperature_out))
-        if self.secondary_relief_open:
-            state.secondary_loop.pressure = max(0.1, saturation_pressure(state.secondary_loop.temperature_out))
 
         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)
-        self.thermal.step_core(state.core, state.primary_loop, total_power, dt)
-        self.thermal.step_secondary(state.secondary_loop, transferred, dt)
-        self._apply_secondary_boiloff(state, dt)
-        self._update_loop_inventory(
-            state.secondary_loop, constants.SECONDARY_LOOP_VOLUME_M3, constants.SECONDARY_INVENTORY_TARGET, dt
+            transferred = heat_transfer(
+                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)
-        if state.core.dnb_margin is not None and state.core.dnb_margin < 0.3:
+        # 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()
-        if state.core.subcooling_margin is not None and state.core.subcooling_margin < 5.0:
+        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,
@@ -426,12 +484,20 @@ class Reactor:
         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 = 0.02 * max(0.0, state.primary_loop.temperature_out - env) * dt
+            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:
-            target_temp = env
+            # 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
             )
@@ -441,6 +507,10 @@ class Reactor:
             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
@@ -472,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)
         ]
@@ -486,19 +557,30 @@ class Reactor:
             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
-                throttle = 0.4 if demand <= 0 else min(1.0, 0.4 + demand / max(1e-6, turbine.rated_output_mw))
-                turbine.throttle = throttle
+                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.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.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
@@ -563,6 +645,83 @@ class Reactor:
         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
@@ -571,11 +730,13 @@ class Reactor:
         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.0
-        return max(0.0, min(1.0, available / constants.NPSH_REQUIRED_MPA))
+            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:
@@ -775,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)):

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,7 +120,21 @@ 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:
+        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(
@@ -128,7 +146,17 @@ def main() -> None:
             dashboard.run()
             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

@@ -21,6 +21,8 @@ class CoreState:
     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
@@ -32,6 +34,8 @@ class CoreState:
     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:
@@ -71,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"
@@ -96,6 +102,10 @@ class PlantState:
     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]:
@@ -129,6 +139,8 @@ class PlantState:
         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.
@@ -142,6 +154,10 @@ class PlantState:
         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(**_with_energy(data["primary_loop"])),
@@ -153,6 +169,10 @@ class PlantState:
             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),
         )
 

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,7 +13,9 @@ 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 primary.mass_flow_rate <= 0.0 or secondary.mass_flow_rate <= 0.0:
         return 0.0
@@ -25,7 +27,8 @@ def heat_transfer(primary: CoolantLoopState, secondary: CoolantLoopState, core_p
         lmtd = delta_t1
     else:
         lmtd = (delta_t1 - delta_t2) / math.log(delta_t1 / delta_t2)
-    ua = constants.STEAM_GENERATOR_UA_MW_PER_K
+    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.
@@ -77,6 +80,32 @@ def saturation_temperature(pressure_mpa: float) -> float:
 class ThermalSolver:
     primary_volume_m3: float = 300.0
 
+    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,
@@ -85,17 +114,39 @@ class ThermalSolver:
         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.
-        heating = 0.005 * max(0.0, residual_power_mw) * dt
+        # Fuel heats from total fission power (even when most is convected) plus any residual left in the coolant.
+        heating = (0.003 * max(0.0, power_mw) + 0.01 * max(0.0, residual_power_mw)) * 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.
-        core.fuel_temperature = min(max(primary.temperature_out, core.fuel_temperature), constants.MAX_CORE_TEMPERATURE)
-        core.clad_temperature = max(primary.temperature_out, core.clad_temperature or primary.temperature_out)
+        # Simple radial split: fuel -> clad conduction with burnup-driven conductivity drop and gap penalty.
+        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())
@@ -127,28 +178,7 @@ class ThermalSolver:
                 mass * cp * constants.ENVIRONMENT_TEMPERATURE, secondary.energy_j - max(0.0, bleed) * mass * cp * dt
             )
 
-        sat_temp = saturation_temperature(max(0.05, secondary.pressure))
-        liquid_energy = mass * cp * sat_temp
-        available = secondary.energy_j
-
-        if available <= liquid_energy:
-            # Subcooled or saturated liquid.
-            temp = available / (mass * cp)
-            secondary.temperature_out = max(temp, constants.ENVIRONMENT_TEMPERATURE)
-            secondary.steam_quality = max(0.0, secondary.steam_quality - 0.01 * dt)
-        else:
-            excess = available - liquid_energy
-            quality = min(1.0, excess / (mass * constants.STEAM_LATENT_HEAT))
-            superheat_energy = max(0.0, excess - quality * mass * constants.STEAM_LATENT_HEAT)
-            superheat_temp = superheat_energy / (mass * cp) if quality >= 1.0 else 0.0
-            secondary.temperature_out = sat_temp + superheat_temp
-            secondary.steam_quality = quality
-            # Re-normalize stored energy to the realized state.
-            secondary.energy_j = liquid_energy + quality * mass * constants.STEAM_LATENT_HEAT + superheat_energy
-
-        secondary.pressure = min(
-            constants.MAX_PRESSURE, max(secondary.pressure, saturation_pressure(secondary.temperature_out))
-        )
+        self._resolve_secondary_state(secondary)
         LOGGER.debug(
             "Secondary loop: transferred=%.1fMW temp_out=%.1fK quality=%.2f energy=%.1eJ",
             transferred_mw,
@@ -157,14 +187,23 @@ class ThermalSolver:
             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:
-        """Rough CHF surrogate using mass flux and pressure."""
-        # Use a coarse mass-flux and pressure scaling to emulate higher CHF with more flow/pressure.
-        mass_flux = max(1.0, primary.mass_flow_rate / 50.0)  # kg/m2-s surrogate
-        flux_factor = max(0.5, min(3.0, (mass_flux / 200.0) ** 0.5))
-        pressure_factor = 0.5 + 0.5 * min(1.5, primary.pressure / max(0.1, constants.MAX_PRESSURE))
-        base_chf = 1.0e7  # W/m2 surrogate
-        return base_chf * flux_factor * pressure_factor
+        """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)

src/reactor_sim/turbine.py

@@ -6,7 +6,7 @@ from dataclasses import dataclass
 import logging
 
 from . import constants
-from .thermal import saturation_temperature
+from .thermal import saturation_temperature, saturation_pressure
 from .state import CoolantLoopState, TurbineState
 
 LOGGER = logging.getLogger(__name__)
@@ -45,7 +45,9 @@ class Turbine:
             state.load_demand_mw = 0.0
             state.load_supplied_mw = 0.0
             state.steam_enthalpy = 0.0
-            state.condenser_temperature = max(305.0, loop.temperature_in - 20.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))
@@ -58,13 +60,24 @@ class Turbine:
         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(loop)
+        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)
@@ -84,11 +97,12 @@ def _ramp(current: float, target: float, dt: float, time_constant: float) -> flo
     return current + (target - current) * alpha
 
 
-def _backpressure_penalty(loop: CoolantLoopState) -> float:
+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, loop.pressure))
+    pressure = max(base, min(max_p, state.condenser_pressure))
     if pressure <= base:
-        return 0.0
+        return min(constants.CONDENSER_BACKPRESSURE_PENALTY, state.fouling_penalty)
     frac = (pressure - base) / max(1e-6, max_p - base)
-    return min(constants.CONDENSER_BACKPRESSURE_PENALTY, frac * constants.CONDENSER_BACKPRESSURE_PENALTY)
+    penalty = frac * constants.CONDENSER_BACKPRESSURE_PENALTY
+    return min(constants.CONDENSER_BACKPRESSURE_PENALTY + state.fouling_penalty, penalty + state.fouling_penalty)

tests/test_simulation.py

@@ -268,10 +268,36 @@ def test_auto_control_resets_shutdown_and_moves_rods():
     assert reactor.control.rod_fraction < 0.95
 
 
-def test_full_power_reaches_steam_and_turbine_output():
-    """Integration: cold start -> pumps/gens on -> ramp to ~3 GW -> steam -> turbines online."""
+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,
@@ -279,18 +305,74 @@ def test_full_power_reaches_steam_and_turbine_output():
             generator_units={1: True, 2: True},
             primary_pumps={1: True, 2: True},
             secondary_pumps={1: True, 2: True},
-            rod_manual=False,
+            rod_manual=True,
+            rod_position=0.55,
         ),
     )
-    for i in range(600):
+    checkpoints = {300, 600, 900, 1800, 2700, 3600}
+    results = {}
+    turbines_started = False
+    for i in range(3600):
         cmd = None
-        if i == 200:
-            cmd = ReactorCommand(secondary_pumps={2: False})
-        if i == 300:
-            cmd = ReactorCommand(secondary_pumps={2: True})
-        if i == 400:
+        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 state.secondary_loop.steam_quality > 0.02
-    assert state.total_electrical_output() > 50.0
+    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

@@ -2,7 +2,7 @@ import pytest
 
 from reactor_sim import constants
 from reactor_sim.state import CoolantLoopState
-from reactor_sim.thermal import ThermalSolver, saturation_temperature
+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:
@@ -36,3 +36,15 @@ def test_secondary_generates_steam_when_energy_exceeds_sensible_heat():
     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