Make keff search control work with depletion.

openmc/deplete/keff_search_control.py

@@ -54,6 +54,15 @@ def run(self, x):
         root : float
             Parameter value that achieves target keff
         """
+        # The keff search happens before the transport operator is called for
+        # this step, so both openmc.lib.materials and the operator's AtomNumber
+        # still hold the compositions from the previous operator call. Push the
+        # current beginning-of-step compositions in first, otherwise the search
+        # is performed on stale materials and _update_vec() below overwrites
+        # `x` with those stale densities, freezing the composition at its
+        # initial state for the entire depletion calculation.
+        self.operator._update_materials_and_nuclides(x)
+
         root = self._search_for_keff()
         self._update_vec(x)
         return root

tests/unit_tests/test_deplete_keff_search_control.py

@@ -1,4 +1,4 @@
-""" Tests for KeffSearchControl class """
+"""Tests for the KeffSearchControl class and openmc.deplete.keff_search_control."""
 
 from pathlib import Path
 
@@ -8,6 +8,7 @@
 import openmc
 import openmc.lib
 from openmc.deplete import CoupledOperator
+from openmc.deplete.keff_search_control import _KeffSearchControl
 
 CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml"
 
@@ -86,6 +87,35 @@ def set_u235_density(u235_density):
     return u235_density
 
 
+class MockOperator:
+    """Minimal operator recording calls to _update_materials_and_nuclides."""
+
+    def __init__(self, calls):
+        self.calls = calls
+
+    def _update_materials_and_nuclides(self, vec):
+        self.calls.append(('update_materials', [v.copy() for v in vec]))
+
+
+@pytest.fixture
+def control_and_calls(monkeypatch):
+    calls = []
+    operator = MockOperator(calls)
+    control = _KeffSearchControl(
+        operator, lambda x: None, x0=0.0, x1=1.0, bracket=[0.0, 2.0])
+
+    def fake_search():
+        calls.append(('search', None))
+        return 0.5
+
+    def fake_update_vec(x):
+        calls.append(('update_vec', None))
+
+    monkeypatch.setattr(control, '_search_for_keff', fake_search)
+    monkeypatch.setattr(control, '_update_vec', fake_update_vec)
+    return control, calls
+
+
 @pytest.mark.parametrize("function, x0, x1, bracket", [
     (translate_cell, -1.0, 1.0, (-5.0, 5.0)),
     (rotate_cell, -45.0, 45.0, (-90.0, 90.0)),
@@ -114,3 +144,28 @@ def test_integrator_add_keff_search_control(run_in_tmpdir, function, x0, x1, bra
     assert integrator._keff_search_control.search_kwargs['x_max'] == bracket[1]
     assert integrator._keff_search_control.search_kwargs['k_tol'] == 0.1
     assert not integrator._keff_search_control.search_kwargs['output']
+
+
+def test_materials_updated_before_search(control_and_calls):
+    """Compositions must be pushed to openmc.lib before the search runs.
+
+    The keff search is executed at the beginning of a depletion step, before
+    the transport operator is called. Without an explicit update, both
+    openmc.lib.materials and the operator's AtomNumber still hold the previous
+    call's compositions, and _update_vec() overwrites the depleted vector with
+    them -- freezing nuclide densities at their initial values.
+    """
+    control, calls = control_and_calls
+    n = [np.array([1.0, 2.0, 3.0]), np.array([4.0, 5.0, 6.0])]
+
+    root = control.run(n)
+
+    assert root == 0.5
+    assert [name for name, _ in calls] == [
+        'update_materials', 'search', 'update_vec']
+
+    # The vector handed to the operator must be the current composition
+    recorded = calls[0][1]
+    assert len(recorded) == len(n)
+    for actual, expected in zip(recorded, n):
+        np.testing.assert_array_equal(actual, expected)