Refine relief venting and pump pressure caps

src/reactor_sim/reactor.py

@@ -281,6 +281,8 @@ class Reactor:
                 state.primary_loop.pressure, saturation_pressure(state.primary_loop.temperature_out), 0.1
             )
             target_pressure = max(0.5, base_head * power_ratio)
+            if self.primary_relief_open:
+                target_pressure = min(target_pressure, 1.0)
             primary_flow_scale = min(
                 self._inventory_flow_scale(state.primary_loop), self._npsh_factor(state.primary_loop)
             )
@@ -335,6 +337,8 @@ class Reactor:
             demand = 0.75
             base_flow, base_head = self.secondary_pump.performance(demand)
             target_pressure = max(0.5, base_head * power_ratio)
+            if self.secondary_relief_open:
+                target_pressure = min(target_pressure, 1.0)
             loop_pressure = max(
                 state.secondary_loop.pressure, saturation_pressure(state.secondary_loop.temperature_out), 0.1
             )
@@ -398,14 +402,25 @@ class Reactor:
         self.thermal.step_core(state.core, state.primary_loop, total_power, dt, residual_power_mw=residual)
         self.thermal.step_secondary(state.secondary_loop, transferred, dt)
         if self.primary_relief_open:
-            state.primary_loop.pressure = max(
-                0.1, min(state.primary_loop.pressure, saturation_pressure(state.primary_loop.temperature_out))
+            self._vent_relief(
+                state.primary_loop,
+                target_pressure=1.0,
+                vent_rate_max=0.02,
+                ramp_time=12.0,
+                dt=dt,
             )
+            for pump_state in state.primary_pumps:
+                pump_state.pressure = state.primary_loop.pressure
         if self.secondary_relief_open:
-            state.secondary_loop.pressure = max(
-                0.1, min(state.secondary_loop.pressure, saturation_pressure(state.secondary_loop.temperature_out))
+            self._vent_relief(
+                state.secondary_loop,
+                target_pressure=1.0,
+                vent_rate_max=0.05,
+                ramp_time=10.0,
+                dt=dt,
             )
-            self._vent_secondary_relief(state.secondary_loop, dt)
+            for pump_state in state.secondary_pumps:
+                pump_state.pressure = state.secondary_loop.pressure
         if not self.control.manual_control and not self.shutdown:
             self.control.safety_backoff(state.core.subcooling_margin, state.core.dnb_margin, dt)
         self._apply_secondary_boiloff(state, dt)
@@ -799,15 +814,21 @@ class Reactor:
             remaining = max(0.0, remaining - delivered)
         return total_power
 
-    def _vent_secondary_relief(self, loop: CoolantLoopState, dt: float) -> None:
+    def _vent_relief(
+        self,
+        loop: CoolantLoopState,
+        target_pressure: float,
+        vent_rate_max: float,
+        ramp_time: float,
+        dt: float,
+    ) -> None:
         """Model relief valve venting: gradual depressurization with mass/enthalpy loss."""
-        target_pressure = 1.0  # MPa surrogate relief sink (not full vacuum)
         if loop.inventory_kg <= 0.0:
             loop.pressure = max(target_pressure, loop.pressure)
             return
         # Vent rate scales with overpressure; capped to keep a multi-second depressurization.
         overfrac = max(0.0, (loop.pressure - target_pressure) / max(1e-6, loop.pressure))
-        vent_rate = min(0.05, 0.01 + 0.1 * overfrac)  # fraction of mass per second
+        vent_rate = min(vent_rate_max, 0.01 + vent_rate_max * overfrac)  # fraction of mass per second
         vent_mass = min(loop.inventory_kg, loop.inventory_kg * vent_rate * dt)
         if vent_mass > 0.0:
             specific_enthalpy = (
@@ -816,9 +837,18 @@ class Reactor:
             )
             loop.inventory_kg = max(0.0, loop.inventory_kg - vent_mass)
             loop.energy_j = max(0.0, loop.energy_j - vent_mass * specific_enthalpy)
-        # Pressure ramps toward target with a ~10s time constant.
-        ramp = min(1.0, dt / 10.0)
+        # Pressure ramps toward target with the requested time constant.
+        ramp = min(1.0, dt / max(1e-6, ramp_time))
         loop.pressure = max(target_pressure, loop.pressure - (loop.pressure - target_pressure) * ramp)
+        # Cool toward saturation at the new pressure to avoid re-pressurizing from superheat.
+        sat_target = saturation_temperature(target_pressure)
+        if loop.temperature_out > sat_target:
+            temp_drop = (loop.temperature_out - sat_target) * ramp
+            loop.temperature_out -= temp_drop
+            loop.temperature_in = min(loop.temperature_in, loop.temperature_out)
+            loop.steam_quality = 0.0
+            cp = constants.COOLANT_HEAT_CAPACITY
+            loop.energy_j = max(0.0, loop.inventory_kg * cp * loop.average_temperature())
         # Re-resolve temperature/quality/pressure to reflect the vented state.
         try:
             self.thermal._resolve_secondary_state(loop)  # type: ignore[attr-defined]