Increase xenon reactivity impact

FEATURES.md

@@ -3,7 +3,8 @@
 - **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.
+- **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 when flow is low.
+- **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.
 - **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.

src/reactor_sim/constants.py

@@ -14,6 +14,7 @@ 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)
+ROD_MANUAL_STEP = 0.025
 STEAM_TURBINE_EFFICIENCY = 0.34
 GENERATOR_EFFICIENCY = 0.96
 ENVIRONMENT_TEMPERATURE = 295.0  # K
@@ -24,6 +25,14 @@ 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
+
+# Turbine/condenser parameters
+TURBINE_THROTTLE_MIN = 0.1
+TURBINE_THROTTLE_MAX = 1.0
+TURBINE_THROTTLE_EFFICIENCY_DROP = 0.15  # efficiency loss when at minimum throttle
+CONDENSER_BASE_PRESSURE_MPA = 0.01
+CONDENSER_MAX_PRESSURE_MPA = 0.3
+CONDENSER_BACKPRESSURE_PENALTY = 0.35  # fractional power loss at max back-pressure
 GENERATOR_SPOOL_TIME = 10.0  # seconds to reach full output
 # Auxiliary power assumptions
 PUMP_POWER_MW = 12.0  # MW draw per pump unit
@@ -35,6 +44,21 @@ 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)
+# 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_HEAT_RATE_MPA_PER_S = 0.08
+PRIMARY_PRESSURIZER_SPRAY_RATE_MPA_PER_S = 0.12
+PRIMARY_PRESSURIZER_LEVEL_DRAW_PER_S = 0.002
+PRIMARY_PRESSURIZER_LEVEL_FILL_PER_S = 0.001
+LOOP_INVENTORY_CORRECTION_RATE = 0.08  # fraction of nominal mass restored per second toward target
+PRIMARY_INVENTORY_TARGET = 0.95
+SECONDARY_INVENTORY_TARGET = 0.9
+SECONDARY_STEAM_LOSS_FRACTION = 0.02  # fraction of steam mass that leaves the loop each second
+NPSH_REQUIRED_MPA = 0.25
+LOW_LEVEL_FLOW_FLOOR = 0.05
 # Threshold inventories (event counts) for flagging common poisons in diagnostics.
 KEY_POISON_THRESHOLDS = {
     "Xe": 1e20,  # xenon

src/reactor_sim/control.py

@@ -45,7 +45,7 @@ class ControlSystem:
         return self.rod_fraction
 
     def set_rods(self, fraction: float) -> float:
-        self.rod_target = clamp(fraction, 0.0, 0.95)
+        self.rod_target = self._quantize_manual(fraction)
         self._advance_banks(self.rod_target, 0.0)
         LOGGER.info("Manual rod target set to %.3f", self.rod_target)
         return self.rod_target
@@ -127,6 +127,11 @@ class ControlSystem:
     def _sync_fraction(self) -> None:
         self.rod_fraction = self.effective_insertion()
 
+    def _quantize_manual(self, fraction: float) -> float:
+        step = constants.ROD_MANUAL_STEP
+        quantized = round(fraction / step) * step
+        return clamp(quantized, 0.0, 0.95)
+
 
     def save_state(
         self,

src/reactor_sim/dashboard.py

@@ -148,9 +148,6 @@ class ReactorDashboard:
                 return
             if ch == ord(" "):
                 self._queue_command(ReactorCommand.scram_all())
-            elif ch in (ord("p"), ord("P")):
-                # Deprecated master toggles ignored.
-                continue
             elif ch in (ord("o"), ord("O")):
                 continue
             elif ch in (ord("g"), ord("G")):
@@ -180,10 +177,10 @@ class ReactorDashboard:
                 self._toggle_turbine_unit(idx)
             elif ch in (ord("+"), ord("=")):
                 # Insert rods (increase fraction)
-                self._queue_command(ReactorCommand(rod_position=self._clamped_rod(0.05)))
+                self._queue_command(ReactorCommand(rod_position=self._clamped_rod(constants.ROD_MANUAL_STEP)))
             elif ch == ord("-"):
                 # Withdraw rods (decrease fraction)
-                self._queue_command(ReactorCommand(rod_position=self._clamped_rod(-0.05)))
+                self._queue_command(ReactorCommand(rod_position=self._clamped_rod(-constants.ROD_MANUAL_STEP)))
             elif ch == ord("["):
                 demand = self._current_demand() - 50.0
                 self._queue_command(ReactorCommand(consumer_demand=max(0.0, demand)))
@@ -387,9 +384,11 @@ class ReactorDashboard:
                     "Flow",
                     f"{state.primary_loop.mass_flow_rate:7.0f}/{self.reactor.primary_pump.nominal_flow * len(self.reactor.primary_pump_units):.0f} kg/s",
                 ),
+                ("Level", f"{state.primary_loop.level*100:6.1f}%"),
                 ("Inlet Temp", f"{state.primary_loop.temperature_in:7.1f} K"),
                 ("Outlet Temp", f"{state.primary_loop.temperature_out:7.1f} K (Target {constants.PRIMARY_OUTLET_TARGET_K:4.0f})"),
                 ("Pressure", f"{state.primary_loop.pressure:5.2f}/{constants.MAX_PRESSURE:4.1f} MPa"),
+                ("Pressurizer", f"{self.reactor.pressurizer_level*100:6.1f}% @ {constants.PRIMARY_PRESSURIZER_SETPOINT_MPA:4.1f} MPa"),
                 ("Relief", "OPEN" if self.reactor.primary_relief_open else "CLOSED"),
             ],
         )
