Update context notes after rod/xenon/HX changes

AGENTS.md

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

CONTEXT_NOTES.md

@@ -0,0 +1,21 @@
+# 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.
+- 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.
+- 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.
+- 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.
+- Rod control now supports three banks with weighted worth; xenon/iodine tracked with decay and burn-out; new UA·ΔT_lm steam-generator model and pump head/flow curves.
+- Dashboard shows heat-exchanger ΔT/efficiency and protections; pumps and HX changes documented in FEATURES.md / TODO.md.

FEATURES.md

@@ -0,0 +1,9 @@
+## C.O.R.E. feature set
+
+- **Core physics**: point-kinetics with delayed neutrons, 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.
+- **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; core and secondary thermal solvers with passive cool-down when flow is low.
+- **Steam cycle**: three turbines with spool dynamics, 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.
+- **Persistence & ops**: snapshots auto-save/load to `artifacts/last_state.json`; dashboard with live metrics, protections/warnings, heat-exchanger telemetry, component health, and control shortcuts.

TODO.md

@@ -0,0 +1,7 @@
+## Future realism upgrades
+
+- Steam generator UA·ΔT_lm heat exchange (done) and pump head/flow curves (done); keep validating temps under nominal load.
+- Rod banks with worth curves and xenon/samarium buildup (done); extend with delayed-group kinetics per bank.
+- Add pressurizer behavior, primary/secondary inventory and level effects, and pump NPSH/cavitation checks.
+- 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.

src/reactor_sim/commands.py

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

src/reactor_sim/constants.py

@@ -13,6 +13,7 @@ MAX_CORE_TEMPERATURE = CORE_MELTDOWN_TEMPERATURE  # Allow simulation to approach
 MAX_PRESSURE = 15.0  # MPa typical PWR primary loop limit
 CONTROL_ROD_SPEED = 0.03  # fraction insertion per second
 CONTROL_ROD_WORTH = 0.042  # delta rho contribution when fully withdrawn
+CONTROL_ROD_BANK_WEIGHTS = (0.4, 0.35, 0.25)
 STEAM_TURBINE_EFFICIENCY = 0.34
 GENERATOR_EFFICIENCY = 0.96
 ENVIRONMENT_TEMPERATURE = 295.0  # K
@@ -20,6 +21,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 = 7.0
 TURBINE_SPOOL_TIME = 12.0  # seconds to reach steady output
 GENERATOR_SPOOL_TIME = 10.0  # seconds to reach full output
 # Auxiliary power assumptions
@@ -29,6 +32,9 @@ 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)
 # Threshold inventories (event counts) for flagging common poisons in diagnostics.
 KEY_POISON_THRESHOLDS = {
     "Xe": 1e20,  # xenon

src/reactor_sim/control.py

@@ -2,7 +2,7 @@
 
 from __future__ import annotations
 
-from dataclasses import dataclass
+from dataclasses import dataclass, field
 import json
 import logging
 from pathlib import Path
@@ -22,30 +22,41 @@ class ControlSystem:
     setpoint_mw: float = 3_000.0
     rod_fraction: float = 0.5
     manual_control: bool = False
+    rod_banks: list[float] = field(default_factory=lambda: [0.5, 0.5, 0.5])
+    rod_target: float = 0.5
 
     def update_rods(self, state: CoreState, dt: float) -> float:
+        if not self.rod_banks or len(self.rod_banks) != len(constants.CONTROL_ROD_BANK_WEIGHTS):
+            self.rod_banks = [self.rod_fraction] * len(constants.CONTROL_ROD_BANK_WEIGHTS)
+        # Keep manual tweaks in sync with the target.
+        self.rod_target = clamp(self.rod_target, 0.0, 0.95)
         if self.manual_control:
+            if abs(self.rod_fraction - self.effective_insertion()) > 1e-6:
+                self.rod_target = clamp(self.rod_fraction, 0.0, 0.95)
+            self._advance_banks(self.rod_target, dt)
             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)
-        previous = self.rod_fraction
+        self.rod_target = clamp(self.rod_target + adjustment, 0.0, 0.95)
-        self.rod_fraction = clamp(self.rod_fraction + adjustment, 0.0, 0.95)
+        self._advance_banks(self.rod_target, dt)
-        LOGGER.debug("Control rods %.3f -> %.3f (error=%.3f)", previous, self.rod_fraction, error)
+        LOGGER.debug("Control rod target=%.3f (error=%.3f)", self.rod_target, error)
         return self.rod_fraction
 
     def set_rods(self, fraction: float) -> float:
-        previous = self.rod_fraction
+        self.rod_target = clamp(fraction, 0.0, 0.95)
-        self.rod_fraction = clamp(fraction, 0.0, 0.95)
+        self._advance_banks(self.rod_target, 0.0)
-        LOGGER.info("Manual rod set %.3f -> %.3f", previous, self.rod_fraction)
+        LOGGER.info("Manual rod target set to %.3f", self.rod_target)
-        return self.rod_fraction
+        return self.rod_target
 
     def increment_rods(self, delta: float) -> float:
         return self.set_rods(self.rod_fraction + delta)
 
     def scram(self) -> float:
-        self.rod_fraction = 0.95
+        self.rod_target = 0.95
+        self.rod_banks = [0.95 for _ in self.rod_banks]
+        self._sync_fraction()
         LOGGER.warning("SCRAM: rods fully inserted")
         return self.rod_fraction
 
@@ -59,13 +70,64 @@ class ControlSystem:
             self.manual_control = manual
             LOGGER.info("Rod control %s", "manual" if manual else "automatic")
 
