From 3e9733b4ff0ec00618fb65a38b1ea42e9f5668c5 Mon Sep 17 00:00:00 2001 From: church89 Date: Mon, 19 Jan 2026 11:51:03 +0100 Subject: [PATCH 1/4] refactoring of reactivity control --- openmc/deplete/abc.py | 122 +- openmc/deplete/reactivity_control.py | 1104 ++--------------- .../ref_depletion_with_refuel.h5 | Bin 39352 -> 28640 bytes .../ref_depletion_with_rotation.h5 | Bin 39352 -> 28640 bytes .../ref_depletion_with_translation.h5 | Bin 39352 -> 28640 bytes .../deplete_with_reactivity_control/test.py | 58 +- 6 files changed, 232 insertions(+), 1052 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 9108396f17..de76a891bc 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -31,11 +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 @@ -626,8 +622,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 @@ -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): diff --git a/openmc/deplete/reactivity_control.py b/openmc/deplete/reactivity_control.py index 08e86d9a6c..a70018375c 100644 --- a/openmc/deplete/reactivity_control.py +++ b/openmc/deplete/reactivity_control.py @@ -1,1031 +1,145 @@ -from abc import ABC, abstractmethod -from copy import deepcopy -from numbers import Real -from warnings import warn import numpy as np -from math import isclose - import openmc.lib from openmc.mpi import comm -from openmc.search import _SCALAR_BRACKETED_METHODS, search_for_keff -from openmc import Material, Cell -from openmc.data import atomic_mass, AVOGADRO -from openmc.checkvalue import ( - check_type, - check_value, - check_less_than, - check_greater_than, - check_iterable_type, - check_length, -) - - -class ReactivityController(ABC): - """Abstract class defining a generalized reactivity control. - - Reactivity control schemes, such as control rod adjustment or material - refuelling to control reactivity and maintain keff constant and equal to a - desired value, usually one. - - A reactivity control scheme can be added here to an integrator instance, - such as :class:`openmc.deplete.CECMIntegrator`, to parametrize one system - variable with the aim of satisfy certain design criteria, such as keeping - keff equal to one, while running transport-depletion calculations. - - This abstract class sets the requirements - for the reactivity control set. Users should instantiate - :class:`openmc.deplete.reactivity_control.CellReactivityController` - or - :class:`openmc.deplete.reactivity_control.MaterialReactivityController` - rather than this class. - .. versionadded:: 0.14.1 +from typing import Callable, Optional +from warnings import warn +class ReactivityController: + """Controller for reactivity adjustment during depletion calculations. + + This class performs keff searches to maintain a target reactivity by adjusting + a model parameter through a provided function. + Parameters ---------- operator : openmc.deplete.Operator - OpenMC operator object - bracket : list of float - Initial bracketing interval to search for the solution, relative to the - solution at previous step. - bracket_limit : list of float - Absolute bracketing interval lower and upper; if during the adaptive - algorithm the search_for_keff solution lies off these limits the closest + Depletion operator instance + function : Callable + Function that modifies the model based on a parameter value + x0 : float + Initial lower bound for the keff search + x1 : float + Initial upper bound for the keff search + bracket : list[float] + Absolute bracketing interval lower and upper + if keff search solution lies off these limits the closest limit will be set as new result. - bracketed_method : {'brentq', 'brenth', 'ridder', 'bisect'}, optional - Solution method to use. - This is equivalent to the `bracket_method` parameter of the - `search_for_keff`. - Defaults to 'brentq'. - tol : float - Tolerance for search_for_keff method. - This is equivalent to the `tol` parameter of the `search_for_keff`. - Default to 0.01 - target : Real, optional - This is equivalent to the `target` parameter of the `search_for_keff`. - Default to 1.0. - print_iterations : Bool, Optional - Print a status message each iteration - Default to True - search_for_keff_output : Bool, Optional - Print full transport logs during iterations (e.g., inactive generations, active - generations). Default to False - - Attributes - ---------- - burn_mats : list of str - List of burnable materials ids - local_mats : list of str - All burnable material IDs being managed by a single process - List of burnable materials ids - + kwargs : dict, optional + Additional keyword arguments to pass to `model.keff_search` """ - - def __init__( - self, - operator, - bracket, - bracket_limit, - bracketed_method="brentq", - tol=0.01, - target=1.0, - print_iterations=True, - search_for_keff_output=True, - ): - + def __init__(self, operator, function: Callable, x0: float, x1: float, bracket: list[float], kwargs: Optional[dict] = None): + if len(bracket) != 2: + raise ValueError(f"bracket must have exactly 2 elements, got {len(bracket)}") + if bracket[0] >= bracket[1]: + raise ValueError(f"bracket[0] must be < bracket[1], got {bracket}") + self.x0 = x0 + self.x1 = x1 self.operator = operator - self.burn_mats = operator.burnable_mats - self.local_mats = operator.local_mats - self.model = operator.model - self.geometry = operator.model.geometry - self.bracket = bracket + self.function = function + self.kwargs = {} if kwargs is None else kwargs + self.kwargs['x_min'] = bracket[0] + self.kwargs['x_max'] = bracket[1] + + def __call__(self, *args, **kwargs): + return self.function(*args, **kwargs) - check_iterable_type("bracket_limit", bracket_limit, Real) - check_length("bracket_limit", bracket_limit, 2) - check_less_than("bracket limit values", bracket_limit[0], bracket_limit[1]) - - self.bracket_limit = bracket_limit - self.bracketed_method = bracketed_method - self.tol = tol - self.target = target - self.print_iterations = print_iterations - self.search_for_keff_output = search_for_keff_output - - @property - def bracketed_method(self): - return self._bracketed_method - - @bracketed_method.setter - def bracketed_method(self, value): - check_value("bracketed_method", value, _SCALAR_BRACKETED_METHODS) - if value != "brentq": - warn("""brentq bracketed method is recommended to ensure - convergence of the method""") - self._bracketed_method = value - - @property - def tol(self): - return self._tol - - @tol.setter - def tol(self, value): - check_type("tol", value, Real) - self._tol = value - - @property - def target(self): - return self._target - - @target.setter - def target(self, value): - check_type("target", value, Real) - self._target = value - - @classmethod - def from_params(cls, obj, operator, **kwargs): - return cls(obj, operator, **kwargs) - - @abstractmethod - def _model_builder(self, param): - """Build the parametric model to be solved. - - Callable function which builds a model according to a passed - parameter. This function must return an openmc.model.Model object. - - Parameters - ---------- - param : parameter - model function variable - - Returns - ------- - openmc.model.Model - OpenMC parametric model - """ - - @abstractmethod def search_for_keff(self, x, step_index): - """Perform the criticality search on parametric material variable. - - The :meth:`openmc.search.search_for_keff` solution is then used to - calculate the new material volume and update the atoms concentrations. - + """Perform keff search and update the atom density vector. + Parameters ---------- x : list of numpy.ndarray - Total atoms concentrations - + Current atom density vector (atoms per material) + step_index : int + Current depletion step index + Returns ------- x : list of numpy.ndarray - Updated total atoms concentrations + Updated atom density vector root : float - Search_for_keff returned root value + Parameter value that achieves target keff """ - - @abstractmethod - def _adjust_volumes(self, root=None): - """Adjust volumes after criticality search when cell or material are - modified. - - Parameters - ---------- - root : float, Optional - :meth:`openmc.search.search_for_keff` volume root for - :class:`openmc.deplete.reactivity_control.CellReactivityController` - instances, in cm3 + root = self._search_for_keff() + x = self._update_vec(x) + return x, root + + def _search_for_keff(self) -> float: + """Perform the keff search using the model's keff_search method. + Returns ------- - volumes : dict - Dictionary of calculated volume, where key mat id and value - material volume, in cm3 + float + Parameter value that achieves target keff + + Raises + ------ + ValueError + If the keff search fails to converge """ - - def _search_for_keff(self, val): - """Perform the criticality search for a given parametric model. - - If the solution lies off the initial bracket, this method iteratively - adapts it until :meth:`openmc.search.search_for_keff` return a valid - solution. - It calculates the ratio between guessed and corresponding keffs values - as the proportional term, so that the bracket can be moved towards the - target. - A bracket limit poses the upper and lower boundaries to the adapting - bracket. If one limit is hit, the algoritm will stop and the closest - limit value will be set. - - Parameters - ---------- - val : float - Previous result value - - Returns - ------- - root : float - Estimated value of the variable parameter where keff is the - targeted value - """ - # make sure we don't modify original bracket - bracket = deepcopy(self.bracket) - - # Run until a search_for_keff root is found or out of limits - root = None - while root == None: - search = search_for_keff( - self._model_builder, - bracket=np.array(bracket) + val, - tol=self.tol, - bracketed_method=self.bracketed_method, - target=self.target, - print_iterations=self.print_iterations, - run_args={"output": self.search_for_keff_output}, - run_in_memory=True, + with openmc.lib.TemporarySession(model=self.operator.model) as session: + # Only pass the first 3 required args plus explicitly provided kwargs + result = self.operator.model.keff_search( + self.function, + self.x0, + self.x1, + **self.kwargs + ) + if not result.converged: + raise ValueError( + f"Search for keff failed to converge. " + f"Termination reason: {result.termination_reason}" + ) + + root = result.root + + # Check if root is outside the bracket bounds and give a warning + if root < self.kwargs['x_min']: + warn( + f"keff search result ({root:.6f}) is below the lower bracket bound " + f"({self.kwargs['x_min']:.6f}). Clamping to bracket lower bound.", + UserWarning ) - # if len(search) is 3 search_for_keff was successful - if len(search) == 3: - res, _, _ = search - - # Check if root is within bracket limits - if self.bracket_limit[0] < res < self.bracket_limit[1]: - root = res - - else: - # Set res with the closest limit and continue - arg_min = abs(np.array(self.bracket_limit) - res).argmin() - warn( - "Search_for_keff returned root out of " - "bracket limit. Set root to {:.2f} and continue.".format( - self.bracket_limit[arg_min] - ) - ) - root = self.bracket_limit[arg_min] - - elif len(search) == 2: - guesses, keffs = search - - # Check if all guesses are within bracket limits - if all( - self.bracket_limit[0] <= guess <= self.bracket_limit[1] - for guess in guesses - ): - # Simple method to iteratively adapt the bracket - warn( - "Search_for_keff returned values below or above " - "target. Trying to iteratively adapt bracket..." - ) - - # if the difference between keffs is smaller than 1 pcm, - # the slope calculation will be overshoot, so let's set the root - # to the closest guess value - if abs(np.diff(keffs)) < 1.0e-5: - arg_min = abs(self.target - np.array(keffs)).argmin() - warn( - "Difference between keff values is " - "smaller than 1 pcm. " - "Set root to guess value with " - "closest keff to target and continue..." - ) - root = guesses[arg_min] - - # Calculate slope as ratio of delta bracket and delta keffs - slope = abs(np.diff(bracket) / np.diff(keffs))[0].n - # Move the bracket closer to presumed keff root by using the - # slope - - # Two cases: both keffs are below or above target - if np.mean(keffs) < self.target: - # direction of moving bracket: +1 is up, -1 is down - if guesses[np.argmax(keffs)] > guesses[np.argmin(keffs)]: - dir = 1 - else: - dir = -1 - bracket[np.argmin(keffs)] = bracket[np.argmax(keffs)] - bracket[np.argmax(keffs)] += ( - slope * (self.target - max(keffs).n) * dir - ) - else: - if guesses[np.argmax(keffs)] > guesses[np.argmin(keffs)]: - dir = -1 - else: - dir = 1 - bracket[np.argmax(keffs)] = bracket[np.argmin(keffs)] - bracket[np.argmin(keffs)] += ( - slope * (min(keffs).n - self.target) * dir - ) - - # check if adapted bracket lies completely outside of limits - msg = ( - "WARNING: Adaptive iterative bracket {} went off " - "bracket limits. Set root to {:.2f} and continue." - ) - if all(np.array(bracket) + val <= self.bracket_limit[0]): - warn(msg.format(bracket, self.bracket_limit[0])) - root = self.bracket_limit[0] - - if all(np.array(bracket) + val >= self.bracket_limit[1]): - warn(msg.format(bracket, self.bracket_limit[1])) - root = self.bracket_limit[1] - - # check if adapted bracket ends are outside bracketing limits - if bracket[1] + val > self.bracket_limit[1]: - bracket[1] = self.bracket_limit[1] - val - - if bracket[0] + val < self.bracket_limit[0]: - bracket[0] = self.bracket_limit[0] - val - - else: - # Set res with closest limit and continue - arg_min = abs(np.array(self.bracket_limit) - guesses).argmin() - warn( - "Search_for_keff returned values off " - "bracket limits. Set root to {:.2f} and continue.".format( - self.bracket_limit[arg_min] - ) - ) - root = self.bracket_limit[arg_min] - - else: - raise ValueError("ERROR: Search_for_keff output is not valid") + elif root > self.kwargs['x_max']: + warn( + f"keff search result ({root:.6f}) is above the upper bracket bound " + f"({self.kwargs['x_max']:.6f}). Clamping to bracket upper bound.", + UserWarning + ) + + #restore the number of initial batches + openmc.lib.settings.set_batches(self.operator.model.settings.batches) return root - def _update_materials(self, x): - """Update number density and material compositions in OpenMC on all - processes. - + def _update_vec(self, x): + """Update the atom density vector from the operator's number object. + + This method synchronizes the atom densities across all MPI ranks by + broadcasting the number object from each rank and updating the x vector + with the current atom densities and volumes. + Parameters ---------- x : list of numpy.ndarray - Total atom concentrations + Atom density vector to update (atoms per material) + + Returns + ------- + list of numpy.ndarray + Updated atom density vector synchronized across all ranks """ - self.operator.number.set_density(x) - for rank in range(comm.size): number_i = comm.bcast(self.operator.number, root=rank) - - for mat in number_i.materials: - nuclides = [] - densities = [] - + + for mat_idx, mat in enumerate(number_i.materials): for nuc in number_i.nuclides: - # get atom density in atoms/b-cm - val = 1.0e-24 * number_i.get_atom_density(mat, nuc) - if nuc in self.operator.nuclides_with_data: - if val > 0.0: - nuclides.append(nuc) - densities.append(val) - - # Update densities on C API side - openmc.lib.materials[int(mat)].set_densities(nuclides, densities) - - def _update_x_and_set_volumes(self, x, volumes): - """Update x vector with new volumes, before assign them to AtomNumber. - - Parameters - ---------- - x : list of numpy.ndarray - Total atoms concentrations - volumes : dict - Updated volumes, where key is material id and value material - volume, in cm3 - - Returns - ------- - x : list of numpy.ndarray - Updated total atoms concentrations - """ - number_i = self.operator.number - - for mat_idx, mat in enumerate(self.local_mats): - - if mat in volumes: - res_vol = volumes[mat] - # Normalize burnable nuclides in x vector without cross section data - for nuc_idx, nuc in enumerate(number_i.burnable_nuclides): - if nuc not in self.operator.nuclides_with_data: - # normalize with new volume - x[mat_idx][nuc_idx] *= number_i.get_mat_volume(mat) / res_vol - - # update all concentration data with the new updated volumes - for nuc, dens in zip( - openmc.lib.materials[int(mat)].nuclides, - openmc.lib.materials[int(mat)].densities, - ): - nuc_idx = number_i.index_nuc[nuc] - n_of_atoms = dens / 1.0e-24 * res_vol - if nuc in number_i.burnable_nuclides: - # convert [#atoms/b-cm] into [#atoms] - x[mat_idx][nuc_idx] = n_of_atoms - # when the nuclide is not in depletion chain update the AtomNumber - else: - # Atom density needs to be in [#atoms/cm3] - number_i[mat, nuc] = n_of_atoms - - # Now we can update the new volume in AtomNumber - number_i.volume[mat_idx] = res_vol - return x - -class CellReactivityController(ReactivityController): - """Abstract class holding reactivity control cell-based functions. - - .. versionadded:: 0.14.1 - - Parameters - ---------- - cell : openmc.Cell or int or str - OpenMC Cell identifier to where apply a reactivty control - operator : openmc.deplete.Operator - OpenMC operator object - attribute : str - openmc.lib.cell attribute. Only support 'translation' and 'rotation' - bracket : list of float - Initial bracketing interval to search for the solution, relative to the - solution at previous step. - bracket_limit : list of float - Absolute bracketing interval lower and upper; if during the adaptive - algorithm the search_for_keff solution lies off these limits the closest - limit will be set as new result. - bracketed_method : {'brentq', 'brenth', 'ridder', 'bisect'}, optional - Solution method to use. - This is equivalent to the `bracket_method` parameter of the - `search_for_keff`. - Defaults to 'brentq'. - tol : float - Tolerance for search_for_keff method. - This is equivalent to the `tol` parameter of the `search_for_keff`. - Default to 0.01 - target : Real, optional - This is equivalent to the `target` parameter of the `search_for_keff`. - Default to 1.0. - print_iterations : Bool, Optional - Whether or not to print `search_for_keff` iterations. - Default to True - search_for_keff_output : Bool, Optional - Whether or not to print transport iterations during `search_for_keff`. - Default to False - - Attributes - ---------- - cell : openmc.Cell or int or str - OpenMC Cell identifier to where apply batch wise scheme - attribute : str - openmc.lib.cell attribute. Only support 'translation' and 'rotation' - depletable_cells : list of openmc.Cell - Cells that fill the openmc.Universe that fills the main cell - to where apply batch wise scheme, if cell materials are set as - depletable. - lib_cell : openmc.lib.Cell - Corresponding openmc.lib.Cell once openmc.lib is initialized - axis : int {0,1,2} - Directional axis for geometrical parametrization, where 0, 1 and 2 stand - for 'x', 'y' and 'z', respectively. - - """ - - def __init__( - self, - cell, - operator, - attribute, - bracket, - bracket_limit, - axis, - bracketed_method="brentq", - tol=0.01, - target=1.0, - samples=1000000, - print_iterations=True, - search_for_keff_output=True, - ): - - super().__init__( - operator, - bracket, - bracket_limit, - bracketed_method, - tol, - target, - print_iterations, - search_for_keff_output, - ) - - self.cell = self._get_cell(cell) - - # Only lib cell attributes are valid - check_value("attribute", attribute, ('translation', 'rotation')) - self.attribute = attribute - - # Initialize to None, as openmc.lib is not initialized yet here - self.lib_cell = None - - # A cell translation or roation must correspond to a geometrical - # variation of its filled materials, thus at least 2 are necessary - if isinstance(self.cell.fill, openmc.Universe): - # check if cell fill universe contains at least 2 cells - check_greater_than("universe cells", - len(self.cell.fill.cells), 2, True) - - if isinstance(self.cell.fill, openmc.Lattice): - # check if cell fill lattice contains at least 2 universes - check_greater_than("lattice universes", - len(self.cell.fill.get_unique_universes()), 2, True) - - self.depletable_cells = [ - cell - for cell in self.cell.fill.cells.values() - if cell.fill.depletable - ] - - # index of cell directional axis - check_value("axis", axis, (0, 1, 2)) - self.axis = axis - - check_type("samples", samples, int) - self.samples = samples - - if self.depletable_cells: - self._initialize_volume_calc() - - def _get_cell_attrib(self): - """Get cell attribute coefficient. - - Returns - ------- - coeff : float - cell coefficient - """ - return getattr(self.lib_cell, self.attribute)[self.axis] - - def _set_cell_attrib(self, val): - """Set cell attribute to the cell instance. - - Parameters - ---------- - val : float - cell coefficient to set - """ - attr = getattr(self.lib_cell, self.attribute) - attr[self.axis] = val - setattr(self.lib_cell, self.attribute, attr) - - def _set_lib_cell(self): - """Set openmc.lib.cell cell to self.lib_cell attribute - """ - self.lib_cell = [ - cell for cell in openmc.lib.cells.values() if cell.id == self.cell.id - ][0] - - def _get_cell(self, val): - """Helper method for getting cell from cell instance or cell name or id. - - Parameters - ---------- - val : Openmc.Cell or str or int representing Cell - - Returns - ------- - val : openmc.Cell - Openmc Cell - """ - cell_bundles = [ - (cell.id, cell.name, cell) - for cell in self.geometry.get_all_cells().values() - ] - - if isinstance(val, Cell): - check_value( - "Cell exists", val, [cell_bundle[2] for cell_bundle in cell_bundles] - ) - - elif isinstance(val, str): - if val.isnumeric(): - check_value( - "Cell id exists", - int(val), - [cell_bundle[0] for cell_bundle in cell_bundles], - ) - - val = [ - cell_bundle[2] - for cell_bundle in cell_bundles - if cell_bundle[0] == int(val) - ][0] - - else: - check_value( - "Cell name exists", - val, - [cell_bundle[1] for cell_bundle in cell_bundles], - ) - - val = [ - cell_bundle[2] - for cell_bundle in cell_bundles - if cell_bundle[1] == val - ][0] - - elif isinstance(val, int): - check_value( - "Cell id exists", val, [cell_bundle[0] for cell_bundle in cell_bundles] - ) - - val = [ - cell_bundle[2] for cell_bundle in cell_bundles if cell_bundle[0] == val - ][0] - - else: - ValueError(f"Cell: {val} is not supported") - - return val - - def _initialize_volume_calc(self): - """Set volume calculation model settings of depletable materials filling - the parametric Cell. - """ - ll, ur = self.geometry.bounding_box - mat_vol = openmc.VolumeCalculation(self.depletable_cells, self.samples, ll, ur) - self.model.settings.volume_calculations = mat_vol - - def _adjust_volumes(self): - """Perform stochastic volume calculation and return new volumes. - - Returns - ------- - volumes : dict - Dictionary of calculated volumes, where key is mat id and value - material volume, in cm3 - """ - openmc.lib.calculate_volumes() - volumes = {} - if comm.rank == 0: - res = openmc.VolumeCalculation.from_hdf5("volume_1.h5") - for cell in self.depletable_cells: - mat_id = cell.fill.id - volumes[str(mat_id)] = res.volumes[cell.id].n - volumes = comm.bcast(volumes) - return volumes - - def _model_builder(self, param): - """Builds the parametric model to be passed to the - :meth:`openmc.search.search_for_keff` method. - - Callable function which builds a model according to a passed - parameter. This function must return an openmc.model.Model object. - - Parameters - ---------- - param : model parametric variable - cell translation coefficient or cell rotation coefficient - - Returns - ------- - self.model : openmc.model.Model - Openmc parametric model - """ - self._set_cell_attrib(param) - - return self.model - - def search_for_keff(self, x, step_index): - """Perform the criticality search on the parametric cell coefficient and - update materials accordingly. - - The :meth:`openmc.search.search_for_keff` solution is then set as the - new cell attribute. - - Parameters - ---------- - x : list of numpy.ndarray - Total atoms concentrations - - Returns - ------- - x : list of numpy.ndarray - Updated total atoms concentrations - root : float - Search_for_keff returned root value - """ - # set _cell argument, once openmc.lib is initialized - if self.lib_cell is None: - self._set_lib_cell() - - # Get cell attribute from previous iteration - val = self._get_cell_attrib() - check_type("Cell coeff", val, Real) - - # Update all material densities from concentration vectors - # before performing the search_for_keff. This is needed before running - # the transport equations in the search_for_keff algorithm - self._update_materials(x) - - # Calculate new cell attribute - root = self._search_for_keff(val) - - # set results value as attribute in the geometry - self._set_cell_attrib(root) - - # if at least one of the cell fill materials is depletable, assign new - # volume to the material and update x and number accordingly - # new volume - if self.depletable_cells: - volumes = self._adjust_volumes() - x = self._update_x_and_set_volumes(x, volumes) - - return x, root - -class MaterialReactivityController(ReactivityController): - """Abstract class holding reactivity control material-based functions. - - .. versionadded:: 0.14.1 - - Parameters - ---------- - material : openmc.Material or int or str - OpenMC Material identifier to where apply reactivity control - operator : openmc.deplete.Operator - OpenMC operator object - material_to_mix : openmc.Material - OpenMC Material to mix with main material for reactivity control. - bracket : list of float - Initial bracketing interval to search for the solution, relative to the - solution at previous step. - bracket_limit : list of float - Absolute bracketing interval lower and upper; if during the adaptive - algorithm the search_for_keff solution lies off these limits the closest - limit will be set as new result. - units : {'grams', 'atoms', 'cc', 'cm3'} str - Units of material parameter to be mixed. - Default to 'cc' - dilute : bool - Whether or not to update material volume after a material mix. - Default to False - bracketed_method : {'brentq', 'brenth', 'ridder', 'bisect'}, optional - Solution method to use. - This is equivalent to the `bracket_method` parameter of the - `search_for_keff`. - Defaults to 'brentq'. - tol : float - Tolerance for search_for_keff method. - This is equivalent to the `tol` parameter of the `search_for_keff`. - Default to 0.01 - target : Real, optional - This is equivalent to the `target` parameter of the `search_for_keff`. - Default to 1.0. - print_iterations : Bool, Optional - Whether or not to print `search_for_keff` iterations. - Default to True - search_for_keff_output : Bool, Optional - Whether or not to print transport iterations during `search_for_keff`. - Default to False - - Attributes - ---------- - material : openmc.Material or int or str - OpenMC Material identifier to where apply batch wise scheme - material_to_mix : openmc.Material - OpenMC Material to mix with main material for reactivity control. - units : {'grams', 'atoms', 'cc', 'cm3'} str - Units of material parameter to be mixed. - Default to 'cc' - dilute : bool - Whether or not to update material volume after a material mix. - Default to False - - """ - - def __init__( - self, - material, - operator, - material_to_mix, - bracket, - bracket_limit, - units='cc', - dilute=False, - bracketed_method="brentq", - tol=0.01, - target=1.0, - print_iterations=True, - search_for_keff_output=True, - ): - - super().__init__( - operator, - bracket, - bracket_limit, - bracketed_method, - tol, - target, - print_iterations, - search_for_keff_output, - ) - - self.material = self._get_material(material) - - check_type("material to mix", material_to_mix, openmc.Material) - for nuc in material_to_mix.get_nuclides(): - check_value("check nuclide exists", nuc, self.operator.nuclides_with_data) - self.material_to_mix = material_to_mix - - check_value('units', units, ('grams', 'atoms', 'cc', 'cm3')) - self.units = units - - check_type("dilute", dilute, bool) - self.dilute = dilute - - def _get_material(self, val): - """Helper method for getting openmc material from Material instance or - material name or id. - - Parameters - ---------- - val : Openmc.Material or str or int representing material name or id - - Returns - ------- - val : openmc.Material - Openmc Material - """ - if isinstance(val, Material): - check_value("Material", str(val.id), self.burn_mats) - - elif isinstance(val, str): - if val.isnumeric(): - check_value("Material id", val, self.burn_mats) - val = [mat for mat in self.model.materials if mat.id == int(val)][0] - - else: - check_value( - "Material name", - val, - [mat.name for mat in self.model.materials if mat.depletable], - ) - val = [mat for mat in self.model.materials if mat.name == val][0] - - elif isinstance(val, int): - check_value("Material id", str(val), self.burn_mats) - val = [mat for mat in self.model.materials if mat.id == val][0] - - else: - ValueError(f"Material: {val} is not supported") - - return val - - def search_for_keff(self, x, step_index): - """Perform the criticality search on parametric material variable. - - The :meth:`openmc.search.search_for_keff` solution is then used to - calculate the new material volume and update the atoms concentrations. - - Parameters - ---------- - x : list of numpy.ndarray - Total atoms concentrations - - Returns - ------- - x : list of numpy.ndarray - Updated total atoms concentrations - root : float - Search_for_keff returned root value - """ - # Update AtomNumber with new conc vectors x. Materials are also updated - # even though they are also re-calculated when running the search_for_kef - self._update_materials(x) - - # Set initial material addition to 0 and let program calculate the - # right amount - root = self._search_for_keff(0) - - # Update concentration vector and volumes with new value for depletion - volumes = self._adjust_volumes(root) - x = self._update_x_and_set_volumes(x, volumes) - - return x, root - - def _set_mix_material_volume(self, param): - """ - Set volume in cc to `self.material_to_mix` as a function of `param` after - converion, based on `self.units` attribute. - - Parameters - ---------- - param : float - grams or atoms or cc of material_to_mix to convert to cc - - """ - if self.units == 'grams': - multiplier = 1 / self.material_to_mix.get_mass_density() - elif self.units == 'atoms': - multiplier = ( - self.material_to_mix.average_molar_mass - / AVOGADRO - / self.material_to_mix.get_mass_density() - ) - elif self.units == 'cc' or self.units == 'cm3': - multiplier = 1 - - self.material_to_mix.volume = param * multiplier - - - def _model_builder(self, param): - """Callable function which builds a model according to a passed - parameter. This function must return an openmc.model.Model object. - - Builds the parametric model to be passed to the - :meth:`openmc.search.search_for_keff` method. - - Parameters - ---------- - param : float - Model function variable, cc of material_to_mix - - Returns - ------- - model : openmc.model.Model - Openmc parametric model - """ - for rank in range(comm.size): - number_i = comm.bcast(self.operator.number, root=rank) - - for mat_id in number_i.materials: - nuclides = [] - densities = [] - - if int(mat_id) == self.material.id: - self._set_mix_material_volume(param) - - if self.dilute: - # increase the total volume of the material - # assuming ideal materials mixing - vol = ( - number_i.get_mat_volume(mat_id) - + self.material_to_mix.volume - ) - else: - # we assume the volume after mix won't change - vol = number_i.get_mat_volume(mat_id) - - # Total atom concentration in [#atoms/cm-b] - #mat_index = number_i.index_mat[mat_id] - #tot_conc = 1.0e-24 * sum(number_i.number[mat_index]) / vol - - for nuc in number_i.index_nuc: - # check only nuclides with cross sections data - if nuc in self.operator.nuclides_with_data: - atoms = number_i[mat_id, nuc] - if nuc in self.material_to_mix.get_nuclides(): - # update atoms number - atoms += self.material_to_mix.get_nuclide_atoms()[nuc] - - atoms_per_bcm = 1.0e-24 * atoms / vol - - if atoms_per_bcm > 0.0: - nuclides.append(nuc) - densities.append(atoms_per_bcm) - - else: - # for all other materials, still check atom density limits - for nuc in number_i.nuclides: - if nuc in self.operator.nuclides_with_data: - # get normalized atoms density in [atoms/b-cm] - atoms_per_bcm = 1.0e-24 * number_i.get_atom_density(mat_id, nuc) - - if atoms_per_bcm > 0.0: - nuclides.append(nuc) - densities.append(atoms_per_bcm) - - # assign nuclides and densities to the in-memory model - openmc.lib.materials[int(mat_id)].set_densities(nuclides, densities) - - # always need to return a model - return self.model - - def _adjust_volumes(self, res): - """Uses :meth:`openmc.deplete.reactivity_control._search_for_keff` - solution as cc to update depletable volume if `self.dilute` attribute - is set to True. - - Parameters - ---------- - res : float - Solution in cc of material_to_mix, coming from - :meth:`openmc.deplete.reactivity_control._search_for_keff` - - Returns - ------- - volumes : dict - Dictionary of calculated volume, where key mat id and value - material volume, in cm3 - """ - number_i = self.operator.number - volumes = {} - - for mat_id in self.local_mats: - if int(mat_id) == self.material.id: - if self.dilute: - volumes[mat_id] = number_i.get_mat_volume(mat_id) + res - else: - volumes[mat_id] = number_i.get_mat_volume(mat_id) - return volumes + nuc_idx = number_i.burnable_nuclides.index(nuc) + atom_density = number_i.get_atom_density(mat, nuc) + volume = number_i.get_mat_volume(mat) + x[mat_idx][nuc_idx] = atom_density * volume + + x = comm.bcast(x, root=rank) + return x \ No newline at end of file diff --git a/tests/regression_tests/deplete_with_reactivity_control/ref_depletion_with_refuel.h5 b/tests/regression_tests/deplete_with_reactivity_control/ref_depletion_with_refuel.h5 index 83a19702b0ef95740453e39a79ae1c44bf364e2f..a508e3d980f2442786be8c1785d0138326eb37b9 100644 GIT binary patch delta 1405 zcmZ{jZ%7ky7{~8!GVd;{sWi9Y++mql+M2d(Ah??@9U{4cZ1}=P)MhVosR$xu^G&^M z*Yif9U{)_IC?;1j5^|$IL=lV#T5oCufe_TR7fI2dyZimlHAV;b=lgv=&-47AAHPdV zRPF;A%k_0HB8!O?n{;=&Av z&3bs)ctmZ!(NM(28`Pq~oyd(gs6`V^sHkaA5uD;xIJz6TabAV}9SC;Vi{M-nSnQWz z!R~<9o&d!3^hF*4#0K5~ZQBbYAwWB1CJt$h6R-~ASzRbhxbW&Vt}qP(1oY}60(j*b zSeOJ6+IgMCb8B22>r|7No|bl%OEf^|2(K@(+Nu=O!e&2TNno zAtfeeXow1&eq`d@sWW5nOgILqwpJMQ21+v$Avt52^ug!=PmiGHKB)9HGQkK^1dQtB zk)_J<6-IK-T3v!zE>xHJG_tF8N;{x*S$ZN7j7ly*S$-rH-A?2Ywl~0 z?}l1mjzf6Xy2aSH=Kd7u{`J|Nhxx-cF;n&5^r`806(-A?dw43ewDsu(T$!_0DCxq? KofrMlb@yLvbt82E literal 39352 zcmeHPeQXrR6`#Gc@evF;3T>_w=(*&>G!EFov=^$44{YF2D?(GK4n-<FkTcXScv^q?bqcZVp-ILPn4Uzdyf(CZSUz8exFGrD<=&3 z8ZZOtVD@l0h%XeD3zQdg&GJ6c+TLj-9*P+#E*6=B3(W}0w-_-~z_b{Q8@+vRVXBjB6n$hRbT1Z)e{g2UN@0P!t~4W>~cg#CwbsKE>Bm%s|jEj;|n z*83VXQ7ZbeXtW`>4$CQ+B9-fN%K;~aID;|F4FN>FR%|mV+6J z0Qsd9)-Q@-|0=j}8vjyjAb$wuo%&<676ke~hk|MOiPy8CMIb5hPCQHD_(^4p#WQuhw)x?XXC+u+TaP=Q zksc?W#g>5~<{9-64s#*k#IpeQBPv@go~iR)n;-6YM)kGzxZ@e=apKuX5G0sqgcOIl z5OCtzfnav2V(~1z9PBirJ)0lycs3WSZR>HzGt%S4v+%88hOkzK4< zJR7P2JHJ7DHb30)j9xR_dff4h^f>Ws=nkki^DMicaBN}(oOrek`w^8Y7SH^*frVGrr^QHP#z>GnC>oZheRnH9$~?vOLw3-m6O4={VrM=^TM>K?HDJN%qurCelzG zr~vCEr_PhDMwlY!alJ=!>Np?9Mr^|hkeuqLA}BYE<1NXJVqD;TrS~%<_&p=Z1@ME; z5Pm0t>ZUEA_c>~PY`^Gk9z*rF^|`5 z_wM*Z`x0zD?)XD`fIrUr3sTaln+JV20x>?%76erO1hL^sIf!Q*mEBdlo4UH1R$^Lf zd#|<6hImEeXGW2nmzwZq?T?}TXk+W381a<-e@+-@eWu!feTYaD{Tcy8{gnC@rMnGb z9_mSdzjqP&LUsaj@@yU<#lJ?r0O=$b!}y>M>9g&e9!C|t+`<2<& z(+X>^lV3?GJS84ezbeeicxmE;>K6*-1Jz6J;j^wWL8JY$C6GY#L|Y!}MB902|F=`V z*PJ=ZJg0pprKX3J`5qYV%~;Sf;PsTA-fUf2*-}LQq*-*~(}gG)q4Rh=;0c9RD$UEC z?{GlOevUzSP<}pKNUm^v;K72|GspK;yDFLZNYymVnU{smyg;dX?lM1Tj<+{;_F_hN zUvp=Bi$T%!ch8gnpHho-dN6J(Ff)zz-A4DGrp`VrnrN51l7#5(>5DfTyQGw5!QM9w z67sw1UClk+iCyuYo`j9Nx-J=i|DgSvcjEg;_d&g1lwSD#30n2{dP<_L)>XRyD)*%k za*KQ(j|cQ0&6AA>wU=<9c;W4bJL2_tlX1KK8Ex+`a*K1TpDZ54Z^DD(Q^4cQJ{u3I z4mmp0@PPZLDgTE+IJj@W$My&{<5| zlc0~vkzdN!%E4r|1eH@mn#cIE#W}_ZFanGKBftnS0*nA7zz8q`i~u9>5hLKVjy&PA z9>chvr1G>LoUFRh0@s!FT#?t8lSP>27y(9r5nu!u0Y-okU<4QeMt~7u1Q>y-BjEHr z!9PjWQC$C%e7@%g2aaVcCV2|^`9L;49Ag9+0Y-okU<4QeMt~7u1ZENfZr2ZOlUz3p zI;z*2sq6f^L-cAkMl;7qvyt|#S;$mMv|xUYB`QG z0*nA7zz8q`i~u9R2+S-53Y|9|$TMFYa+oLb{BUO7cdiX1zz8q`i~u802m*!97Y6do z7g9W5CFWxn#SH0Q&hR1WU?(*o4KY#E}Ec;(91DC>9J5C!` zSa{_{k`AGY42AI-q&fg_Km{`yNB zr*+Y{fA#4a8>iM|uPqC`9ZI)*N5yN;eK_C7>9z+uPbPk1LVY@Ggm>Rj#4 zp|kokS62?i7ah_Mzx>Q2$A@3h%YN59(6Rbm{q(_A>uR4rtna^a!IL=BsL%cB;);}s w(~2|wRnOfqrdNFYg4T8QUv=}p^uWtcZ{7nq5k3@V;B@F)zyIRp(QnlK9}P}QB>(^b diff --git a/tests/regression_tests/deplete_with_reactivity_control/ref_depletion_with_rotation.h5 b/tests/regression_tests/deplete_with_reactivity_control/ref_depletion_with_rotation.h5 index 3831319390ece1aebcc4ee14e36bc11b913ad9e0..9c772c6187bdadafca834e7cb71637f3c207730c 100644 GIT binary patch delta 1427 zcmZ{jPiWIn9LJM%c1ir#Io#UXuAy$O+F84;j>2F~x8mv;t_sdFu`A+MJ-E6}Srya@ zqDNbQmklP=>0x>p9Z@I_GzAeaLh)ix8a&8OS|{RE5R<%@)lv*3?|r|&&-?w}@9~lo zDCq;4a!^t{z3&V}JZ3tjC)5Y>)aY_yjs%XSAxqJ00{Mj0na)Rm(?!8BO@qP13U_r& z6zKuK{`eM}2xF|!k=X0aqOC}J(`M1~Rjxdt@RI0oKL%n~2aE*5)hn{&@`@EQ2$RDOdIA*>LS3+xiN=v4U{WWy zEYwZEU<5Z^B0=mo!i8~;abUtA^EPTRF8rGM%%~A>MI@(0^ehv^E-mwo)a$q~a(f%i zgeu$ucEkU=Oc5(W@L#Q^Lths5zd7;89eD057`{EZ7x=xh_%M5aEs*#$zYt?fzQNA+ zm0dHFKHLBrjNf^6B919za&^;J-;%{Tdf9$hOy}LL|I4{{>g3EkS@QV^@~2apf;)4T Pt?J^7HzUoP&4c1!4^JP? literal 39352 zcmeHPeQX@X6`#GcorAd4CN8)L=$V$rDOTf;1dfl5FSe7MMb#K6 z(YZ}4ml6>~pl%WEMHCd&h#09tM?s0*R>BFWX;if&AVroS#Uivq^q(}PksHyZJNw?d zy&HS)d}|-E;I>#h{?@OVU#C^IkVbP>I_l|Hh7U}O3Wgde@W{w~7 zHDCtP!RX;|5ML-P6(}#RFv@#GSEMIAI6M$Wak0t}EH@$~-x(e-1PqJbXgD_1GZ;fb z&5)2X#5bz{4DziKc?5jZZ26V|kAQ8VT5z~r5FowoA>yAyU0DyBu&*h~vlus;|lm;>|!Q=p(rV#>cOMbRU`vZpd!; zQsvr>yVEP7Mx$6U@B_lu4_vDHv)!eEfa(W`M}8b1sFuVN;v=Oy4Z#n{(G`xV6%djg z%J3@naBA;zEl1E&@X6E}`!B0&9@ z`Ynxf3G{0uP2=5S`8609A!o8&wdkWqaNgq0{4xE3x* zZUEz>pNDiW$|bMMl$*l%%=IX+!~+^vRlKZq{MF{T8XD~H_4W<*ekUCD(ikh7W;f2_ zuxXU_n+RyeC@aV3XgJgX%l*(`I89S7J$8O@AKo|iS0>DS%7K8LzY1VKqT@ltn6m6 zP>c2~emLV9zE)lE^G4ER$Fs6qpx(?g>LDCvL%@z_<6p@vRVbdttHDkJ?JbmN9;~*d z$JuyBdhB>MPz#2bXPNzkV{;>5$1^X!N2gMY%`?*59IIJZ)gps zbwy(4J{#f{jh|^nvR-P!ov}ZL_M^?*0mX=??EkaExa>34{_8AnK>ouPEJV z2=h=^>ifMj$QQCxDJRdy;Zyu;;7&*X5$N&Z%)!!OQ7)@pgX2gZqQ%{V%BpE!h`a2;X-ox;{y*C+^!PO16Ebi@sX-&m@_W(oq2&$_1tTG&K!+|dSaN- zH`LJ+=?qge_1!Ziz^BwAmF|ri3XIHfEG-d-%IR>{rTE8xfM?#AqL)s(I7 z<5Bjnb)#IMnF`DuwqvrdY3}?>ojRLm1?yM7pF(5PI3*M44qj>`z`k5hXmk8;PQJ64 z?q7mFDo1`PS1U&onG#e^4QU?Z$rR@pBftnS0*nA7zz8q`i~u9R2rvSSz^9CW-8yog z!+MP4dXmc1dT_q#rt@4^(sM;#U(Od{mSY4M0Y-okU<4QeMt~7u1Q-EEfDvE>7LI`3 z^MtZ_st)4%pX75rKN#7YshH#`IGaEMt~7u1Q-EEfDvE>awA|jZ_M>wSRCh#BuCGURg3#<=8Ysz<<)W= zX9O4lMt~7u1Q-EEfDu?)2;@6&9LX_X9J84x^89dV-FL1HBftnS0*nA7kPia+&KGv) zm@g!GzK{?6EXxQm0*nA7zzAHH2vmK+KXL1Y+Lm{#{;E&B`fSVQEyL}5omPmckbY#>S=}=>`QFnj&oth6%oly=#9-#XSXTbWZ?Y+2j9 zf4J^%C)OWr5+8i@!{OJbnhqa2xKrXZQ~SV?=M9`{{&4Z$bex{JtMk5Nhb^4`(Ns5b z@Glllr}p?x)@-zJdf_Ku8`mqIP2up^iNA-RxzLfqY2wa*zWB5CFQoE6Zg}Ousoxye zf8KIe>(zhWrBC@vK6lM)lX~p>iMBm8r}cl`^W4D=#~#r~n{GRIcNH&Hl$oVZQ>`j-GreMO8C{5@tAB&q5OJPb{Nw?533`od7%0Ei774E5BbIZ}?WJS}6cYu5ui2!>(x8dKY%%-uad?INAc2_%Py$*#a-x zj%&>~Tg$jat5!6$3%k))t!TO(7d7uKgJnU3qkFKM5HvW@i(#ju3@)^T)o~qC4kvu{ z1|eo-uL>BTHpvF)+g=(81NtE|^H}SgNOXwI>c(NpO;)dSrD-rA;#WTvBrDg!(j=JB z&+8nSTjvr)%kA;(73C*kxsrqTz{P^#bQoL5%z*Hp;X?CgPhgyFd3PrG@V9KD3EY;4u zLXwj$(jX6-;o=031VWe;-p-^V(y!_7NQ=13A{jMevoD;bu9Wb^bT^+nx8+-REnb$NVpf3$hk{ki79J)rC=fiX_>i4@KA&2P?>}hKg!YN2Ug_tcym%vp36MA)Nd2|_P zKHcaDbhNc}$~IDe$zX%iLAt>^zv@e)7KwOjtoXoT&&gwSLzXxKNua<1$2x4M063dKxFF*5oHn185NH348#kg^-7y;_P z)Ng5=%b;H)X&Ue5%df$>T0_HusQ)ZjB0FJ>r?DPhtjHwxZ<2GbfsDGFA@reK$1QL{ za)THjy%o~^D3@DXEH{DiE4QJ*0uN|hRq-;{@vDq+HIj@E_+zPo1JQ(^##q@jt8o^G zO}%W;K)|7ovT}S*L?hj>+>a!q1)9p~vGRjw@V>FXGGOLY83e5SRR{JXDo(tgFMmZm zQ~iy}4|_bT#R{8x?D35BSn)#8~t-!=JRk7rb0Q;$8Kksd3aP542Ac}7Ta zSPB6vo*nfUm#P-ef_H$OF0^Oz!yeCW!fKm(?D35BSn({l9t<(h5O5AkB4EX{_di!$ ztXe$Ft^+%7pgofx_IO6GnN21cyL(?`rvh!fWHu!B2vapOaja1;b(D9ckl>Ak9yosR?Fn=TRP7K9Jim1Iv{XCe*7 z(RE;*(7psyk3Ig79^jAl{(_uz>iR+7jX;dgiv z=^HZk*$}U2{46L^^3ow(Mf+oDKiZWqLNVeg`~Q+KF8WNh|GE*8F#0tDi25n@D@wN; z!aUTO|9R1Qh?;u^rM$ZV=<68l*3_bAB9E@Ur_|yp>;ZV&-!U!h`b5;X-ni;{y*CT+Su#y=GMk@sX-&nA0ySoq2&$^*o?|&Yb9r^bcW1 zEY;oLcQ8uP{CCfk0H0Eee7ZlOE6_8eeTSm4;Yfc9izfTzt|TFb;;BS;bib4`ESUSI zK|+34y}vsiOYTp^<4F^DbzL(3{z3a4u9@#2-3Rr4QF`U~$8Xf%<*W(!7+3lJtK63^ z$gT2uIv&t}9L{1qsJ(;(#Vc<=+!2@4l}VWGFKBzZ$}KLfeztfJzX=bD&jF8jj+l5z zb;!|$h6mg~hw^_2goFL|JKf(m+h5?Hj+ee+pK>nD+iL~YD%W1R0zu5~Z7wZWN7?E= z9%cVpFUoaXK+YW*tkLuhOPr(^=%!Rw6#*q7@G#|;0QQ|T_S69i^Vy{2rvSS03*N%FanGKBftnS0*nA7@F^o;wT?Vv zvmRr(o}}`$9-OVZ$qLt%^jwkGm$OBf7y(9r5nu!u0Y-okU<4R}xg%ip zJi$9l)nQ!!lYF`72S-m9D<*jg`T0OGJ{)5N7y(9r5nu!u0Y-okU<4Kt0(RFAy|Y|5 zj9RQ0c%872dI8sk5nu!u0Y-okU<4R}@(5VX8_Rtcmd1G_$ed}M>{Xsi_?&;R{og-aOp1s^kjop;n z^~}D%oXh_!^wkZ!mTmdP&+y;Ih4@$3cVBqoKx00C_-)_v7qh4HanFzb(Z1*R-Ff+E z4;`Lr`RKWP{_U{=&m+C>YTM8Ja(MTx$F!&3cy#^p)fcr3BM07#wM}WMy| Date: Mon, 19 Jan 2026 14:35:00 +0100 Subject: [PATCH 2/4] remove obsolete bits and update unit test --- openmc/deplete/abc.py | 4 +- openmc/search.py | 39 +-- .../test_deplete_reactivity_control.py | 312 ++++++++++-------- 3 files changed, 177 insertions(+), 178 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index de76a891bc..47ee997c4a 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -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, \ diff --git a/openmc/search.py b/openmc/search.py index e372a40a0f..cff4ef14d5 100644 --- a/openmc/search.py +++ b/openmc/search.py @@ -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,11 +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 All remaining keyword arguments are passed to the root-finding method. @@ -202,26 +191,12 @@ 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 \ No newline at end of file diff --git a/tests/unit_tests/test_deplete_reactivity_control.py b/tests/unit_tests/test_deplete_reactivity_control.py index a202226db8..506ef42641 100644 --- a/tests/unit_tests/test_deplete_reactivity_control.py +++ b/tests/unit_tests/test_deplete_reactivity_control.py @@ -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() From 5737c39e415914c9da1c064d3c3d84d9113e9fa2 Mon Sep 17 00:00:00 2001 From: church89 Date: Mon, 19 Jan 2026 14:40:08 +0100 Subject: [PATCH 3/4] minor fix --- openmc/deplete/abc.py | 2 +- openmc/search.py | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 47ee997c4a..8d971107ba 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -34,7 +34,6 @@ from .transfer_rates import TransferRates, ExternalSourceRates from .reactivity_control import ReactivityController - __all__ = [ "OperatorResult", "TransportOperator", "ReactionRateHelper", "NormalizationHelper", "FissionYieldHelper", @@ -909,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: diff --git a/openmc/search.py b/openmc/search.py index cff4ef14d5..7566f57cb7 100644 --- a/openmc/search.py +++ b/openmc/search.py @@ -105,6 +105,8 @@ 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. + + .. versionadded:: 0.13.1 **kwargs All remaining keyword arguments are passed to the root-finding method. From 22b1fe29eb24eb33a8211d4861808225352ff89d Mon Sep 17 00:00:00 2001 From: church89 Date: Mon, 19 Jan 2026 14:41:43 +0100 Subject: [PATCH 4/4] minor fix --- openmc/search.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openmc/search.py b/openmc/search.py index 7566f57cb7..cd626e4f77 100644 --- a/openmc/search.py +++ b/openmc/search.py @@ -201,4 +201,5 @@ def search_for_keff(model_builder, initial_guess=None, target=1.0, # Perform the search zero_value = root_finder(**args) - return zero_value, guesses, results \ No newline at end of file + return zero_value, guesses, results + \ No newline at end of file