@@ -404,6 +403,7 @@ 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}%"),
                 ("Inlet Temp", f"{state.secondary_loop.temperature_in:7.1f} K (Target {constants.PRIMARY_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"),
@@ -435,6 +435,7 @@ class ReactorDashboard:
                     "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}"),
@@ -548,15 +549,19 @@ class ReactorDashboard:
         inventory = state.core.fission_product_inventory or {}
         particles = state.core.emitted_particles or {}
         lines: list[tuple[str, str]] = []
-        def fmt(symbol: str, label: str) -> tuple[str, str]:
-            qty = inventory.get(symbol, 0.0)
+        def fmt(symbol: str, label: str, qty: float) -> tuple[str, str]:
             threshold = constants.KEY_POISON_THRESHOLDS.get(symbol)
             flag = " !" if threshold is not None and qty >= threshold else ""
             return (f"{label}{flag}", f"{qty:9.2e}")
 
-        lines.append(fmt("Xe", "Xe (xenon)"))
-        lines.append(fmt("Sm", "Sm (samarium)"))
-        lines.append(fmt("I", "I (iodine)"))
+        lines.append(fmt("Xe", "Xe (xenon)", getattr(state.core, "xenon_inventory", 0.0)))
+        lines.append(fmt("Sm", "Sm (samarium)", inventory.get("Sm", 0.0)))
+        lines.append(fmt("I", "I (iodine)", getattr(state.core, "iodine_inventory", 0.0)))
+        try:
+            xe_penalty = -self.reactor.neutronics.xenon_penalty(state.core)
+            lines.append(("Xe Δρ", f"{xe_penalty:+.4f}"))
+        except Exception:
+            pass
         lines.append(("Neutrons (src)", f"{particles.get('n', 0.0):9.2e}"))
         lines.append(("Gammas", f"{particles.get('gamma', 0.0):9.2e}"))
         lines.append(("Alphas", f"{particles.get('alpha', 0.0):9.2e}"))
@@ -710,7 +715,9 @@ class ReactorDashboard:
 
     def _clamped_rod(self, delta: float) -> float:
         new_fraction = self.reactor.control.rod_fraction + delta
-        return max(0.0, min(0.95, new_fraction))
+        step = constants.ROD_MANUAL_STEP
+        quantized = round(new_fraction / step) * step
+        return max(0.0, min(0.95, quantized))
 
     def _install_log_capture(self) -> None:
         if self._log_handler:

src/reactor_sim/neutronics.py

@@ -33,7 +33,7 @@ class NeutronDynamics:
     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