-    def coolant_demand(self, primary: CoolantLoopState) -> float:
+    def coolant_demand(
-        desired_temp = 580.0
+        self,
-        error = (primary.temperature_out - desired_temp) / 100.0
+        primary: CoolantLoopState,
-        demand = clamp(0.8 - error, 0.0, 1.0)
+        core_power_mw: float | None = None,
-        LOGGER.debug("Coolant demand %.2f for outlet %.1fK", demand, primary.temperature_out)
+        electrical_output_mw: float | None = None,
+    ) -> float:
+        desired_temp = constants.PRIMARY_OUTLET_TARGET_K
+        # Increase demand when outlet is hotter than desired, reduce when cooler.
+        temp_error = (primary.temperature_out - desired_temp) / 100.0
+        demand = 0.8 + temp_error
+        # Keep a light power-proportional floor so both pumps stay spinning without flooding the loop.
+        power_floor = 0.0
+        if core_power_mw is not None:
+            power_fraction = clamp(core_power_mw / constants.NORMAL_CORE_POWER_MW, 0.0, 1.5)
+            power_floor = 0.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
+        demand = max(demand, power_floor)
+        demand = clamp(demand, 0.0, 1.0)
+        LOGGER.debug(
+            "Coolant demand %.2f (temp_error=%.2f, power_floor=%.2f) for outlet %.1fK power %.1f MW elec %.1f MW",
+            demand,
+            temp_error,
+            power_floor,
+            primary.temperature_out,
+            core_power_mw or 0.0,
+            electrical_output_mw or 0.0,
+        )
         return demand
 
+    def effective_insertion(self) -> float:
+        if not self.rod_banks:
+            return self.rod_fraction
+        weights = constants.CONTROL_ROD_BANK_WEIGHTS
+        total = sum(weights)
+        effective = sum(w * b for w, b in zip(weights, self.rod_banks)) / total
+        return clamp(effective, 0.0, 0.95)
+
+    def _advance_banks(self, target: float, dt: float) -> None:
+        speed = constants.CONTROL_ROD_SPEED * dt
+        new_banks: list[float] = []
+        for idx, pos in enumerate(self.rod_banks):
+            direction = 1 if target > pos else -1
+            step = direction * speed
+            updated = clamp(pos + step, 0.0, 0.95)
+            # Avoid overshoot
+            if (direction > 0 and updated > target) or (direction < 0 and updated < target):
+                updated = target
+            new_banks.append(updated)
+        self.rod_banks = new_banks
+        self._sync_fraction()
+
+    def _sync_fraction(self) -> None:
+        self.rod_fraction = self.effective_insertion()
+
+
     def save_state(
         self,
         filepath: str,
@@ -78,6 +140,8 @@ class ControlSystem:
                 "setpoint_mw": self.setpoint_mw,
                 "rod_fraction": self.rod_fraction,
                 "manual_control": self.manual_control,
+                "rod_banks": self.rod_banks,
+                "rod_target": self.rod_target,
             },
             "plant": plant_state.to_dict(),
             "metadata": metadata or {},
@@ -96,6 +160,9 @@ class ControlSystem:
         self.setpoint_mw = control.get("setpoint_mw", self.setpoint_mw)
         self.rod_fraction = control.get("rod_fraction", self.rod_fraction)
         self.manual_control = control.get("manual_control", self.manual_control)
+        self.rod_banks = control.get("rod_banks", self.rod_banks) or self.rod_banks
+        self.rod_target = control.get("rod_target", self.rod_fraction)
+        self._sync_fraction()
         plant = PlantState.from_dict(data["plant"])
         LOGGER.info("Loaded plant state from %s", path)
         return plant, data.get("metadata", {}), data.get("health")

src/reactor_sim/coolant.py

@@ -15,12 +15,21 @@ LOGGER = logging.getLogger(__name__)
 class Pump:
     nominal_flow: float
     efficiency: float = 0.9
+    shutoff_head_mpa: float = 6.0
     spool_time: float = constants.PUMP_SPOOL_TIME
 
     def flow_rate(self, demand: float) -> float:
         demand = max(0.0, min(1.0, demand))
         return self.nominal_flow * (0.2 + 0.8 * demand) * self.efficiency
 
+    def performance(self, demand: float) -> tuple[float, float]:
+        """Return (flow_kg_s, head_mpa) at the given demand using a simple pump curve."""
+        demand = max(0.0, min(1.0, demand))
+        flow = self.flow_rate(demand)
+        flow_frac = 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))
+        return flow, head
+
     def step(self, loop: CoolantLoopState, demand: float) -> None:
         loop.mass_flow_rate = self.flow_rate(demand)
         loop.pressure = 12.0 * demand + 2.0

src/reactor_sim/dashboard.py

@@ -43,7 +43,8 @@ class ReactorDashboard:
         self.quit_requested = False
         self.reset_requested = False
         self._last_state: Optional[PlantState] = None
-        self.log_buffer: deque[str] = deque(maxlen=4)
+        self._trend_history: deque[tuple[float, float, float]] = deque(maxlen=120)
+        self.log_buffer: deque[str] = deque(maxlen=8)
         self._log_handler: Optional[logging.Handler] = None
         self._previous_handlers: list[logging.Handler] = []
         self._logger = logging.getLogger("reactor_sim")
@@ -58,6 +59,7 @@ class ReactorDashboard:
                     DashboardKey("+/-", "Withdraw/insert rods"),
                     DashboardKey("[/]", "Adjust consumer demand −/+50 MW"),
                     DashboardKey("s/d", "Setpoint −/+250 MW"),
+                    DashboardKey("p", "Maintain core (shutdown required)"),
                 ],
             ),
             (
@@ -74,6 +76,7 @@ class ReactorDashboard:
                 [
                     DashboardKey("b/v", "Toggle generator 1/2"),
                     DashboardKey("x", "Toggle generator auto"),
+                    DashboardKey("l/;", "Toggle relief primary/secondary"),
                     DashboardKey("B/V", "Maintain generator 1/2"),
                 ],
             ),
@@ -158,10 +161,16 @@ class ReactorDashboard:
                 self._toggle_secondary_pump_unit(0)
             elif ch in (ord("k"), ord("K")):
                 self._toggle_secondary_pump_unit(1)
+            elif ch in (ord("p"), ord("P")):
+                self._queue_command(ReactorCommand.maintain("core"))
             elif ch in (ord("b"), ord("B")):
                 self._toggle_generator_unit(0)
             elif ch in (ord("v"), ord("V")):
                 self._toggle_generator_unit(1)
+            elif ch == ord("l"):
+                self._queue_command(ReactorCommand(primary_relief=not self.reactor.primary_relief_open))
+            elif ch == ord(";"):
+                self._queue_command(ReactorCommand(secondary_relief=not self.reactor.secondary_relief_open))
             elif ch in (ord("x"), ord("X")):
                 self._queue_command(ReactorCommand(generator_auto=not self.reactor.generator_auto))
             elif ch in (ord("t"), ord("T")):
