From 3e9733b4ff0ec00618fb65a38b1ea42e9f5668c5 Mon Sep 17 00:00:00 2001 From: church89 Date: Mon, 19 Jan 2026 11:51:03 +0100 Subject: [PATCH] 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 9108396f1..de76a891b 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 08e86d9a6..a70018375 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|