+    xenon_reactivity_coeff: float = 0.05
 
     def reactivity(self, state: CoreState, control_fraction: float, rod_banks: list[float] | None = None) -> float:
         if rod_banks:
@@ -50,7 +50,7 @@ class NeutronDynamics:
             rod_term
             + temperature_feedback(state.fuel_temperature)
             - fuel_reactivity_penalty(state.burnup)
-            - self._xenon_penalty(state)
+            - self.xenon_penalty(state)
         )
         return rho
 
@@ -96,5 +96,9 @@ class NeutronDynamics:
         decay_Xe = state.xenon_inventory * self.xenon_decay_const
         state.xenon_inventory = max(0.0, state.xenon_inventory + (prod_Xe - decay_Xe - burn_Xe) * dt)
 
+    def xenon_penalty(self, state: CoreState) -> float:
+        """Return delta-rho penalty from xenon inventory (positive magnitude)."""
+        return self._xenon_penalty(state)
+
     def _xenon_penalty(self, state: CoreState) -> float:
-        return min(0.03, state.xenon_inventory * self.xenon_reactivity_coeff)
+        return min(0.05, state.xenon_inventory * self.xenon_reactivity_coeff)

src/reactor_sim/reactor.py

@@ -36,6 +36,7 @@ class Reactor:
     atomic_model: AtomicPhysics
     consumer: ElectricalConsumer | None = None
     health_monitor: HealthMonitor = field(default_factory=HealthMonitor)
+    pressurizer_level: float = 0.6
     primary_pump_active: bool = True
     secondary_pump_active: bool = True
     primary_pump_units: list[bool] = field(default_factory=lambda: [True, True])
@@ -85,6 +86,9 @@ class Reactor:
 
     def initial_state(self) -> PlantState:
         ambient = constants.ENVIRONMENT_TEMPERATURE