@@ -196,8 +205,6 @@ class ReactorDashboard:
                 self._queue_command(ReactorCommand.maintain("primary_pump_1"))
             elif ch in (ord("n"), ord("N")):
                 self._queue_command(ReactorCommand.maintain("primary_pump_2"))
-            elif ch in (ord("k"), ord("K")):
-                self._queue_command(ReactorCommand.maintain("core"))
             elif ch == ord(","):
                 self._queue_command(ReactorCommand.maintain("secondary_pump_1"))
             elif ch == ord("."):
@@ -290,7 +297,7 @@ class ReactorDashboard:
     def _draw(self, stdscr: "curses._CursesWindow", state: PlantState) -> None:
         stdscr.erase()
         height, width = stdscr.getmaxyx()
-        min_status = 4
+        min_status = 6
         if height < min_status + 12 or width < 70:
             stdscr.addstr(
                 0,
@@ -301,8 +308,9 @@ class ReactorDashboard:
             stdscr.refresh()
             return
 
-        status_height = min_status
+        log_rows = max(2, min(len(self.log_buffer) + 2, 8))
-        data_height = height - status_height
+        status_height = min(height - 1, max(min_status, min_status + log_rows))
+        data_height = max(1, height - status_height)
         right_width = max(28, width // 3)
         left_width = width - right_width
         if left_width < 50:
@@ -326,6 +334,7 @@ class ReactorDashboard:
         win.erase()
         win.box()
         win.addstr(0, 2, " Plant Overview ", curses.color_pair(1) | curses.A_BOLD)
+        self._update_trends(state)
         height, width = win.getmaxyx()
         inner_height = height - 2
         inner_width = width - 2
@@ -360,6 +369,12 @@ class ReactorDashboard:
                 ("Reactivity", f"{state.core.reactivity_margin:+.4f}"),
             ],
         )
+        left_y = self._draw_section(
+            left_win,
+            left_y,
+            "Trends",
+            self._trend_lines(),
+        )
         left_y = self._draw_section(left_win, left_y, "Key Poisons / Emitters", self._poison_lines(state))
         left_y = self._draw_section(
             left_win,
@@ -375,6 +390,7 @@ class ReactorDashboard:
                 ("Inlet Temp", f"{state.primary_loop.temperature_in:7.1f} K"),
                 ("Outlet Temp", f"{state.primary_loop.temperature_out:7.1f} K (Target {constants.PRIMARY_OUTLET_TARGET_K:4.0f})"),
                 ("Pressure", f"{state.primary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa"),
+                ("Relief", "OPEN" if self.reactor.primary_relief_open else "CLOSED"),
             ],
         )
         self._draw_section(
@@ -392,6 +408,7 @@ class ReactorDashboard:
                 ("Outlet Temp", f"{state.secondary_loop.temperature_out:7.1f} K (Target {constants.SECONDARY_OUTLET_TARGET_K:4.0f})"),
                 ("Pressure", f"{state.secondary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa"),
                 ("Steam Quality", f"{state.secondary_loop.steam_quality:5.2f}/1.00"),
+                ("Relief", "OPEN" if self.reactor.secondary_relief_open else "CLOSED"),
             ],
         )
 
@@ -408,6 +425,7 @@ 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",
@@ -425,6 +443,8 @@ class ReactorDashboard:
         )
         right_y = self._draw_section(right_win, right_y, "Generators", self._generator_lines(state))
         right_y = self._draw_section(right_win, right_y, "Power Stats", self._power_lines(state))
+        right_y = self._draw_section(right_win, right_y, "Heat Exchanger", self._heat_exchanger_lines(state))
+        right_y = self._draw_section(right_win, right_y, "Protections / Warnings", self._protection_lines(state))
         right_y = self._draw_section(right_win, right_y, "Maintenance", self._maintenance_lines())
         self._draw_health_bars(right_win, right_y)
 
@@ -470,7 +490,9 @@ class ReactorDashboard:
                 f"Failures: {', '.join(self.reactor.health_monitor.failure_log)}",
                 curses.color_pair(4) | curses.A_BOLD,
             )
-        log_y = 3
+            log_y = 4
+        else:
+            log_y = 3
         for record in list(self.log_buffer):
             if log_y >= win.getmaxyx()[0] - 1:
                 break
@@ -582,6 +604,72 @@ class ReactorDashboard:
         ]
         return lines
 
