diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index 8a9509e9009..5e0063c90e5 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -17,6 +17,7 @@ from .results import * from .integrators import * from .transfer_rates import * +from .reactivity_control import * from . import abc from . import cram from . import helpers diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 5d9b8a2403b..8c6dda6babc 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -18,7 +18,7 @@ import numpy as np from uncertainties import ufloat -from openmc.checkvalue import check_type, check_greater_than, PathLike +from openmc.checkvalue import check_value, check_type, check_greater_than, PathLike from openmc.mpi import comm from openmc.utility_funcs import change_directory from openmc import Material @@ -28,7 +28,12 @@ from .pool import deplete from .reaction_rates import ReactionRates from .transfer_rates import TransferRates - +from openmc import Material, Cell +from .reactivity_control import ( + ReactivityController, + CellReactivityController, + MaterialReactivityController +) __all__ = [ "OperatorResult", "TransportOperator", @@ -543,6 +548,9 @@ class Integrator(ABC): transfer_rates : openmc.deplete.TransferRates Instance of TransferRates class to perform continuous transfer during depletion + reactivity_control : openmc.deplete.ReactivityController + Instance of ReactivityController class to perform reactivity control during + transport-depletion simulation. .. versionadded:: 0.14.0 @@ -633,6 +641,7 @@ def __init__( self.source_rates = np.asarray(source_rates) self.transfer_rates = None + self._reactivity_control = None if isinstance(solver, str): # Delay importing of cram module, which requires this file @@ -678,6 +687,15 @@ def solver(self, func): self._solver = func + @property + def reactivity_control(self): + return self._reactivity_control + + @reactivity_control.setter + def reactivity_control(self, reactivity_control): + check_type('reactivity control', reactivity_control, ReactivityController) + self._reactivity_control = reactivity_control + def _timed_deplete(self, n, rates, dt, matrix_func=None): start = time.time() results = deplete( @@ -764,6 +782,13 @@ def _get_start_data(self): return (self.operator.prev_res[-1].time[-1], len(self.operator.prev_res) - 1) + def _get_bos_from_reactivity_control(self, step_index, bos_conc): + """Get BOS from reactivity control.""" + x = deepcopy(bos_conc) + # Get new vector after keff criticality control + x, root = self._reactivity_control.search_for_keff(x, step_index) + return x, root + def integrate( self, final_step: bool = True, @@ -798,10 +823,18 @@ def integrate( # Solve transport equation (or obtain result from restart) if i > 0 or self.operator.prev_res is None: + # Update geometry/material according to reactivity control + if self._reactivity_control is not None and source_rate != 0.0: + n, root = self._get_bos_from_reactivity_control(i, n) + else: + root = None n, res = self._get_bos_data_from_operator(i, source_rate, n) else: n, res = self._get_bos_data_from_restart(source_rate, n) - + if self._reactivity_control: + root = self.operator.prev_res[-1].reac_cont + else: + root = None # Solve Bateman equations over time interval proc_time, n_list, res_list = self(n, res.rates, dt, source_rate, i) @@ -811,9 +844,8 @@ def integrate( # Remove actual EOS concentration for next step n = n_list.pop() - StepResult.save(self.operator, n_list, res_list, [t, t + dt], - source_rate, self._i_res + i, proc_time, path) + source_rate, self._i_res + i, proc_time, root, path) t += dt @@ -823,9 +855,13 @@ def integrate( # solve) if output and final_step and comm.rank == 0: print(f"[openmc.deplete] t={t} (final operator evaluation)") + if self._reactivity_control is not None and source_rate != 0.0: + n, root = self._get_bos_from_reactivity_control(i+1, n) + else: + root = None res_list = [self.operator(n, source_rate if final_step else 0.0)] StepResult.save(self.operator, [n], res_list, [t, t], - source_rate, self._i_res + len(self), proc_time, path) + source_rate, self._i_res + len(self), proc_time, root, path) self.operator.write_bos_data(len(self) + self._i_res) self.operator.finalize() @@ -864,6 +900,29 @@ def add_transfer_rate( self.transfer_rates.set_transfer_rate(material, components, transfer_rate, transfer_rate_units, destination_material) + def add_reactivity_control( + self, + obj: Union[Cell, Material], + **kwargs): + """Add pre-defined Cell based or Material bases reactivity control to + integrator scheme. + + Parameters + ---------- + obj : openmc.Cell or openmc.Material + Cell or Material identifier to where add reactivity control + **kwargs + keyword arguments that are passed to the specific ReactivityController + class. + + """ + 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) + @add_params class SIIntegrator(Integrator): r"""Abstract class for the Stochastic Implicit Euler integrators @@ -1023,13 +1082,13 @@ def integrate( n = n_list.pop() StepResult.save(self.operator, n_list, res_list, [t, t + dt], - p, self._i_res + i, proc_time, path) + p, self._i_res + i, proc_time, path=path) t += dt # No final simulation for SIE, use last iteration results StepResult.save(self.operator, [n], [res_list[-1]], [t, t], - p, self._i_res + len(self), proc_time, path) + p, self._i_res + len(self), proc_time, path=path) self.operator.write_bos_data(self._i_res + len(self)) self.operator.finalize() diff --git a/openmc/deplete/reactivity_control.py b/openmc/deplete/reactivity_control.py new file mode 100644 index 00000000000..08e86d9a6ce --- /dev/null +++ b/openmc/deplete/reactivity_control.py @@ -0,0 +1,1031 @@ +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 + + 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 + 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 + + """ + + def __init__( + self, + operator, + bracket, + bracket_limit, + bracketed_method="brentq", + tol=0.01, + target=1.0, + print_iterations=True, + search_for_keff_output=True, + ): + + 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 + + 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. + + 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 + """ + + @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 + Returns + ------- + volumes : dict + Dictionary of calculated volume, where key mat id and value + material volume, in cm3 + """ + + 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, + ) + + # 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") + + return root + + def _update_materials(self, x): + """Update number density and material compositions in OpenMC on all + processes. + + Parameters + ---------- + x : list of numpy.ndarray + Total atom concentrations + """ + 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 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 diff --git a/openmc/deplete/stepresult.py b/openmc/deplete/stepresult.py index 9cf33898f3b..bc112f5f2b2 100644 --- a/openmc/deplete/stepresult.py +++ b/openmc/deplete/stepresult.py @@ -56,6 +56,8 @@ class StepResult: Number of stages in simulation. data : numpy.ndarray Atom quantity, stored by stage, mat, then by nuclide. + reac_cont : float + The root returned by the reactivity controller. proc_time : int Average time spent depleting a material across all materials and processes @@ -74,6 +76,7 @@ def __init__(self): self.mat_to_hdf5_ind = None self.data = None + self.reac_cont = None def __repr__(self): t = self.time[0] @@ -349,6 +352,10 @@ def _write_hdf5_metadata(self, handle): "depletion time", (1,), maxshape=(None,), dtype="float64") + handle.create_dataset( + "reac_cont_root", (1,), maxshape=(None,), + dtype="float64") + def _to_hdf5(self, handle, index, parallel=False): """Converts results object into an hdf5 object. @@ -379,6 +386,7 @@ def _to_hdf5(self, handle, index, parallel=False): time_dset = handle["/time"] source_rate_dset = handle["/source_rate"] proc_time_dset = handle["/depletion time"] + root_dset = handle["/reac_cont_root"] # Get number of results stored number_shape = list(number_dset.shape) @@ -412,6 +420,10 @@ def _to_hdf5(self, handle, index, parallel=False): proc_shape[0] = new_shape proc_time_dset.resize(proc_shape) + root_shape = list(root_dset.shape) + root_shape[0] = new_shape + root_dset.resize(root_shape) + # If nothing to write, just return if len(self.index_mat) == 0: return @@ -435,6 +447,7 @@ def _to_hdf5(self, handle, index, parallel=False): proc_time_dset[index] = ( self.proc_time / (comm.size * self.n_hdf5_mats) ) + root_dset[index] = self.reac_cont @classmethod def from_hdf5(cls, handle, step): @@ -469,6 +482,10 @@ def from_hdf5(cls, handle, step): if step < proc_time_dset.shape[0]: results.proc_time = proc_time_dset[step] + if "reac_cont_root" in handle: + root_dset = handle["/reac_cont_root"] + results.reac_cont = root_dset[step] + if results.proc_time is None: results.proc_time = np.array([np.nan]) @@ -509,7 +526,7 @@ def from_hdf5(cls, handle, step): @staticmethod def save(op, x, op_results, t, source_rate, step_ind, proc_time=None, - path: PathLike = "depletion_results.h5"): + root=None, path: PathLike = "depletion_results.h5"): """Creates and writes depletion results to disk Parameters @@ -530,7 +547,8 @@ def save(op, x, op_results, t, source_rate, step_ind, proc_time=None, Total process time spent depleting materials. This may be process-dependent and will be reduced across MPI processes. - + root : float + The root returned by the reactivity controller. path : PathLike Path to file to write. Defaults to 'depletion_results.h5'. @@ -564,6 +582,7 @@ def save(op, x, op_results, t, source_rate, step_ind, proc_time=None, results.proc_time = proc_time if results.proc_time is not None: results.proc_time = comm.reduce(proc_time, op=MPI.SUM) + results.reac_cont = root if not Path(path).is_file(): Path(path).parent.mkdir(parents=True, exist_ok=True) diff --git a/openmc/search.py b/openmc/search.py index 7f6b3b5b9a4..9b18e870025 100644 --- a/openmc/search.py +++ b/openmc/search.py @@ -1,5 +1,6 @@ from collections.abc import Callable from numbers import Real +from warnings import warn import scipy.optimize as sopt @@ -12,7 +13,7 @@ def _search_keff(guess, target, model_builder, model_args, print_iterations, - run_args, guesses, results): + run_args, guesses, results, run_in_memory): """Function which will actually create our model, run the calculation, and obtain the result. This function will be passed to the root finding algorithm @@ -39,7 +40,8 @@ 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 @@ -51,7 +53,11 @@ def _search_keff(guess, target, model_builder, model_args, print_iterations, model = model_builder(guess, **model_args) # Run the model and obtain keff - sp_filepath = model.run(**run_args) + if run_in_memory: + openmc.lib.run(**run_args) + sp_filepath = f'statepoint.{model.settings.batches}.h5' + else: + sp_filepath = model.run(**run_args) with openmc.StatePoint(sp_filepath) as sp: keff = sp.keff @@ -70,7 +76,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, **kwargs): + run_args=None, run_in_memory=False, **kwargs): """Function to perform a keff search by modifying a model parametrized by a single independent variable. @@ -105,6 +111,9 @@ 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 @@ -193,12 +202,26 @@ 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_args, guesses, results, run_in_memory) # Create a new dictionary with the arguments from args and kwargs args.update(kwargs) # Perform the search - zero_value = root_finder(**args) + if not run_in_memory: + zero_value = root_finder(**args) + return zero_value, guesses, results - return zero_value, guesses, results + 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 diff --git a/tests/regression_tests/deplete_with_reactivity_control/__init__.py b/tests/regression_tests/deplete_with_reactivity_control/__init__.py new file mode 100644 index 00000000000..e69de29bb2d 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 new file mode 100644 index 00000000000..83a19702b0e Binary files /dev/null and b/tests/regression_tests/deplete_with_reactivity_control/ref_depletion_with_refuel.h5 differ 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 new file mode 100644 index 00000000000..3831319390e Binary files /dev/null and b/tests/regression_tests/deplete_with_reactivity_control/ref_depletion_with_rotation.h5 differ diff --git a/tests/regression_tests/deplete_with_reactivity_control/ref_depletion_with_translation.h5 b/tests/regression_tests/deplete_with_reactivity_control/ref_depletion_with_translation.h5 new file mode 100644 index 00000000000..6ec626a2ffe Binary files /dev/null and b/tests/regression_tests/deplete_with_reactivity_control/ref_depletion_with_translation.h5 differ diff --git a/tests/regression_tests/deplete_with_reactivity_control/test.py b/tests/regression_tests/deplete_with_reactivity_control/test.py new file mode 100644 index 00000000000..c9d0023f209 --- /dev/null +++ b/tests/regression_tests/deplete_with_reactivity_control/test.py @@ -0,0 +1,140 @@ +""" Tests for ReactivityController class """ + +from pathlib import Path +import shutil +import sys + +import pytest +import numpy as np + +import openmc +import openmc.lib +from openmc.deplete import CoupledOperator +from openmc.deplete import ( + CellReactivityController, + MaterialReactivityController +) + +from tests.regression_tests import config + +@pytest.fixture +def model(): + f = openmc.Material(name='f') + f.set_density('g/cm3', 10.29769) + f.add_element('U', 1., enrichment=2.4) + f.add_element('O', 2.) + + h = openmc.Material(name='h') + h.set_density('g/cm3', 0.001598) + h.add_element('He', 2.4044e-4) + + w = openmc.Material(name='w') + w.set_density('g/cm3', 0.740582) + w.add_element('H', 2) + w.add_element('O', 1) + + # Define overall material + materials = openmc.Materials([f, h, w]) + + # Define surfaces + radii = [0.5, 0.8, 1] + height = 80 + surf_in = openmc.ZCylinder(r=radii[0]) + surf_mid = openmc.ZCylinder(r=radii[1]) + surf_out = openmc.ZCylinder(r=radii[2], boundary_type='reflective') + surf_top = openmc.ZPlane(z0=height/2, boundary_type='vacuum') + surf_bot = openmc.ZPlane(z0=-height/2, boundary_type='vacuum') + + surf_trans = openmc.ZPlane(z0=0) + surf_rot1 = openmc.XPlane(x0=0) + surf_rot2 = openmc.YPlane(y0=0) + + # Define cells + cell_f = openmc.Cell(name='fuel_cell', fill=f, + region=-surf_in & -surf_top & +surf_bot) + cell_g = openmc.Cell(fill=h, + region = +surf_in & -surf_mid & -surf_top & +surf_bot & +surf_rot2) + + # Define unbounded cells for rotation universe + cell_w = openmc.Cell(fill=w, region = -surf_rot1) + cell_h = openmc.Cell(fill=h, region = +surf_rot1) + universe_rot = openmc.Universe(cells=(cell_w, cell_h)) + cell_rot = openmc.Cell(name="rot_cell", fill=universe_rot, + region = +surf_in & -surf_mid & -surf_top & +surf_bot & -surf_rot2) + + # Define unbounded cells for translation universe + cell_w = openmc.Cell(fill=w, region=+surf_in & -surf_trans ) + cell_h = openmc.Cell(fill=h, region=+surf_in & +surf_trans) + universe_trans = openmc.Universe(cells=(cell_w, cell_h)) + cell_trans = openmc.Cell(name="trans_cell", fill=universe_trans, + region=+surf_mid & -surf_out & -surf_top & +surf_bot) + + # Define overall geometry + geometry = openmc.Geometry([cell_f, cell_g, cell_rot, cell_trans]) + + # Set material volume for depletion fuel. + f.volume = np.pi * radii[0]**2 * height + + settings = openmc.Settings() + settings.particles = 1000 + settings.inactive = 10 + settings.batches = 50 + + return openmc.Model(geometry, materials, settings) + +@pytest.fixture +def mix_mat(): + mix_mat = openmc.Material() + mix_mat.add_element("U", 1, percent_type="ao", enrichment=90) + mix_mat.add_element("O", 2) + mix_mat.set_density("g/cc", 10.4) + return mix_mat + +@pytest.mark.parametrize("obj, attribute, bracket, bracket_limit, axis, ref_result", [ + ('trans_cell', 'translation', [-5,5], [-40,40], 2, 'depletion_with_translation'), + ('rot_cell', 'rotation', [-5,5], [-90,90], 2, 'depletion_with_rotation'), + ('f', None, [-1,2], [-100,100], None, 'depletion_with_refuel') + ]) +def test_reactivity_control(run_in_tmpdir, model, obj, attribute, bracket, + bracket_limit, axis, ref_result, mix_mat): + + chain_file = Path(__file__).parents[2] / 'chain_simple.xml' + op = CoupledOperator(model, chain_file) + + integrator = openmc.deplete.PredictorIntegrator( + op, [1], 174., timestep_units = 'd') + + if attribute: + integrator.reactivity_control = CellReactivityController( + cell=obj, + operator=op, + attribute=attribute, + bracket=bracket, + bracket_limit=bracket_limit, + axis=axis + ) + else: + integrator.reactivity_control = MaterialReactivityController( + material=obj, + operator=op, + material_to_mix=mix_mat, + bracket=bracket, + bracket_limit=bracket_limit, + ) + integrator.integrate() + + # Get path to test and reference results + path_test = op.output_dir / 'depletion_results.h5' + path_reference = Path(__file__).with_name(f'ref_{ref_result}.h5') + + # If updating results, do so and return + if config['update']: + shutil.copyfile(str(path_test), str(path_reference)) + return + + # Load the reference/test results + res_test = openmc.deplete.Results(path_test) + res_ref = openmc.deplete.Results(path_reference) + + # Use high tolerance here + assert res_test[0].reac_cont == pytest.approx(res_ref[0].reac_cont, rel=2) diff --git a/tests/unit_tests/test_deplete_reactivity_control.py b/tests/unit_tests/test_deplete_reactivity_control.py new file mode 100644 index 00000000000..a202226db89 --- /dev/null +++ b/tests/unit_tests/test_deplete_reactivity_control.py @@ -0,0 +1,222 @@ +""" Tests for ReactivityController class """ + +from pathlib import Path + +import pytest +import numpy as np + +import openmc +import openmc.lib +from openmc.deplete import CoupledOperator +from openmc.deplete import ( + CellReactivityController, + MaterialReactivityController +) + +CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml" + +@pytest.fixture +def model(): + f = openmc.Material(name="fuel") + f.add_element("U", 1, percent_type="ao", enrichment=4.25) + f.add_element("O", 2) + f.set_density("g/cc", 10.4) + f.temperature = 293.15 + + w = openmc.Material(name="water") + w.add_element("O", 1) + w.add_element("H", 2) + w.set_density("g/cc", 1.0) + w.temperature = 293.15 + w.depletable = True + + h = openmc.Material(name='helium') + h.add_element('He', 1) + h.set_density('g/cm3', 0.001598) + + radii = [0.42, 0.45] + height = 0.5 + + f.volume = np.pi * radii[0] ** 2 * height + w.volume = np.pi * (radii[1]**2 - radii[0]**2) * height/2 + + materials = openmc.Materials([f, w, h]) + + surf_interface = openmc.ZPlane(z0=0) + surf_top = openmc.ZPlane(z0=height/2) + surf_bot = openmc.ZPlane(z0=-height/2) + surf_in = openmc.Sphere(r=radii[0]) + surf_out = openmc.Sphere(r=radii[1], boundary_type='vacuum') + + cell_water = openmc.Cell(fill=w, region=-surf_interface) + cell_helium = openmc.Cell(fill=h, region=+surf_interface) + universe = openmc.Universe(cells=(cell_water, cell_helium)) + cell_fuel = openmc.Cell(name='fuel_cell', fill=f, + region=-surf_in & -surf_top & +surf_bot) + cell_universe = openmc.Cell(name='universe_cell',fill=universe, + region=+surf_in & -surf_out & -surf_top & +surf_bot) + geometry = openmc.Geometry([cell_fuel, cell_universe]) + + settings = openmc.Settings() + settings.particles = 1000 + settings.inactive = 10 + settings.batches = 50 + + return openmc.Model(geometry, materials, settings) + +@pytest.fixture +def operator(model): + return CoupledOperator(model, CHAIN_PATH) + +@pytest.fixture +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 + ) + + model.export_to_xml() + openmc.lib.init() + + # 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)) + + # 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) + + 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()