Merge pull request #54 from openmsr/batchwise_pr_keff_search_update

Batchwise pr keff search update
This commit is contained in:
Lorenzo Chierici 2026-01-19 14:50:01 +01:00 committed by GitHub
commit 1258e263b0
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
8 changed files with 411 additions and 1229 deletions

View file

@ -20,10 +20,10 @@ from warnings import warn
import numpy as np
from uncertainties import ufloat
from openmc.checkvalue import check_value, check_type, check_greater_than, PathLike
from openmc.checkvalue import check_type, check_greater_than, PathLike
from openmc.mpi import comm
from openmc.utility_funcs import change_directory
from openmc import Material, Cell
from openmc import Material
from .stepresult import StepResult
from .chain import _get_chain
from .results import Results, _SECONDS_PER_MINUTE, _SECONDS_PER_HOUR, \
@ -31,12 +31,7 @@ from .results import Results, _SECONDS_PER_MINUTE, _SECONDS_PER_HOUR, \
from .pool import deplete
from .reaction_rates import ReactionRates
from .transfer_rates import TransferRates, ExternalSourceRates
from .reactivity_control import (
ReactivityController,
CellReactivityController,
MaterialReactivityController
)
from .reactivity_control import ReactivityController
__all__ = [
@ -626,8 +621,6 @@ class Integrator(ABC):
reactivity_control : openmc.deplete.ReactivityController
Instance of ReactivityController class to perform reactivity control during
transport-depletion simulation.
Transfer rates for the depletion system used to model continuous
removal/feed between materials.
.. versionadded:: 0.15.4
@ -915,6 +908,7 @@ class Integrator(ABC):
n, res = self._get_bos_data_from_operator(i, source_rate, n)
else:
n, res = self._get_bos_data_from_restart(source_rate, n)
# Get root from reactivity control
if self._reactivity_control:
root = self.operator.prev_res[-1].reac_cont
else:
@ -1093,27 +1087,111 @@ class Integrator(ABC):
self.transfer_rates.set_redox(material, buffer, oxidation_states, timesteps)
def add_reactivity_control(
self,
obj: Union[Cell, Material],
**kwargs):
"""Add pre-defined Cell based or Material bases reactivity control to
integrator scheme.
self,
function: Callable,
x0: float,
x1: float,
bracket: list[float],
**kwargs
):
"""Add reactivity control to the integrator scheme.
This method creates a :class:`openmc.deplete.ReactivityController` that
performs keff searches during depletion to maintain a target reactivity
by adjusting a model parameter through the provided function.
.. important::
The function **must** modify the model through ``openmc.lib`` (e.g.,
``openmc.lib.cells``, ``openmc.lib.materials``) and **NOT** through
``openmc.model``. The function is called within a
:class:`openmc.lib.TemporarySession` context where only the C API
(``openmc.lib``) is available for modifications.
Parameters
----------
obj : openmc.Cell or openmc.Material
Cell or Material identifier to where add reactivity control
function : Callable
Function that modifies the model through ``openmc.lib`` based on a
parameter value. The function should take a single float parameter
and modify ``openmc.lib`` objects accordingly (e.g., adjust control
rod position via ``openmc.lib.cells[...].translation``, material
density via ``openmc.lib.materials[...].set_densities(...)``, etc.).
**Important**: The function must modify ``openmc.lib`` objects, not
``openmc.model`` objects.
x0: float
Initial lower bound for the keff search.
x1: float
Initial upper bound for the keff search.
bracket : list[float]
Bracket interval [x_min, x_max] that constrains the allowed parameter
values during the keff search. This is a required parameter
that defines the absolute bounds for the search. The bracket must contain
exactly 2 elements with bracket[0] < bracket[1]. These values are passed
directly to the ``x_min`` and ``x_max`` optional arguments in
:meth:`openmc.Model.keff_search`, which enforce hard limits on the
parameter range. If the keff search converges to a value outside this
bracket, it will be clamped to the nearest bracket bound with a warning.
**kwargs
keyword arguments that are passed to the specific ReactivityController
class.
Additional keyword arguments passed to
:meth:`openmc.Model.keff_search`. Common options include:
* ``target`` : float, optional
Target keff value to search for. Defaults to 1.0.
* ``k_tol`` : float, optional
Stopping criterion on the function value. Defaults to 1e-4.
* ``sigma_final`` : float, optional
Maximum accepted k-effective uncertainty. Defaults to 3e-4.
* ``maxiter`` : int, optional
Maximum number of iterations. Defaults to 50.
See :meth:`openmc.Model.keff_search` for a complete list of
available options.
Examples
--------
Add reactivity control that adjusts a control rod position:
>>> def adjust_rod_position(position):
... openmc.lib.cells[rod_cell.id].translation = [0, 0, position]
>>> integrator.add_reactivity_control(
... adjust_rod_position,
... x0=0.0,
... x1=5.0,
... bracket=[-10,10],
... target=1.0,
... k_tol=1e-4
... )
Add reactivity control that adjusts material density:
>>> def adjust_material_density(density_factor):
... # Get the material from openmc.lib
... lib_mat = openmc.lib.materials[material_id]
... # Get current nuclides and densities
... nuclides = lib_mat.nuclides
... current_densities = lib_mat.densities
... # Scale all densities by the factor
... new_densities = densities * density_factor
... # Update the material densities
... lib_mat.set_densities(nuclides, new_densities)
>>> integrator.add_reactivity_control(
... adjust_material_density,
... x0=0.8,
... x1=1.2,
... bracket=[0.2, 1.5],
... target=1.0
... )
.. versionadded:: 0.15.4
"""
if isinstance(obj, Cell):
reactivity_control = CellReactivityController
elif isinstance(obj, Material):
reactivity_control = MaterialReactivityController
self._reactivity_control = reactivity_control.from_params(obj,
self.operator, **kwargs)
self._reactivity_control = ReactivityController(
self.operator,
function,
x0,
x1,
bracket,
**kwargs)
@add_params
class SIIntegrator(Integrator):

File diff suppressed because it is too large Load diff

View file

@ -1,6 +1,5 @@
from collections.abc import Callable
from numbers import Real
from warnings import warn
import scipy.optimize as sopt
@ -13,7 +12,7 @@ _SCALAR_BRACKETED_METHODS = {'brentq', 'brenth', 'ridder', 'bisect'}
def _search_keff(guess, target, model_builder, model_args, print_iterations,
run_args, guesses, results, run_in_memory):
run_args, guesses, results):
"""Function which will actually create our model, run the calculation, and
obtain the result. This function will be passed to the root finding
algorithm
@ -40,8 +39,7 @@ def _search_keff(guess, target, model_builder, model_args, print_iterations,
results : Iterable of Real
Running list of results thus far, to be updated during the execution of
this function.
run_in_memory : bool
Whether or not to run the openmc model in memory.
Returns
-------
float
@ -53,11 +51,7 @@ def _search_keff(guess, target, model_builder, model_args, print_iterations,
model = model_builder(guess, **model_args)
# Run the model and obtain keff
if run_in_memory:
openmc.lib.run(**run_args)
sp_filepath = f'statepoint.{model.settings.batches}.h5'
else:
sp_filepath = model.run(**run_args)
sp_filepath = model.run(**run_args)
with openmc.StatePoint(sp_filepath) as sp:
keff = sp.keff
@ -76,7 +70,7 @@ def _search_keff(guess, target, model_builder, model_args, print_iterations,
def search_for_keff(model_builder, initial_guess=None, target=1.0,
bracket=None, model_args=None, tol=None,
bracketed_method='bisect', print_iterations=False,
run_args=None, run_in_memory=False, **kwargs):
run_args=None, **kwargs):
"""Function to perform a keff search by modifying a model parametrized by a
single independent variable.
@ -111,9 +105,6 @@ def search_for_keff(model_builder, initial_guess=None, target=1.0,
run_args : dict, optional
Keyword arguments to pass to :meth:`openmc.Model.run`. Defaults to no
arguments.
run_in_memory : bool
Whether or not to run the openmc model in memory.
Defaults to False.
.. versionadded:: 0.13.1
**kwargs
@ -202,26 +193,13 @@ def search_for_keff(model_builder, initial_guess=None, target=1.0,
# Add information to be passed to the searching function
args['args'] = (target, model_builder, model_args, print_iterations,
run_args, guesses, results, run_in_memory)
run_args, guesses, results)
# Create a new dictionary with the arguments from args and kwargs
args.update(kwargs)
# Perform the search
if not run_in_memory:
zero_value = root_finder(**args)
return zero_value, guesses, results
zero_value = root_finder(**args)
else:
try:
zero_value, root_res = root_finder(**args, full_output=True, disp=False)
if root_res.converged:
return zero_value, guesses, results
else:
warn(f'{root_res.flag}')
return guesses, results
# In case the root finder is not successful
except Exception as e:
warn(f'{e}')
return guesses, results
return zero_value, guesses, results

View file

@ -10,10 +10,6 @@ import numpy as np
import openmc
import openmc.lib
from openmc.deplete import CoupledOperator
from openmc.deplete import (
CellReactivityController,
MaterialReactivityController
)
from tests.regression_tests import config
@ -82,45 +78,37 @@ def model():
return openmc.Model(geometry, materials, settings)
@pytest.fixture
def mix_mat():
mix_mat = openmc.Material()
mix_mat.add_element("U", 1, percent_type="ao", enrichment=90)
mix_mat.add_element("O", 2)
mix_mat.set_density("g/cc", 10.4)
return mix_mat
def translate_cell(position):
cell_trans = [cell for cell in openmc.lib.cells.values() if cell.name == 'trans_cell'][0]
openmc.lib.cells[cell_trans.id].translation = [0, 0, position]
@pytest.mark.parametrize("obj, attribute, bracket, bracket_limit, axis, ref_result", [
('trans_cell', 'translation', [-5,5], [-40,40], 2, 'depletion_with_translation'),
('rot_cell', 'rotation', [-5,5], [-90,90], 2, 'depletion_with_rotation'),
('f', None, [-1,2], [-100,100], None, 'depletion_with_refuel')
def rotate_cell(angle):
cell_rot = [cell for cell in openmc.lib.cells.values() if cell.name == 'rot_cell'][0]
openmc.lib.cells[cell_rot.id].rotation = [0, 0, angle]
def adjust_material_density(density_factor):
f = [material for material in openmc.lib.materials.values() if material.name == 'f'][0]
nuclides = openmc.lib.materials[f.id].nuclides
current_densities = openmc.lib.materials[f.id].densities
new_densities = [d * density_factor for d in current_densities]
openmc.lib.materials[f.id].set_densities(nuclides, new_densities)
@pytest.mark.parametrize("function, x0, x1, bracket, ref_result", [
(translate_cell, -11, -5, [-15,0], 'depletion_with_translation'),
(rotate_cell, -80, -50, [-90,0], 'depletion_with_rotation'),
(adjust_material_density, 0.5, 1.2, [0.2, 1.5], 'depletion_with_refuel')
])
def test_reactivity_control(run_in_tmpdir, model, obj, attribute, bracket,
bracket_limit, axis, ref_result, mix_mat):
def test_reactivity_control(run_in_tmpdir, model, function, x0, x1, bracket, ref_result):
chain_file = Path(__file__).parents[2] / 'chain_simple.xml'
op = CoupledOperator(model, chain_file)
integrator = openmc.deplete.PredictorIntegrator(
op, [1], 174., timestep_units = 'd')
if attribute:
integrator.reactivity_control = CellReactivityController(
cell=obj,
operator=op,
attribute=attribute,
bracket=bracket,
bracket_limit=bracket_limit,
axis=axis
)
else:
integrator.reactivity_control = MaterialReactivityController(
material=obj,
operator=op,
material_to_mix=mix_mat,
bracket=bracket,
bracket_limit=bracket_limit,
)
integrator.add_reactivity_control(function, x0, x1, bracket,
kwargs={'output': True, 'k_tol': 0.1, 'sigma_final': 1e-3})
integrator.integrate()
# Get path to test and reference results

View file

@ -8,10 +8,7 @@ import numpy as np
import openmc
import openmc.lib
from openmc.deplete import CoupledOperator
from openmc.deplete import (
CellReactivityController,
MaterialReactivityController
)
from openmc.deplete import ReactivityController
CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml"
@ -73,150 +70,177 @@ def integrator(operator):
return openmc.deplete.PredictorIntegrator(
operator, [1,1], 0.0, timestep_units = 'd')
@pytest.fixture
def mix_mat():
mix_mat = openmc.Material()
mix_mat.add_element("U", 1, percent_type="ao", enrichment=90)
mix_mat.add_element("O", 2)
mix_mat.set_density("g/cc", 10.4)
return mix_mat
@pytest.mark.parametrize("case_name, obj, attribute, bracket, limit, axis", [
('cell translation','universe_cell', 'translation', [-1,1], [-10,10], 2),
('cell rotation', 'universe_cell', 'rotation', [-1,1], [-10,10], 2),
('material mix', 'fuel', None, [0,5], [-10,10], None),
])
def test_attributes(case_name, model, operator, integrator, obj, attribute,
bracket, limit, axis, mix_mat):
"""
Test classes attributes are set correctly
"""
if attribute:
integrator.reactivity_control = CellReactivityController(
cell=obj,
operator=operator,
attribute=attribute,
bracket=bracket,
bracket_limit=limit,
axis=axis
)
assert integrator.reactivity_control.depletable_cells == [cell for cell in \
model.geometry.get_cells_by_name(obj)[0].fill.cells.values() \
if cell.fill.depletable]
assert integrator.reactivity_control.axis == axis
else:
integrator.reactivity_control = MaterialReactivityController(
material=obj,
operator=operator,
material_to_mix=mix_mat,
bracket=bracket,
bracket_limit=limit
)
assert integrator.reactivity_control.material_to_mix == mix_mat
assert integrator.reactivity_control.bracket == bracket
assert integrator.reactivity_control.bracket_limit == limit
assert integrator.reactivity_control.burn_mats == operator.burnable_mats
assert integrator.reactivity_control.local_mats == operator.local_mats
@pytest.mark.parametrize("obj, attribute, value_to_set", [
('universe_cell', 'translation', 0),
('universe_cell', 'rotation', 0)
])
def test_cell_methods(run_in_tmpdir, model, operator, integrator, obj, attribute,
value_to_set):
"""
Test cell base class internal method
"""
integrator.reactivity_control = CellReactivityController(
cell=obj,
operator=operator,
attribute=attribute,
bracket=[-1,1],
bracket_limit=[-10,10],
axis=2
)
model.export_to_xml()
openmc.lib.init()
integrator.reactivity_control._set_lib_cell()
integrator.reactivity_control._set_cell_attrib(value_to_set)
assert integrator.reactivity_control._get_cell_attrib() == value_to_set
vol = integrator.reactivity_control._adjust_volumes()
integrator.reactivity_control._update_x_and_set_volumes(operator.number.number, vol)
for cell in integrator.reactivity_control.depletable_cells:
mat_id = str(cell.fill.id)
index_mat = operator.number.index_mat[mat_id]
assert vol[mat_id] == operator.number.volume[index_mat]
openmc.lib.finalize()
@pytest.mark.parametrize("units, value, dilute", [
('cc', 1.0, True),
('cm3', 1.0, False),
('grams', 1.0, True),
('grams', 1.0, False),
('atoms', 1.0, True),
('atoms', 1.0, False),
])
def test_internal_methods(run_in_tmpdir, model, operator, integrator, mix_mat,
units, value, dilute):
"""
Method to update volume in AtomNumber after depletion step.
Method inheritated by all derived classes so one check is enough.
"""
integrator.reactivity_control = MaterialReactivityController(
material='fuel',
operator=operator,
material_to_mix=mix_mat,
bracket=[-1,1],
bracket_limit=[-10,10],
units=units,
dilute=dilute
def test_reactivity_controller_init(operator):
"""Test ReactivityController initialization"""
def dummy_function(x):
return x
# Valid initialization
controller = ReactivityController(
operator=operator,
function=dummy_function,
x0=0.0,
x1=1.0,
bracket=[0.0, 2.0]
)
assert controller.operator == operator
assert controller.function == dummy_function
assert controller.x0 == 0.0
assert controller.x1 == 1.0
assert controller.kwargs['x_min'] == 0.0
assert controller.kwargs['x_max'] == 2.0
# Test invalid bracket with wrong length
with pytest.raises(ValueError, match="bracket must have exactly 2 elements"):
ReactivityController(operator, dummy_function, 0.0, 1.0, [0.0])
# Test invalid bracket with wrong order
with pytest.raises(ValueError, match="bracket\\[0\\] must be < bracket\\[1\\]"):
ReactivityController(operator, dummy_function, 0.0, 1.0, [2.0, 1.0])
def test_reactivity_controller_call(operator):
"""Test the __call__ method acts as a proxy to the function"""
call_log = []
def test_function(x):
call_log.append(x)
return x * 2
controller = ReactivityController(
operator=operator,
function=test_function,
x0=0.0,
x1=1.0,
bracket=[0.0, 2.0]
)
# Test that calling the controller calls the underlying function
result = controller(5.0)
assert result == 10.0
assert call_log == [5.0]
# Test with multiple arguments
def multi_arg_function(x, y, z=3):
return x + y + z
controller2 = ReactivityController(
operator=operator,
function=multi_arg_function,
x0=0.0,
x1=1.0,
bracket=[0.0, 2.0]
)
result = controller2(1, 2, z=4)
assert result == 7
def translate_cell(position):
"""Helper function to translate a cell"""
cell = [c for c in openmc.lib.cells.values() if c.name == 'universe_cell'][0]
openmc.lib.cells[cell.id].translation = [0, 0, position]
return position
def rotate_cell(angle):
"""Helper function to rotate a cell"""
cell = [c for c in openmc.lib.cells.values() if c.name == 'universe_cell'][0]
openmc.lib.cells[cell.id].rotation = [0, 0, angle]
return angle
def adjust_fuel_density(density_factor):
"""Helper function to adjust fuel density"""
fuel = [m for m in openmc.lib.materials.values() if m.name == 'fuel'][0]
nuclides = openmc.lib.materials[fuel.id].nuclides
current_densities = openmc.lib.materials[fuel.id].densities
new_densities = [d * density_factor for d in current_densities]
openmc.lib.materials[fuel.id].set_densities(nuclides, new_densities)
return density_factor
@pytest.mark.parametrize("function, x0, x1, bracket, test_value", [
(translate_cell, -1.0, 1.0, [-5.0, 5.0], 0.5),
(rotate_cell, -45.0, 45.0, [-90.0, 90.0], 10.0),
(adjust_fuel_density, 0.8, 1.2, [0.5, 1.5], 1.0)
])
def test_reactivity_controller_with_openmc_functions(
run_in_tmpdir, model, operator, function, x0, x1, bracket, test_value
):
"""Test ReactivityController with actual OpenMC geometry modification functions"""
controller = ReactivityController(
operator=operator,
function=function,
x0=x0,
x1=x1,
bracket=bracket,
kwargs={'output': False, 'k_tol': 0.1}
)
# Export model and initialize OpenMC library
model.export_to_xml()
openmc.lib.init()
try:
# Test that the controller can be called to modify the geometry
result = controller(test_value)
assert result == test_value
# Verify the function was actually called by checking it doesn't raise
# (actual verification would require checking geometry state)
finally:
openmc.lib.finalize()
# check material volume are adjusted correctly
mat = integrator.reactivity_control.material
mat_index = operator.number.index_mat[str(mat.id)]
nuc_index = operator.number.index_nuc['U235']
vol_before_mix = operator.number.get_mat_volume(str(mat.id))
def test_controller_kwargs_handling(operator):
"""Test that additional kwargs are properly stored and accessible"""
def dummy_func(x):
return x
custom_kwargs = {
'target': 1.0,
'k_tol': 1e-4,
'output': True
}
controller = ReactivityController(
operator=operator,
function=dummy_func,
x0=0.0,
x1=1.0,
bracket=[0.0, 2.0],
kwargs=custom_kwargs
)
# Check that bracket limits were added to kwargs
assert controller.kwargs['x_min'] == 0.0
assert controller.kwargs['x_max'] == 2.0
# Check that custom kwargs were preserved
assert controller.kwargs['target'] == 1.0
assert controller.kwargs['k_tol'] == 1e-4
assert controller.kwargs['output'] == True
# set mix_mat volume
integrator.reactivity_control._set_mix_material_volume(value)
mix_mat_vol = mix_mat.volume
#extract nuclide in atoms
mix_nuc_atoms = mix_mat.get_nuclide_atoms()['U235']
operator.number.number[mat_index][nuc_index] += mix_nuc_atoms
#update volume
vol_after_mix = integrator.reactivity_control._adjust_volumes(mix_mat_vol)
def test_integrator_add_reactivity_control(run_in_tmpdir, model, operator, integrator):
"""Test adding reactivity control to integrator using the add_reactivity_control method"""
def adjust_density(factor):
# This would modify material density in practice
return factor
# Test that add_reactivity_control method exists and works
integrator.add_reactivity_control(
function=adjust_density,
x0=0.8,
x1=1.2,
bracket=[0.5, 1.5],
kwargs={'k_tol': 0.1, 'output': False}
)
assert integrator.reactivity_control is not None
assert isinstance(integrator.reactivity_control, ReactivityController)
assert integrator.reactivity_control.x0 == 0.8
assert integrator.reactivity_control.x1 == 1.2
assert integrator.reactivity_control.kwargs['x_min'] == 0.5
assert integrator.reactivity_control.kwargs['x_max'] == 1.5
if dilute:
assert vol_before_mix + mix_mat_vol == vol_after_mix[str(mat.id)]
else:
assert vol_before_mix == vol_after_mix[str(mat.id)]
#check nuclide densities get assigned correctly in memory
x = [i[:operator.number.n_nuc_burn] for i in operator.number.number]
integrator.reactivity_control._update_materials(x)
nuc_index_lib = openmc.lib.materials[mat.id].nuclides.index('U235')
dens_to_compare = 1.0e-24 * operator.number.get_atom_density(str(mat.id), 'U235')
assert openmc.lib.materials[mat.id].densities[nuc_index_lib] == pytest.approx(dens_to_compare)
# check volumes get assigned correctly in AtomNumber
new_x = integrator.reactivity_control._update_x_and_set_volumes(x, vol_after_mix)
dens_to_compare = 1.0e24 * vol_after_mix[str(mat.id)] *\
openmc.lib.materials[mat.id].densities[nuc_index_lib]
assert new_x[mat_index][nuc_index] == pytest.approx(dens_to_compare)
# assert volume in AtomNumber is set correctly
assert operator.number.get_mat_volume(str(mat.id)) == vol_after_mix[str(mat.id)]
openmc.lib.finalize()