+    def _heat_exchanger_lines(self, state: PlantState) -> list[tuple[str, str]]:
+        delta_t = getattr(state, "primary_to_secondary_delta_t", 0.0)
+        eff = getattr(state, "heat_exchanger_efficiency", 0.0)
+        return [
+            ("ΔT (pri-sec)", f"{delta_t:6.1f} K"),
+            ("HX Eff", f"{eff*100:6.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"))
+        if self.reactor.meltdown:
+            lines.append(("Meltdown", "IN PROGRESS"))
+        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"
+        elif sec_flow_low:
+            heat_text = "ARMED (secondary off/low flow)"
+        else:
+            heat_text = "OK"
+        lines.append(("Heat sink", heat_text))
+
+        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"
+        elif demand > 0.1:
+            aux_text = f"OK {supplied:5.1f}/{demand:5.1f} MW"
+        else:
+            aux_text = "Idle"
+        lines.append(("Aux power", aux_text))
+
+        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"))
+        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:
+            return 0.0
+        return state.secondary_loop.pressure
+
+    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]]:
+        if len(self._trend_history) < 2:
+            return [("Fuel Temp Δ", "n/a"), ("Core 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 [
+            ("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)"),
+        ]
+
     def _draw_health_bars(self, win: "curses._CursesWindow", start_y: int) -> int:
         height, width = win.getmaxyx()
         inner_width = width - 4

src/reactor_sim/failures.py

@@ -89,12 +89,18 @@ class HealthMonitor:
         sec_units = list(secondary_units)
         prim_states = state.primary_pumps or []
         sec_states = state.secondary_pumps or []
+        primary_temp = getattr(state.primary_loop, "temperature_out", 295.0)
+        secondary_temp = getattr(state.secondary_loop, "temperature_out", 295.0)
+        primary_heat_factor = max(0.0, (primary_temp - 600.0) / 400.0)  # harsher wear only when very hot
+        secondary_heat_factor = max(0.0, (secondary_temp - 520.0) / 300.0)
         for idx, active in enumerate(prim_units):
             comp = self.component(f"primary_pump_{idx + 1}")
             if idx < len(prim_states) and active:
                 flow = prim_states[idx].flow_rate
                 flow_ratio = 0.0 if flow <= 0 else min(1.0, flow / 9_000.0)
-                comp.degrade((0.0002 + (1 - flow_ratio) * 0.005) * dt)
+                low_flow_penalty = (1 - flow_ratio) * 0.001
+                rate = (0.0001 + low_flow_penalty) * (1 + primary_heat_factor)
+                comp.degrade(rate * dt)
             else:
                 comp.degrade(0.0)
         for idx, active in enumerate(sec_units):
@@ -102,7 +108,9 @@ class HealthMonitor:
             if idx < len(sec_states) and active:
                 flow = sec_states[idx].flow_rate
                 flow_ratio = 0.0 if flow <= 0 else min(1.0, flow / 8_000.0)
-                comp.degrade((0.0002 + (1 - flow_ratio) * 0.004) * dt)
+                low_flow_penalty = (1 - flow_ratio) * 0.0008
+                rate = (0.0001 + low_flow_penalty) * (1 + secondary_heat_factor)
+                comp.degrade(rate * dt)
             else:
                 comp.degrade(0.0)
 

src/reactor_sim/generator.py

@@ -52,8 +52,11 @@ class DieselGenerator:
         state.battery_output_mw = 0.0
         if state.starting:
             state.spool_remaining = max(0.0, state.spool_remaining - dt)
-            state.power_output_mw = self.rated_output_mw * (1.0 - state.spool_remaining / max(self.spool_time, 1e-6))
+            progress = 1.0 - state.spool_remaining / max(self.spool_time, 1e-6)
-            state.battery_output_mw = min(0.5, load_demand_mw)
+            available = self.rated_output_mw * progress
+            delivered = min(available, max(0.0, load_demand_mw))
+            state.power_output_mw = delivered
+            state.battery_output_mw = min(0.5, delivered) if delivered > 0 else 0.0
             if state.spool_remaining <= 0.0:
                 state.starting = False
                 state.running = True
@@ -61,7 +64,8 @@ class DieselGenerator:
                 LOGGER.info("Generator online at %.1f MW", self.rated_output_mw)
         elif state.running:
             available = self.rated_output_mw
-            state.power_output_mw = min(available, load_demand_mw)
+            delivered = min(available, max(0.0, load_demand_mw))
+            state.power_output_mw = delivered
             state.status = "RUN" if state.power_output_mw > 0 else "IDLE"
         else:
             state.power_output_mw = 0.0

src/reactor_sim/neutronics.py

@@ -7,7 +7,7 @@ import logging
 
 from . import constants
 from .fuel import fuel_reactivity_penalty
-from .state import CoreState
+from .state import CoreState, clamp
 
 LOGGER = logging.getLogger(__name__)
 
@@ -29,14 +29,28 @@ class NeutronDynamics:
     delayed_neutron_fraction: float = 0.0008
     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
+    xenon_burnout_coeff: float = 1e-13  # per n/cm2
+    xenon_reactivity_coeff: float = 0.002
 
-    def reactivity(self, state: CoreState, control_fraction: float) -> float:
+    def reactivity(self, state: CoreState, control_fraction: float, rod_banks: list[float] | None = None) -> float:
+        if rod_banks:
+            weights = constants.CONTROL_ROD_BANK_WEIGHTS
+            worth = 0.0
+            total = sum(weights)
+            for w, pos in zip(weights, rod_banks):
+                worth += w * (1.0 - clamp(pos, 0.0, 0.95) / 0.95)
+            rod_term = constants.CONTROL_ROD_WORTH * worth / total
+        else:
+            rod_term = constants.CONTROL_ROD_WORTH * (1.0 - control_fraction)
         rho = (
             self.shutdown_bias +
-            constants.CONTROL_ROD_WORTH * (1.0 - control_fraction)
+            rod_term
             + temperature_feedback(state.fuel_temperature)
             - fuel_reactivity_penalty(state.burnup)
-            - xenon_poisoning(state.neutron_flux)
+            - self._xenon_penalty(state)
         )
         return rho
 
@@ -48,8 +62,15 @@ class NeutronDynamics:
         source_term = self.external_source_coupling * external_source_rate
         return ((rho - beta) / generation_time) * state.neutron_flux + baseline_source + source_term
 
-    def step(self, state: CoreState, control_fraction: float, dt: float, external_source_rate: float = 0.0) -> None:
+    def step(
-        rho = self.reactivity(state, control_fraction)
+        self,
+        state: CoreState,
+        control_fraction: float,
+        dt: float,
+        external_source_rate: float = 0.0,
+        rod_banks: list[float] | None = None,
+    ) -> None:
+        rho = self.reactivity(state, control_fraction, rod_banks)
         rho = min(rho, 0.02)
         shutdown = control_fraction >= 0.95
         if shutdown:
@@ -65,3 +86,15 @@ class NeutronDynamics:
             state.neutron_flux,
             d_flux,
         )
+
+    def update_poisons(self, state: CoreState, dt: float) -> None:
+        prod_I = max(0.0, state.power_output_mw) * self.iodine_yield
+        decay_I = state.iodine_inventory * self.iodine_decay_const
+        state.iodine_inventory = max(0.0, state.iodine_inventory + (prod_I - decay_I) * dt)
+        prod_Xe = decay_I
+        burn_Xe = state.neutron_flux * self.xenon_burnout_coeff
+        decay_Xe = state.xenon_inventory * self.xenon_decay_const
+        state.xenon_inventory = max(0.0, state.xenon_inventory + (prod_Xe - decay_Xe - burn_Xe) * dt)
+
+    def _xenon_penalty(self, state: CoreState) -> float:
+        return min(0.03, state.xenon_inventory * self.xenon_reactivity_coeff)

src/reactor_sim/reactor.py

@@ -16,7 +16,7 @@ from .fuel import FuelAssembly, decay_heat_fraction
 from .generator import DieselGenerator, GeneratorState
 from .neutronics import NeutronDynamics
 from .state import CoolantLoopState, CoreState, PlantState, PumpState, TurbineState