+        primary_nominal_mass = constants.PRIMARY_LOOP_VOLUME_M3 * constants.COOLANT_DENSITY
+        secondary_nominal_mass = constants.SECONDARY_LOOP_VOLUME_M3 * constants.COOLANT_DENSITY
+        self.pressurizer_level = 0.6
         core = CoreState(
             fuel_temperature=ambient,
             neutron_flux=1e5,
@@ -119,6 +123,8 @@ class Reactor:
             pressure=0.5,
             mass_flow_rate=0.0,
             steam_quality=0.0,
+            inventory_kg=primary_nominal_mass * constants.PRIMARY_INVENTORY_TARGET,
+            level=constants.PRIMARY_INVENTORY_TARGET,
         )
         secondary = CoolantLoopState(
             temperature_in=ambient,
@@ -126,6 +132,8 @@ class Reactor:
             pressure=0.5,
             mass_flow_rate=0.0,
             steam_quality=0.0,
+            inventory_kg=secondary_nominal_mass * constants.SECONDARY_INVENTORY_TARGET,
+            level=constants.SECONDARY_INVENTORY_TARGET,
         )
         primary_pumps = [
             PumpState(active=self.primary_pump_active and self.primary_pump_units[idx], flow_rate=0.0, pressure=0.5)
@@ -202,6 +210,13 @@ class Reactor:
         state.core.add_emitted_particles(particles)
         self._check_poison_alerts(state)
 
+        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",
             self.control.coolant_demand(
@@ -237,15 +252,16 @@ 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
-        power_ratio = 1.0 if aux_demand <= 0 else 1.0
-        if aux_demand > 0 and aux_available < 0.5 * aux_demand:
+        supplied = aux_available if aux_demand <= 0 else min(aux_available, aux_demand)
+        power_ratio = 1.0 if aux_demand <= 0 else min(1.0, supplied / max(1e-6, aux_demand))
+        if aux_demand > 0 and aux_available < 0.99 * aux_demand:
             LOGGER.warning("Aux power deficit: available %.1f/%.1f MW", aux_available, aux_demand)
         state.aux_draws = {
             "base": aux_base * power_ratio,
             "primary_pumps": aux_pump_primary * power_ratio,
             "secondary_pumps": aux_pump_secondary * power_ratio,
             "total_demand": aux_demand,
-            "supplied": aux_available,
+            "supplied": supplied,
             "generator_output": generator_power,
             "turbine_output": turbine_electrical,
         }
@@ -256,12 +272,15 @@ class Reactor:
             target_flow = base_flow * power_ratio
             loop_pressure = max(0.1, saturation_pressure(state.primary_loop.temperature_out))
             target_pressure = max(0.5, base_head * power_ratio)
+            primary_flow_scale = min(
+                self._inventory_flow_scale(state.primary_loop), self._npsh_factor(state.primary_loop)
+            )
             for idx, pump_state in enumerate(state.primary_pumps):
                 unit_enabled = (
                     self.primary_pump_active and idx < len(self.primary_pump_units) and self.primary_pump_units[idx]
                 )
                 powered = power_ratio > 0.1
-                desired_flow = target_flow if unit_enabled else 0.0
+                desired_flow = target_flow * primary_flow_scale if unit_enabled else 0.0
                 desired_pressure = target_pressure if unit_enabled else 0.5
                 if not powered:
                     desired_flow = 0.0
@@ -275,6 +294,8 @@ class Reactor:
                 pump_state.active = unit_enabled and powered and pump_state.flow_rate > 1.0
                 if not powered or not unit_enabled:
                     pump_state.status = "STOPPING" if pump_state.flow_rate > 1.0 else "OFF"
+                elif primary_flow_scale < 0.99:
+                    pump_state.status = "CAV"
                 elif pump_state.flow_rate < max(1.0, desired_flow * 0.8):
                     pump_state.status = "STARTING"
                 else:
@@ -289,8 +310,9 @@ class Reactor:
             state.primary_loop.mass_flow_rate = self._ramp_value(
                 state.primary_loop.mass_flow_rate, 0.0, dt, self.primary_pump.spool_time
             )
+            target_pressure = max(0.5, saturation_pressure(state.primary_loop.temperature_out))
             state.primary_loop.pressure = self._ramp_value(
-                state.primary_loop.pressure, max(0.1, saturation_pressure(state.primary_loop.temperature_out)), dt, self.primary_pump.spool_time
+                state.primary_loop.pressure, target_pressure, dt, self.primary_pump.spool_time
             )
             for pump_state in state.primary_pumps:
                 pump_state.active = False
@@ -305,12 +327,15 @@ class Reactor:
             target_pressure = max(0.5, base_head * power_ratio)
             loop_pressure = max(0.1, saturation_pressure(state.secondary_loop.temperature_out))
             target_flow = base_flow * power_ratio
+            secondary_flow_scale = min(
+                self._inventory_flow_scale(state.secondary_loop), self._npsh_factor(state.secondary_loop)
+            )
             for idx, pump_state in enumerate(state.secondary_pumps):
                 unit_enabled = (
                     self.secondary_pump_active and idx < len(self.secondary_pump_units) and self.secondary_pump_units[idx]
                 )
                 powered = power_ratio > 0.1
-                desired_flow = target_flow if unit_enabled else 0.0
+                desired_flow = target_flow * secondary_flow_scale if unit_enabled else 0.0
                 desired_pressure = target_pressure if unit_enabled else 0.5
                 if not powered:
                     desired_flow = 0.0
@@ -324,6 +349,8 @@ class Reactor:
                 pump_state.active = unit_enabled and powered and pump_state.flow_rate > 1.0
                 if not powered or not unit_enabled:
                     pump_state.status = "STOPPING" if pump_state.flow_rate > 1.0 else "OFF"
+                elif secondary_flow_scale < 0.99:
+                    pump_state.status = "CAV"
                 elif pump_state.flow_rate < max(1.0, desired_flow * 0.8):
                     pump_state.status = "STARTING"
                 else:
@@ -349,6 +376,7 @@ class Reactor:
                 pump_state.pressure = state.secondary_loop.pressure
                 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:
@@ -358,9 +386,12 @@ class Reactor:
             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_core(state.core, state.primary_loop, total_power, dt)
         self.thermal.step_secondary(state.secondary_loop, transferred)
+        self._apply_secondary_boiloff(state, dt)
+        self._update_loop_inventory(
+            state.secondary_loop, constants.SECONDARY_LOOP_VOLUME_M3, constants.SECONDARY_INVENTORY_TARGET, dt
+        )
 
         self._step_turbine_bank(state, transferred, dt)
         self._maintenance_tick(state, dt)
@@ -430,6 +461,10 @@ class Reactor:
                 break
             turbine_state = state.turbines[idx]
             if idx in active_indices:
+                # Simple throttle map: reduce throttle when electrical demand is low, open as demand rises.
+                demand = turbine_state.load_demand_mw
+                throttle = 0.4 if demand <= 0 else min(1.0, 0.4 + demand / max(1e-6, turbine.rated_output_mw))
+                turbine.throttle = throttle
                 turbine.step(state.secondary_loop, turbine_state, steam_power_mw=power_per_unit, dt=dt)
                 if power_per_unit <= 0.0 and turbine_state.electrical_output_mw < 0.1:
                     turbine_state.status = "OFF"
@@ -488,6 +523,61 @@ class Reactor:
             turbine_state.load_demand_mw = demand_per_unit
             turbine_state.load_supplied_mw = share if demand_per_unit <= 0 else min(share, demand_per_unit)
 
+    def _nominal_inventory(self, volume_m3: float) -> float:
+        return volume_m3 * constants.COOLANT_DENSITY
+
+    def _update_loop_inventory(
+        self, loop: CoolantLoopState, volume_m3: float, target_level: float, dt: float
+    ) -> None:
+        nominal_mass = self._nominal_inventory(volume_m3)
+        if nominal_mass <= 0.0:
+            loop.level = 0.0
+            return
+        if loop.inventory_kg <= 0.0:
+            loop.inventory_kg = nominal_mass * target_level
+        current_level = loop.inventory_kg / nominal_mass
+        correction = (target_level - current_level) * constants.LOOP_INVENTORY_CORRECTION_RATE
+        loop.inventory_kg = max(0.0, loop.inventory_kg + correction * nominal_mass * dt)
+        loop.level = min(1.2, max(0.0, loop.inventory_kg / nominal_mass))
+
+    def _inventory_flow_scale(self, loop: CoolantLoopState) -> float:
+        if loop.level <= constants.LOW_LEVEL_FLOW_FLOOR:
+            return 0.0
+        if loop.level <= 0.25:
+            return max(0.0, (loop.level - constants.LOW_LEVEL_FLOW_FLOOR) / (0.25 - constants.LOW_LEVEL_FLOW_FLOOR))
+        return 1.0
+
+    def _npsh_factor(self, loop: CoolantLoopState) -> float:
+        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))
+
+    def _apply_pressurizer(self, primary: CoolantLoopState, dt: float) -> None:
+        if self.shutdown and primary.mass_flow_rate <= 100.0:
+            return
+        target = constants.PRIMARY_PRESSURIZER_SETPOINT_MPA
+        band = constants.PRIMARY_PRESSURIZER_DEADBAND_MPA
+        heat_rate = constants.PRIMARY_PRESSURIZER_HEAT_RATE_MPA_PER_S
+        spray_rate = constants.PRIMARY_PRESSURIZER_SPRAY_RATE_MPA_PER_S
+        if primary.pressure < target - band and self.pressurizer_level > 0.05:
+            primary.pressure = min(target, primary.pressure + heat_rate * dt)
+            self.pressurizer_level = max(0.0, self.pressurizer_level - constants.PRIMARY_PRESSURIZER_LEVEL_DRAW_PER_S * dt)
+        elif primary.pressure > target + band:
+            primary.pressure = max(target - band, primary.pressure - spray_rate * dt)
+            self.pressurizer_level = min(1.0, self.pressurizer_level + constants.PRIMARY_PRESSURIZER_LEVEL_FILL_PER_S * dt)
+        primary.pressure = min(constants.MAX_PRESSURE, max(saturation_pressure(primary.temperature_out), primary.pressure))
+
+    def _apply_secondary_boiloff(self, state: PlantState, dt: float) -> None:
+        loop = state.secondary_loop
+        if loop.mass_flow_rate <= 0.0 or loop.steam_quality <= 0.0:
+            return
+        steam_mass = loop.mass_flow_rate * loop.steam_quality * constants.SECONDARY_STEAM_LOSS_FRACTION * dt
+        loop.inventory_kg = max(0.0, loop.inventory_kg - steam_mass)
+        nominal = self._nominal_inventory(constants.SECONDARY_LOOP_VOLUME_M3)
+        loop.level = min(1.2, max(0.0, loop.inventory_kg / nominal)) if nominal > 0 else 0.0
+
     def _handle_failure(self, component: str) -> None:
         if component == "core":
             LOGGER.critical("Core failure detected. Initiating SCRAM.")
