Restore full power when rods are withdrawn

src/reactor_sim/constants.py

@@ -11,6 +11,7 @@ COOLANT_DENSITY = 700.0  # kg/m^3 averaged between phases
 MAX_CORE_TEMPERATURE = 1_800.0  # K
 MAX_PRESSURE = 15.0  # MPa typical PWR primary loop limit
 CONTROL_ROD_SPEED = 0.03  # fraction insertion per second
+CONTROL_ROD_WORTH = 0.03  # delta rho contribution when fully withdrawn
 STEAM_TURBINE_EFFICIENCY = 0.34
 GENERATOR_EFFICIENCY = 0.96
 ENVIRONMENT_TEMPERATURE = 295.0  # K

src/reactor_sim/neutronics.py

@@ -33,7 +33,7 @@ class NeutronDynamics:
     def reactivity(self, state: CoreState, control_fraction: float) -> float:
         rho = (
             self.shutdown_bias +
-            0.02 * (1.0 - control_fraction)
+            constants.CONTROL_ROD_WORTH * (1.0 - control_fraction)
             + temperature_feedback(state.fuel_temperature)
             - fuel_reactivity_penalty(state.burnup)
             - xenon_poisoning(state.neutron_flux)

tests/test_simulation.py

@@ -144,3 +144,19 @@ def test_primary_pumps_spool_up_over_seconds():
         reactor.step(state, dt=1.0, command=ReactorCommand(coolant_demand=1.0))
 
     assert state.primary_loop.mass_flow_rate == pytest.approx(target_flow, rel=0.1)
+
+
+def test_full_rod_withdrawal_reaches_gigawatt_power():
+    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.secondary_pump_active = True
+
+    for _ in range(60):
+        reactor.step(state, dt=1.0)
+
+    assert state.core.power_output_mw > 2_000.0
+    assert state.core.fuel_temperature > 400.0