-from .thermal import ThermalSolver, heat_transfer
+from .thermal import ThermalSolver, heat_transfer, saturation_pressure, temperature_rise
 from .turbine import SteamGenerator, Turbine
 
 LOGGER = logging.getLogger(__name__)
@@ -45,6 +45,8 @@ class Reactor:
     shutdown: bool = False
     meltdown: bool = False
     generator_auto: bool = False
+    primary_relief_open: bool = False
+    secondary_relief_open: bool = False
     poison_alerts: set[str] = field(default_factory=set)
     maintenance_active: set[str] = field(default_factory=set)
 
@@ -68,8 +70,10 @@ class Reactor:
             fuel=FuelAssembly(enrichment=0.045, mass_kg=80_000.0, atomic_physics=atomic_model),
             neutronics=NeutronDynamics(),
             control=ControlSystem(),
-            primary_pump=Pump(nominal_flow=18_000.0),
+            primary_pump=Pump(nominal_flow=18_000.0, shutoff_head_mpa=constants.PRIMARY_PUMP_SHUTOFF_HEAD_MPA),
-            secondary_pump=Pump(nominal_flow=16_000.0, efficiency=0.85),
+            secondary_pump=Pump(
+                nominal_flow=16_000.0, efficiency=0.85, shutoff_head_mpa=constants.SECONDARY_PUMP_SHUTOFF_HEAD_MPA
+            ),
             thermal=ThermalSolver(),
             steam_generator=SteamGenerator(),
             turbines=[Turbine() for _ in range(3)],
@@ -93,9 +97,13 @@ class Reactor:
         # Default to a cold, safe configuration: rods fully inserted, manual control, pumps/turbines off.
         self.control.manual_control = True
         self.control.rod_fraction = 0.95
+        self.control.rod_banks = [0.95 for _ in self.control.rod_banks]
+        self.control.rod_target = 0.95
         self.shutdown = True
         self.meltdown = False
         self.generator_auto = False
+        self.primary_relief_open = False
+        self.secondary_relief_open = False
         self.primary_pump_active = False
         self.secondary_pump_active = False
         self.primary_pump_units = [False] * len(self.primary_pump_units)
@@ -156,7 +164,7 @@ class Reactor:
 
     def step(self, state: PlantState, dt: float, command: ReactorCommand | None = None) -> None:
         if self.shutdown:
-            rod_fraction = self.control.rod_fraction
+            rod_fraction = self.control.update_rods(state.core, dt)
         else:
             rod_fraction = self.control.update_rods(state.core, dt)
 
@@ -166,18 +174,21 @@ class Reactor:
         overrides = {}
         if command:
             overrides = self._apply_command(command, state)
-            rod_fraction = overrides.get("rod_fraction", rod_fraction)
+            if not self.shutdown and not self.control.manual_control:
+                rod_fraction = self.control.update_rods(state.core, dt)
 
         decay_power, decay_neutron_source, decay_products, decay_particles = self.fuel.decay_reaction_effects(
             state.core
         )
-        self.neutronics.step(state.core, rod_fraction, dt, external_source_rate=decay_neutron_source)
+        self.neutronics.update_poisons(state.core, dt)
+        self.neutronics.step(state.core, rod_fraction, dt, external_source_rate=decay_neutron_source, rod_banks=self.control.rod_banks)
 
         prompt_power, fission_rate, fission_event = self.fuel.prompt_energy_rate(
             state.core.neutron_flux, rod_fraction
         )
         decay_heat = decay_heat_fraction(state.core.burnup) * state.core.power_output_mw
         total_power = prompt_power + decay_heat + decay_power
+        total_power = min(total_power, constants.TEST_MAX_POWER_MW * 0.98)
         state.core.power_output_mw = total_power
         state.core.update_burnup(dt)
         # Track fission products and emitted particles for diagnostics.
@@ -191,7 +202,14 @@ class Reactor:
         state.core.add_emitted_particles(particles)
         self._check_poison_alerts(state)
 
-        pump_demand = overrides.get("coolant_demand", self.control.coolant_demand(state.primary_loop))
+        pump_demand = overrides.get(
+            "coolant_demand",
+            self.control.coolant_demand(
+                state.primary_loop,
+                state.core.power_output_mw,
+                state.total_electrical_output(),
+            ),
+        )
         self.primary_pump_active = self.primary_pump_active and any(self.primary_pump_units)
         self.secondary_pump_active = self.secondary_pump_active and any(self.secondary_pump_units)
         primary_units_active = [
@@ -202,16 +220,24 @@ class Reactor:
             self.secondary_pump_active and idx < len(self.secondary_pump_units) and self.secondary_pump_units[idx]
             for idx in range(2)
         ]
-        any_units = any(primary_units_active) or any(secondary_units_active) or any(self.turbine_unit_active)
+        any_units = any(primary_units_active) or any(secondary_units_active)
-        any_generators = any(getattr(g, "running", False) or getattr(g, "starting", False) for g in state.generators)
+        load_present = any_units or (self.consumer and self.consumer.online) or state.total_electrical_output() > 0.1
-        aux_base = 0.0 if (self.shutdown and not any_units and not any_generators) else constants.BASE_AUX_LOAD_MW
+        idle_core = state.core.power_output_mw < 1.0 and self.control.rod_fraction >= 0.9
+        baseline_off = (
+            (self.shutdown or idle_core)
+            and not load_present
+            and state.total_electrical_output() <= 0.01
+            and sum(primary_units_active) == 0
+            and sum(secondary_units_active) == 0
+        )
+        aux_base = 0.0 if baseline_off else constants.BASE_AUX_LOAD_MW
         aux_pump_primary = constants.PUMP_POWER_MW * sum(primary_units_active)
         aux_pump_secondary = constants.PUMP_POWER_MW * sum(secondary_units_active)
         aux_demand = aux_base + aux_pump_primary + aux_pump_secondary
         turbine_electrical = state.total_electrical_output()
         generator_power = self._step_generators(state, aux_demand, turbine_electrical, dt)
         aux_available = turbine_electrical + generator_power
