Add pressurizer and coolant inventory controls

FEATURES.md

@@ -4,6 +4,7 @@
 - **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 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, 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

@@ -35,6 +35,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/dashboard.py

@@ -384,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"),
             ],
         )
@@ -401,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"),

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(
@@ -256,12 +271,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 +293,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 +309,10 @@ 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
             )
+            pressurizer_floor = constants.PRIMARY_PRESSURIZER_SETPOINT_MPA - 0.5 if self.pressurizer_level > 0.05 else 0.5
+            target_pressure = max(pressurizer_floor, 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:
@@ -360,6 +388,10 @@ class Reactor:
             transferred = heat_transfer(state.primary_loop, state.secondary_loop, total_power)
         self.thermal.step_core(state.core, state.primary_loop, total_power, dt)
         self.thermal.step_secondary(state.secondary_loop, transferred)
+        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)
@@ -487,6 +519,59 @@ 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:
+        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.")
@@ -765,6 +850,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": [
                 {
@@ -800,6 +886,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)

tests/test_pressurizer.py

@@ -0,0 +1,35 @@
+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 = 5.0
+    reactor.pressurizer_level = 0.8
+    reactor.primary_pump_active = False
+    reactor.secondary_pump_active = False
+
+    reactor.step(state, dt=1.0, command=ReactorCommand.scram_all())
+
+    assert state.primary_loop.pressure > 5.0
+    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