@@ -766,6 +856,7 @@ class Reactor:
             "generator_auto": self.generator_auto,
             "primary_relief_open": self.primary_relief_open,
             "secondary_relief_open": self.secondary_relief_open,
+            "pressurizer_level": self.pressurizer_level,
             "maintenance_active": list(self.maintenance_active),
             "generators": [
                 {
@@ -801,6 +892,7 @@ class Reactor:
         self.generator_auto = metadata.get("generator_auto", self.generator_auto)
         self.primary_relief_open = metadata.get("primary_relief_open", self.primary_relief_open)
         self.secondary_relief_open = metadata.get("secondary_relief_open", self.secondary_relief_open)
+        self.pressurizer_level = metadata.get("pressurizer_level", self.pressurizer_level)
         maint = metadata.get("maintenance_active")
         if maint is not None:
             self.maintenance_active = set(maint)

src/reactor_sim/state.py

@@ -43,6 +43,8 @@ class CoolantLoopState:
     pressure: float  # MPa
     mass_flow_rate: float  # kg/s
     steam_quality: float  # fraction of vapor
+    inventory_kg: float = 0.0  # bulk mass of coolant
+    level: float = 1.0  # fraction full relative to nominal volume
 
     def average_temperature(self) -> float:
         return 0.5 * (self.temperature_in + self.temperature_out)

src/reactor_sim/thermal.py

@@ -15,7 +15,7 @@ LOGGER = logging.getLogger(__name__)
 
 def heat_transfer(primary: CoolantLoopState, secondary: CoolantLoopState, core_power_mw: float) -> float:
     """Return MW transferred to the secondary loop."""
-    if secondary.mass_flow_rate <= 0.0:
+    if primary.mass_flow_rate <= 0.0 or secondary.mass_flow_rate <= 0.0:
         return 0.0
     delta_t1 = max(1e-3, primary.temperature_out - secondary.temperature_in)
     delta_t2 = max(1e-3, primary.temperature_in - secondary.temperature_out)
@@ -26,7 +26,20 @@ def heat_transfer(primary: CoolantLoopState, secondary: CoolantLoopState, core_p
     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))
+    ua_limited = ua * lmtd
+
+    # Prevent the heat exchanger from over-transferring and inverting the outlet temperatures.
+    primary_capacity = primary.mass_flow_rate * constants.COOLANT_HEAT_CAPACITY
+    secondary_capacity = secondary.mass_flow_rate * constants.COOLANT_HEAT_CAPACITY
+    approach = 2.0  # K minimum approach between loop outlets
+    pinch_limited = 0.0
+    if primary_capacity > 0.0 and secondary_capacity > 0.0:
+        temp_gap = primary.temperature_out - secondary.temperature_in - approach
+        if temp_gap > 0.0:
+            pinch_watts = temp_gap / ((1.0 / primary_capacity) + (1.0 / secondary_capacity))
+            pinch_limited = max(0.0, pinch_watts / constants.MEGAWATT)
+
+    transferred = max(0.0, min(core_power_mw, ua_limited, pinch_limited))
     LOGGER.debug("Heat transfer %.2f MW with LMTD=%.1fK (ΔT1=%.1f ΔT2=%.1f)", transferred, lmtd, delta_t1, delta_t2)
     return transferred
 
@@ -51,11 +64,20 @@ def saturation_pressure(temp_k: float) -> float:
 class ThermalSolver:
     primary_volume_m3: float = 300.0
 
-    def step_core(self, core: CoreState, primary: CoolantLoopState, power_mw: float, dt: float) -> None:
+    def step_core(
+        self,
+        core: CoreState,
+        primary: CoolantLoopState,
+        power_mw: float,
+        dt: float,
+        residual_power_mw: float | None = None,
+    ) -> None:
+        if residual_power_mw is None:
+            residual_power_mw = power_mw
         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, power_mw - temp_rise) * dt