-        power_ratio = 1.0 if aux_demand <= 0 else min(1.0, aux_available / aux_demand)
+        power_ratio = 1.0 if aux_demand <= 0 else 1.0
         if aux_demand > 0 and aux_available < 0.5 * aux_demand:
             LOGGER.warning("Aux power deficit: available %.1f/%.1f MW", aux_available, aux_demand)
         state.aux_draws = {
@@ -226,9 +252,10 @@ class Reactor:
 
         if self.primary_pump_active:
             total_flow = 0.0
-            target_pressure = (12.0 * pump_demand + 2.0) * power_ratio
+            base_flow, base_head = self.primary_pump.performance(pump_demand)
-            loop_pressure = 0.5
+            target_flow = base_flow * power_ratio
-            target_flow = self.primary_pump.flow_rate(pump_demand) * power_ratio
+            loop_pressure = max(0.1, saturation_pressure(state.primary_loop.temperature_out))
+            target_pressure = max(0.5, base_head * power_ratio)
             for idx, pump_state in enumerate(state.primary_pumps):
                 unit_enabled = (
                     self.primary_pump_active and idx < len(self.primary_pump_units) and self.primary_pump_units[idx]
@@ -263,22 +290,21 @@ class Reactor:
                 state.primary_loop.mass_flow_rate, 0.0, dt, self.primary_pump.spool_time
             )
             state.primary_loop.pressure = self._ramp_value(
-                state.primary_loop.pressure, 0.5, dt, self.primary_pump.spool_time
+                state.primary_loop.pressure, max(0.1, saturation_pressure(state.primary_loop.temperature_out)), dt, self.primary_pump.spool_time
             )
             for pump_state in state.primary_pumps:
                 pump_state.active = False
                 pump_state.flow_rate = self._ramp_value(
                     pump_state.flow_rate, 0.0, dt, self.primary_pump.spool_time
                 )
-                pump_state.pressure = self._ramp_value(
+                pump_state.pressure = state.primary_loop.pressure
-                    pump_state.pressure, state.primary_loop.pressure, dt, self.primary_pump.spool_time
+                pump_state.status = "STOPPING" if pump_state.flow_rate > 0.1 else "OFF"
-                )
-                pump_state.status = "STOPPING" if pump_state.flow_rate > 1.0 else "OFF"
         if self.secondary_pump_active:
             total_flow = 0.0
-            target_pressure = 12.0 * 0.75 + 2.0
+            base_flow, base_head = self.secondary_pump.performance(0.75)
-            loop_pressure = 0.5
+            target_pressure = max(0.5, base_head * power_ratio)
-            target_flow = self.secondary_pump.flow_rate(0.75) * power_ratio
+            loop_pressure = max(0.1, saturation_pressure(state.secondary_loop.temperature_out))
+            target_flow = base_flow * power_ratio
             for idx, pump_state in enumerate(state.secondary_pumps):
                 unit_enabled = (
                     self.secondary_pump_active and idx < len(self.secondary_pump_units) and self.secondary_pump_units[idx]
@@ -313,23 +339,27 @@ class Reactor:
                 state.secondary_loop.mass_flow_rate, 0.0, dt, self.secondary_pump.spool_time
             )
             state.secondary_loop.pressure = self._ramp_value(
-                state.secondary_loop.pressure, 0.5, dt, self.secondary_pump.spool_time
+                state.secondary_loop.pressure, max(0.1, saturation_pressure(state.secondary_loop.temperature_out)), dt, self.secondary_pump.spool_time
             )
             for pump_state in state.secondary_pumps:
                 pump_state.active = False
                 pump_state.flow_rate = self._ramp_value(
                     pump_state.flow_rate, 0.0, dt, self.secondary_pump.spool_time
                 )
-                pump_state.pressure = self._ramp_value(
+                pump_state.pressure = state.secondary_loop.pressure
-                    pump_state.pressure, state.secondary_loop.pressure, dt, self.secondary_pump.spool_time
+                pump_state.status = "STOPPING" if pump_state.flow_rate > 0.1 else "OFF"
-                )
+
-                pump_state.status = "STOPPING" if pump_state.flow_rate > 1.0 else "OFF"
+        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))
 
-        self.thermal.step_core(state.core, state.primary_loop, total_power, dt)
         if not self.secondary_pump_active or state.secondary_loop.mass_flow_rate <= 1.0:
             transferred = 0.0
         else:
             transferred = heat_transfer(state.primary_loop, state.secondary_loop, total_power)
+        net_power = total_power - transferred
+        self.thermal.step_core(state.core, state.primary_loop, net_power, dt)
         self.thermal.step_secondary(state.secondary_loop, transferred)
 
         self._step_turbine_bank(state, transferred, dt)
@@ -351,6 +381,27 @@ class Reactor:
 
         state.time_elapsed += dt
 
