Bias shutdown reactivity and keep cold start subcritical

src/reactor_sim/neutronics.py

@@ -28,9 +28,11 @@ class NeutronDynamics:
     beta_effective: float = 0.0065
     delayed_neutron_fraction: float = 0.0008
     external_source_coupling: float = 1e-6
+    shutdown_bias: float = -0.02
 
     def reactivity(self, state: CoreState, control_fraction: float) -> float:
         rho = (
+            self.shutdown_bias +
             0.02 * (1.0 - control_fraction)
             + temperature_feedback(state.fuel_temperature)
             - fuel_reactivity_penalty(state.burnup)
@@ -38,15 +40,22 @@ class NeutronDynamics:
         )
         return rho
 
-    def flux_derivative(self, state: CoreState, rho: float, external_source_rate: float = 0.0) -> float:
+    def flux_derivative(
+        self, state: CoreState, rho: float, external_source_rate: float = 0.0, baseline_source: float = 1e5
+    ) -> float:
         generation_time = constants.NEUTRON_LIFETIME
         beta = self.beta_effective
         source_term = self.external_source_coupling * external_source_rate
-        return ((rho - beta) / generation_time) * state.neutron_flux + 1e5 + source_term
+        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:
         rho = self.reactivity(state, control_fraction)
-        d_flux = self.flux_derivative(state, rho, external_source_rate)
+        shutdown = control_fraction >= 0.95
+        if shutdown:
+            rho = min(rho, -0.04)
+        baseline = 0.0 if shutdown else 1e5
+        source = 0.0 if shutdown else external_source_rate
+        d_flux = self.flux_derivative(state, rho, source, baseline_source=baseline)
         state.neutron_flux = max(0.0, state.neutron_flux + d_flux * dt)
         state.reactivity_margin = rho
         LOGGER.debug(