+        heating = 0.005 * 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.

src/reactor_sim/turbine.py

@@ -27,6 +27,7 @@ class Turbine:
     mechanical_efficiency: float = constants.STEAM_TURBINE_EFFICIENCY
     rated_output_mw: float = 400.0  # cap per unit electrical output
     spool_time: float = constants.TURBINE_SPOOL_TIME
+    throttle: float = 1.0  # 0-1 valve position
 
     def step(
         self,
@@ -46,11 +47,15 @@ class Turbine:
             state.condenser_temperature = max(305.0, loop.temperature_in - 20.0)
             return
 
+        throttle = min(constants.TURBINE_THROTTLE_MAX, max(constants.TURBINE_THROTTLE_MIN, self.throttle))
+        throttle_eff = 1.0 - constants.TURBINE_THROTTLE_EFFICIENCY_DROP * (constants.TURBINE_THROTTLE_MAX - throttle)
+
         enthalpy = 2_700.0 + loop.steam_quality * 600.0
-        mass_flow = effective_mass_flow * 0.6
+        mass_flow = effective_mass_flow * 0.6 * throttle
         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
+        backpressure_loss = 1.0 - _backpressure_penalty(loop)
+        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
@@ -71,5 +76,15 @@ class Turbine:
 def _ramp(current: float, target: float, dt: float, time_constant: float) -> float:
     if time_constant <= 0.0:
         return target
-    alpha = min(1.0, max(0.0, dt / time_constant))
+    alpha = min(1.0, max(0.0, dt / max(1e-6, time_constant)))
     return current + (target - current) * alpha
+
+
+def _backpressure_penalty(loop: CoolantLoopState) -> float:
+    base = constants.CONDENSER_BASE_PRESSURE_MPA
+    max_p = constants.CONDENSER_MAX_PRESSURE_MPA
+    pressure = max(base, min(max_p, loop.pressure))
+    if pressure <= base:
+        return 0.0
+    frac = (pressure - base) / max(1e-6, max_p - base)
+    return min(constants.CONDENSER_BACKPRESSURE_PENALTY, frac * constants.CONDENSER_BACKPRESSURE_PENALTY)

tests/test_aux_power.py

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

tests/test_control.py

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

tests/test_neutronics.py

@@ -1,3 +1,5 @@
+import pytest
+
 from reactor_sim.neutronics import NeutronDynamics
 from reactor_sim.state import CoreState
 
@@ -23,3 +25,21 @@ def test_reactivity_increases_with_rod_withdrawal():
     rho_full_out = dynamics.reactivity(state, control_fraction=0.0)
     rho_half = dynamics.reactivity(state, control_fraction=0.5)
     assert rho_full_out > rho_half
+
+
+def test_poisons_accumulate_under_power():
+    dynamics = NeutronDynamics()
+    state = _core_state(power=800.0, flux=1e6)
+    dynamics.update_poisons(state, dt=100.0)
+    dynamics.update_poisons(state, dt=100.0)
+    assert state.iodine_inventory > 0.0
+    assert state.xenon_inventory > 0.0
+
+
+def test_xenon_penalty_caps():
+    dynamics = NeutronDynamics()
+    state = _core_state()
+    state.xenon_inventory = 50.0
+    assert dynamics.xenon_penalty(state) == 0.05
+    state.xenon_inventory = 0.2
+    assert dynamics.xenon_penalty(state) == pytest.approx(0.01)

tests/test_pressurizer.py

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

tests/test_turbine.py

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