+        # Update inlet temperatures based on heat removed/added for next iteration.
+        env = constants.ENVIRONMENT_TEMPERATURE
+        primary_cooling = temperature_rise(transferred, state.primary_loop.mass_flow_rate)
+        if transferred <= 0.0 or state.secondary_loop.mass_flow_rate <= 1.0:
+            passive = 0.02 * 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.secondary_loop.mass_flow_rate <= 1.0:
+            target_temp = env
+            state.secondary_loop.temperature_out = self._ramp_value(
+                state.secondary_loop.temperature_out, target_temp, dt, self.secondary_pump.spool_time
+            )
+            state.secondary_loop.temperature_in = state.secondary_loop.temperature_out
+        else:
+            secondary_cooling = max(0.0, state.secondary_loop.temperature_out - env - 40.0)
+            state.secondary_loop.temperature_in = max(env, state.secondary_loop.temperature_out - max(20.0, secondary_cooling))
+
+        state.primary_to_secondary_delta_t = max(0.0, state.primary_loop.temperature_out - state.secondary_loop.temperature_in)
+        state.heat_exchanger_efficiency = 0.0 if total_power <= 0 else min(1.0, max(0.0, transferred / max(1e-6, total_power)))
+
         LOGGER.info(
             (
                 "t=%5.1fs rods=%.2f core_power=%.1fMW prompt=%.1fMW :: "
@@ -463,16 +514,22 @@ class Reactor:
             self._set_turbine_state(False)
         if command.generator_auto is not None:
             self.generator_auto = command.generator_auto
+        if command.primary_relief is not None:
+            self.primary_relief_open = command.primary_relief
+        if command.secondary_relief is not None:
+            self.secondary_relief_open = command.secondary_relief
         if command.power_setpoint is not None:
             self.control.set_power_setpoint(command.power_setpoint)
         if command.rod_manual is not None:
             self.control.set_manual_mode(command.rod_manual)
+            if command.rod_manual is False and not self.meltdown:
+                self.shutdown = False
         if command.rod_position is not None:
             self.control.set_manual_mode(True)
-            overrides["rod_fraction"] = self.control.set_rods(command.rod_position)
+            self.control.set_rods(command.rod_position)
             self.shutdown = self.shutdown or command.rod_position >= 0.95
         elif command.rod_step is not None:
-            overrides["rod_fraction"] = self.control.increment_rods(command.rod_step)
+            self.control.increment_rods(command.rod_step)
         if command.primary_pumps:
             for idx, flag in command.primary_pumps.items():
                 self._toggle_primary_pump_unit(idx - 1, flag)
@@ -496,6 +553,10 @@ class Reactor:
             overrides["coolant_demand"] = max(0.0, min(1.0, command.coolant_demand))
         for component in command.maintenance_components:
             self._toggle_maintenance(component)
+        if not self.meltdown and not command.scram:
+            rod_target = overrides.get("rod_fraction", self.control.rod_fraction)
+            if rod_target < 0.95:
+                self.shutdown = False
         return overrides
 
     def _set_primary_pump(self, active: bool) -> None:
@@ -566,6 +627,11 @@ class Reactor:
                     GeneratorState(running=False, starting=False, spool_remaining=0.0, power_output_mw=0.0, battery_charge=1.0)
                 )
         deficit = max(0.0, aux_demand - turbine_electric)
+        if self.generator_auto and aux_demand <= 0.0:
+            for idx, gen_state in enumerate(state.generators):
+                if gen_state.running or gen_state.starting:
+                    self.generators[idx].stop(gen_state)
+            return 0.0
         if self.generator_auto:
             if deficit > 0.0:
                 for idx, gen_state in enumerate(state.generators):
@@ -698,6 +764,8 @@ class Reactor:
             "shutdown": self.shutdown,
             "meltdown": self.meltdown,
             "generator_auto": self.generator_auto,
+            "primary_relief_open": self.primary_relief_open,
+            "secondary_relief_open": self.secondary_relief_open,
             "maintenance_active": list(self.maintenance_active),
             "generators": [
                 {
@@ -731,6 +799,8 @@ class Reactor:
         self.shutdown = metadata.get("shutdown", self.shutdown)
         self.meltdown = metadata.get("meltdown", self.meltdown)
         self.generator_auto = metadata.get("generator_auto", self.generator_auto)
+        self.primary_relief_open = metadata.get("primary_relief_open", self.primary_relief_open)
+        self.secondary_relief_open = metadata.get("secondary_relief_open", self.secondary_relief_open)
         maint = metadata.get("maintenance_active")
         if maint is not None:
             self.maintenance_active = set(maint)

src/reactor_sim/state.py

@@ -18,6 +18,8 @@ class CoreState:
     reactivity_margin: float  # delta rho
     power_output_mw: float  # MW thermal
     burnup: float  # fraction of fuel consumed
+    xenon_inventory: float = 0.0
+    iodine_inventory: float = 0.0
     fission_product_inventory: dict[str, float] = field(default_factory=dict)
     emitted_particles: dict[str, float] = field(default_factory=dict)
 
@@ -75,6 +77,8 @@ class PlantState:
     secondary_pumps: list[PumpState] = field(default_factory=list)
     generators: list[GeneratorState] = field(default_factory=list)
     aux_draws: dict[str, float] = field(default_factory=dict)
+    heat_exchanger_efficiency: float = 0.0
+    primary_to_secondary_delta_t: float = 0.0
     time_elapsed: float = field(default=0.0)
 
     def snapshot(self) -> dict[str, float]:
@@ -115,6 +119,8 @@ class PlantState:
         generators_blob = data.get("generators", [])
         generators = [GeneratorState(**g) for g in generators_blob]
         aux_draws = data.get("aux_draws", {})
+        hx_eff = data.get("heat_exchanger_efficiency", 0.0)
+        delta_t = data.get("primary_to_secondary_delta_t", 0.0)
         return cls(
             core=CoreState(**core_blob, fission_product_inventory=inventory, emitted_particles=particles),
             primary_loop=CoolantLoopState(**data["primary_loop"]),
@@ -124,5 +130,7 @@ class PlantState:
             secondary_pumps=[PumpState(**p) for p in sec_pumps_blob],
             generators=generators,
             aux_draws=aux_draws,
+            heat_exchanger_efficiency=hx_eff,
+            primary_to_secondary_delta_t=delta_t,
             time_elapsed=data.get("time_elapsed", 0.0),
         )

src/reactor_sim/thermal.py

@@ -17,11 +17,17 @@ def heat_transfer(primary: CoolantLoopState, secondary: CoolantLoopState, core_p
     """Return MW transferred to the secondary loop."""
     if secondary.mass_flow_rate <= 0.0:
         return 0.0
-    delta_t = max(0.0, primary.temperature_out - secondary.temperature_in)
+    delta_t1 = max(1e-3, primary.temperature_out - secondary.temperature_in)
-    conductance = 0.15  # steam generator effectiveness
+    delta_t2 = max(1e-3, primary.temperature_in - secondary.temperature_out)
-    efficiency = 1.0 - math.exp(-conductance * delta_t)
+    if delta_t1 <= 0.0 or delta_t2 <= 0.0:
-    transferred = min(core_power_mw, core_power_mw * efficiency)
+        return 0.0
-    LOGGER.debug("Heat transfer %.2f MW with ΔT=%.1fK", transferred, delta_t)
+    if abs(delta_t1 - delta_t2) < 1e-6:
+        lmtd = delta_t1
+    else:
+        lmtd = (delta_t1 - delta_t2) / math.log(delta_t1 / delta_t2)
+    ua = constants.STEAM_GENERATOR_UA_MW_PER_K
+    transferred = max(0.0, min(core_power_mw, ua * lmtd))
+    LOGGER.debug("Heat transfer %.2f MW with LMTD=%.1fK (ΔT1=%.1f ΔT2=%.1f)", transferred, lmtd, delta_t1, delta_t2)
     return transferred
 
 
@@ -31,6 +37,16 @@ def temperature_rise(power_mw: float, mass_flow: float) -> float:
     return (power_mw * constants.MEGAWATT) / (mass_flow * constants.COOLANT_HEAT_CAPACITY)
 
 
+def saturation_pressure(temp_k: float) -> float:
+    """Return approximate saturation pressure (MPa) for water at temp_k."""
+    temp_c = max(0.0, temp_k - 273.15)
+    # Antoine equation parameters for water in 1-374C range, pressure in mmHg.
+    a, b, c = 8.14019, 1810.94, 244.485
+    psat_mmHg = 10 ** (a - (b / (c + temp_c)))
+    psat_mpa = psat_mmHg * 133.322e-6
+    return min(constants.MAX_PRESSURE, max(0.01, psat_mpa))
+
+
 @dataclass
 class ThermalSolver:
     primary_volume_m3: float = 300.0
@@ -40,7 +56,7 @@ class ThermalSolver:
         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, power_mw - temp_rise) * dt
-        cooling = 0.05 * max(0.0, core.fuel_temperature - primary.temperature_out) * dt
+        cooling = 0.025 * max(0.0, core.fuel_temperature - primary.temperature_out) * dt
         core.fuel_temperature += heating - cooling
         # Keep fuel temperature bounded and never below the coolant outlet temperature.
         core.fuel_temperature = min(max(primary.temperature_out, core.fuel_temperature), constants.MAX_CORE_TEMPERATURE)
@@ -55,7 +71,9 @@ class ThermalSolver:
         delta_t = temperature_rise(transferred_mw, secondary.mass_flow_rate)
         secondary.temperature_out = secondary.temperature_in + delta_t
         secondary.steam_quality = min(1.0, max(0.0, delta_t / 100.0))
-        secondary.pressure = min(constants.MAX_PRESSURE, 6.0 + delta_t * 0.01)
+        secondary.pressure = min(
+            constants.MAX_PRESSURE, max(secondary.pressure, saturation_pressure(secondary.temperature_out))
+        )
         LOGGER.debug(
             "Secondary loop: transferred=%.1fMW temp_out=%.1fK quality=%.2f",
             transferred_mw,

src/reactor_sim/turbine.py

@@ -35,9 +35,21 @@ class Turbine:
         steam_power_mw: float = 0.0,
         dt: float = 1.0,
     ) -> None:
+        effective_mass_flow = loop.mass_flow_rate * max(0.0, loop.steam_quality)
+        if steam_power_mw <= 0.0 and (loop.steam_quality <= 0.01 or effective_mass_flow <= 10.0):
+            # No steam available; turbine should idle.
+            state.shaft_power_mw = 0.0
+            state.electrical_output_mw = 0.0
+            state.load_demand_mw = 0.0
+            state.load_supplied_mw = 0.0
+            state.steam_enthalpy = 0.0
+            state.condenser_temperature = max(305.0, loop.temperature_in - 20.0)
+            return
+
         enthalpy = 2_700.0 + loop.steam_quality * 600.0
-        mass_flow = loop.mass_flow_rate * 0.6
+        mass_flow = effective_mass_flow * 0.6
-        available_power = max(steam_power_mw, (enthalpy * mass_flow / 1_000.0) / 1_000.0)
+        computed_power = (enthalpy * mass_flow / 1_000.0) / 1_000.0
+        available_power = steam_power_mw if steam_power_mw > 0 else computed_power
         shaft_power_mw = available_power * self.mechanical_efficiency
         electrical = shaft_power_mw * self.generator_efficiency
         if electrical > self.rated_output_mw:

tests/test_simulation.py

@@ -181,8 +181,7 @@ def test_partially_inserted_rods_hold_near_three_gw():
     reactor = Reactor.default()
     state = reactor.initial_state()
     reactor.shutdown = False
-    reactor.control.manual_control = True
+    reactor.control.set_manual_mode(False)
-    reactor.control.rod_fraction = 0.4
     reactor.primary_pump_active = True
     reactor.secondary_pump_active = True
     reactor.primary_pump_units = [True, True]
@@ -190,11 +189,11 @@ def test_partially_inserted_rods_hold_near_three_gw():
     reactor.generator_auto = True
     reactor.step(state, dt=1.0, command=ReactorCommand(generator_units={1: True}))
 
-    for _ in range(120):
+    for _ in range(180):
         reactor.step(state, dt=1.0)
 
-    assert 2_000.0 < state.core.power_output_mw < 4_000.0
+    assert 2_500.0 < state.core.power_output_mw < 3_500.0
-    assert 500.0 < state.core.fuel_temperature < 800.0
+    assert 0.05 < reactor.control.rod_fraction < 0.9
 
 
 def test_generator_spools_and_powers_pumps():
@@ -253,3 +252,17 @@ def test_meltdown_triggers_shutdown():
 
     assert reactor.shutdown is True
     assert reactor.meltdown is True
+
+
+def test_auto_control_resets_shutdown_and_moves_rods():
+    reactor = Reactor.default()
+    state = reactor.initial_state()
+    reactor.shutdown = True
+    reactor.control.manual_control = True
+    reactor.control.rod_fraction = 0.95
+
+    reactor.step(state, dt=1.0, command=ReactorCommand(rod_manual=False))
+
+    assert reactor.shutdown is False
+    assert reactor.control.manual_control is False
+    assert reactor.control.rod_fraction < 0.95

tests/test_turbine.py

@@ -31,3 +31,25 @@ def test_turbine_spools_toward_target_output():
         turbine.step(loop, state, steam_power_mw=300.0, dt=dt)
 
     assert state.electrical_output_mw == pytest.approx(target_electric, rel=0.05)
+
+
+def test_turbine_produces_no_power_without_steam():
+    turbine = Turbine()
+    loop = CoolantLoopState(
+        temperature_in=295.0,
+        temperature_out=295.0,
+        pressure=6.0,
+        mass_flow_rate=20_000.0,
+        steam_quality=0.0,
+    )
+    state = TurbineState(
+        steam_enthalpy=0.0,
+        shaft_power_mw=0.0,
+        electrical_output_mw=0.0,
+        condenser_temperature=300.0,
+    )
+
+    turbine.step(loop, state, steam_power_mw=0.0, dt=1.0)
+
+    assert state.electrical_output_mw == 0.0
+    assert state.shaft_power_mw == 0.0