From 7583f78750624945a8294eac5517256bc048f9b0 Mon Sep 17 00:00:00 2001 From: church89 Date: Wed, 5 Jun 2024 09:27:29 +0200 Subject: [PATCH 1/2] implement suggested changes to reactivity control classes --- openmc/deplete/abc.py | 60 +-- openmc/deplete/reactivity_control.py | 764 +++++++-------------------- 2 files changed, 213 insertions(+), 611 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index a1ab198c29..b8d2e1ef91 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -31,11 +31,9 @@ from .reaction_rates import ReactionRates from .transfer_rates import TransferRates from openmc import Material, Cell from .reactivity_control import ( + ReactivityController, CellReactivityController, - GeometricalCellReactivityController, - TemperatureCellReactivityController, - MaterialReactivityController, - RefuelMaterialReactivityController + MaterialReactivityController ) __all__ = [ @@ -718,25 +716,7 @@ class Integrator(ABC): @reactivity_control.setter def reactivity_control(self, reactivity_control): - if isinstance(reactivity_control, CellReactivityController): - if not isinstance(reactivity_control.cell, Cell): - raise TypeError( - "Reactivity control attribute must be Cell " - "not {}".format(type(reactivity_control.cell))) - elif isinstance(reactivity_control, MaterialReactivityController): - if not isinstance(reactivity_control.material, Material): - raise TypeError( - "Reactivity control attribute must be Material " - "not {}".format(type(reactivity_control.material))) - else: - raise TypeError( - "Reactivity Control must either be " - "CellReactivityController, or MaterialReactivityController, " - "not {}".format(type(reactivity_control))) - - check_value('attribute', reactivity_control.attrib_name, - ('translation', 'rotation', 'temperature', 'refuel')) - + check_type('reactivity control', reactivity_control, ReactivityController) self._reactivity_control = reactivity_control def _timed_deplete(self, n, rates, dt, matrix_func=None): @@ -945,34 +925,26 @@ class Integrator(ABC): def add_reactivity_control( self, - obj: Union[Cell, Material, int, str], - attr: str, + obj: Union[Cell, Material], **kwargs): - """Add reactivity control to integrator scheme. + """Add pre-defined Cell based or Material bases reactivity control to + integrator scheme. Parameters ---------- - obj : openmc.Cell or openmc.Material object or id or str name - Cell or Materials identifier to where add reactivity control - attr : str - Attribute to specify the type of reactivity control. Accepted values - are: 'translation', 'rotation', 'temperature' for an openmc.Cell - object; 'refuel' for an openmc.Material object. + obj : openmc.Cell or openmc.Material + Cell or Material identifier to where add reactivity control **kwargs - keyword arguments that are passed to the ReactivityController class. + keyword arguments that are passed to the specific ReactivityController + class. """ - check_value('attribute', attr, ('translation', 'rotation', - 'temperature', 'refuel')) - if attr in ('translation', 'rotation'): - reactivity_control = GeometricalCellReactivityController - elif attr == 'temperature': - reactivity_control = TemperatureCellReactivityController - elif attr == 'refuel': - reactivity_control = RefuelMaterialReactivityController - - self._reactivity_control = reactivity_control.from_params(obj, attr, - self.operator, **kwargs) + 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): diff --git a/openmc/deplete/reactivity_control.py b/openmc/deplete/reactivity_control.py index e36f64b79e..08e86d9a6c 100644 --- a/openmc/deplete/reactivity_control.py +++ b/openmc/deplete/reactivity_control.py @@ -14,6 +14,7 @@ from openmc.checkvalue import ( check_type, check_value, check_less_than, + check_greater_than, check_iterable_type, check_length, ) @@ -33,11 +34,9 @@ class ReactivityController(ABC): This abstract class sets the requirements for the reactivity control set. Users should instantiate - :class:`openmc.deplete.reactivity_control.GeometricalCellReactivityController` + :class:`openmc.deplete.reactivity_control.CellReactivityController` or - :class:`openmc.deplete.reactivity_control.TemperatureCellReactivityController` - or - :class:`openmc.deplete.reactivity_control.RefuelMaterialReactivityController` + :class:`openmc.deplete.reactivity_control.MaterialReactivityController` rather than this class. .. versionadded:: 0.14.1 @@ -52,10 +51,6 @@ class ReactivityController(ABC): 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. - density_treatment : str, optional - Whether or not to keep constant volume or density after a depletion step - before the next one. - Default to 'constant-volume' bracketed_method : {'brentq', 'brenth', 'ridder', 'bisect'}, optional Solution method to use. This is equivalent to the `bracket_method` parameter of the @@ -90,7 +85,6 @@ class ReactivityController(ABC): operator, bracket, bracket_limit, - density_treatment="constant-volume", bracketed_method="brentq", tol=0.01, target=1.0, @@ -103,18 +97,12 @@ class ReactivityController(ABC): self.local_mats = operator.local_mats self.model = operator.model self.geometry = operator.model.geometry - - check_value( - "density_treatment", - density_treatment, - ("constant-density", "constant-volume"), - ) - self.density_treatment = density_treatment 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 @@ -130,7 +118,8 @@ class ReactivityController(ABC): def bracketed_method(self, value): check_value("bracketed_method", value, _SCALAR_BRACKETED_METHODS) if value != "brentq": - warn("brentq bracketed method is recommended") + warn("""brentq bracketed method is recommended to ensure + convergence of the method""") self._bracketed_method = value @property @@ -152,8 +141,8 @@ class ReactivityController(ABC): self._target = value @classmethod - def from_params(cls, obj, attr, operator, **kwargs): - return cls(obj, attr, operator, **kwargs) + def from_params(cls, obj, operator, **kwargs): + return cls(obj, operator, **kwargs) @abstractmethod def _model_builder(self, param): @@ -193,17 +182,36 @@ class ReactivityController(ABC): 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 - adapt it until :meth:`openmc.search.search_for_keff` return a valid + 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 to move the bracket towards the target. - A bracket limit pose the upper and lower boundaries to the adapting - bracket. If one limit is hit, the algorithm will stop and the closest - limit value will be used. + 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 ---------- @@ -341,36 +349,9 @@ class ReactivityController(ABC): return root - def _update_volumes(self): - """Update volumes stored in AtomNumber. - - After a depletion step, both material volume and density change, due to - changes in nuclides composition. - At present we lack an implementation to calculate density and volume - changes due to the different molecules speciation. Therefore, OpenMC - assumes by default that the depletable volume does not change and only - updates the nuclide densities and consequently the total material density. - This method, vice-versa, assumes that the total material density does not - change and update the material volumes instead. - """ - number_i = self.operator.number - for mat_idx, mat in enumerate(self.local_mats): - # Total number of atoms-gram per mol - agpm = 0 - for nuc in number_i.nuclides: - agpm += number_i[mat, nuc] * atomic_mass(nuc) - # Get mass dens from beginning, intended to be held constant - density = openmc.lib.materials[int(mat)].get_density("g/cm3") - number_i.volume[mat_idx] = agpm / AVOGADRO / density - def _update_materials(self, x): - """Update number density and material compositions in OpenMC on all processes. - - If density_treatment is set to 'constant-density' - :meth:`openmc.deplete.reactivity_control._update_volumes` is called to update - material volumes in AtomNumber, keeping the material total density - constant, before re-normalizing the atom densities and assigning them - to the model in memory. + """Update number density and material compositions in OpenMC on all + processes. Parameters ---------- @@ -379,9 +360,6 @@ class ReactivityController(ABC): """ self.operator.number.set_density(x) - if self.density_treatment == "constant-density": - self._update_volumes() - for rank in range(comm.size): number_i = comm.bcast(self.operator.number, root=rank) @@ -448,26 +426,19 @@ class ReactivityController(ABC): number_i.volume[mat_idx] = res_vol return x - class CellReactivityController(ReactivityController): """Abstract class holding reactivity control cell-based functions. - Specific classes for running reactivity control depletion calculations are - implemented as derived class of CellReactivityController. - Users should instantiate - :class:`openmc.deplete.reactivity_control.GeometricalCellReactivityController` - or - :class:`openmc.deplete.reactivity_control.TemperatureCellReactivityController` - rather than this class. - .. versionadded:: 0.14.1 Parameters ---------- cell : openmc.Cell or int or str - OpenMC Cell identifier to where apply batch wise scheme + 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. @@ -475,10 +446,6 @@ class CellReactivityController(ReactivityController): 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. - density_treatment : str - Whether or not to keep constant volume or density after a depletion step - before the next one. - Default to 'constant-volume' bracketed_method : {'brentq', 'brenth', 'ridder', 'bisect'}, optional Solution method to use. This is equivalent to the `bracket_method` parameter of the @@ -502,7 +469,9 @@ class CellReactivityController(ReactivityController): ---------- cell : openmc.Cell or int or str OpenMC Cell identifier to where apply batch wise scheme - universe_cells : list of openmc.Cell + 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. @@ -518,13 +487,14 @@ class CellReactivityController(ReactivityController): self, cell, operator, + attribute, bracket, bracket_limit, - axis=None, - density_treatment="constant-volume", + axis, bracketed_method="brentq", tol=0.01, target=1.0, + samples=1000000, print_iterations=True, search_for_keff_output=True, ): @@ -533,7 +503,6 @@ class CellReactivityController(ReactivityController): operator, bracket, bracket_limit, - density_treatment, bracketed_method, tol, target, @@ -543,19 +512,41 @@ class CellReactivityController(ReactivityController): 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 - if axis is not None: - # index of cell directional axis - check_value("axis", axis, (0, 1, 2)) + # 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 - # Initialize container of universe cells, to populate only if materials - # are set as depletables - self.universe_cells = None + check_type("samples", samples, int) + self.samples = samples + + if self.depletable_cells: + self._initialize_volume_calc() - @abstractmethod def _get_cell_attrib(self): """Get cell attribute coefficient. @@ -564,8 +555,8 @@ class CellReactivityController(ReactivityController): coeff : float cell coefficient """ + return getattr(self.lib_cell, self.attribute)[self.axis] - @abstractmethod def _set_cell_attrib(self, val): """Set cell attribute to the cell instance. @@ -574,9 +565,13 @@ class CellReactivityController(ReactivityController): 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""" + """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] @@ -644,6 +639,33 @@ class CellReactivityController(ReactivityController): 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. @@ -703,347 +725,28 @@ class CellReactivityController(ReactivityController): # set results value as attribute in the geometry self._set_cell_attrib(root) - # if at least one of the cell materials is depletable, calculate new - # volume and update x and number accordingly + # 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.universe_cells: - volumes = self._calculate_volumes() + if self.depletable_cells: + volumes = self._adjust_volumes() x = self._update_x_and_set_volumes(x, volumes) return x, root - -class GeometricalCellReactivityController(CellReactivityController): - """Reactivity control cell-based with geometrical-attribute class. - - A user doesn't need to call this class directly. - Instead an instance of this class is automatically created by calling - :meth:`openmc.deplete.Integrator.add_reactivity_control` method from an - integrator class, such as :class:`openmc.deplete.CECMIntegrator`. - - .. versionadded:: 0.14.1 - - Parameters - ---------- - cell : openmc.Cell or int or str - OpenMC Cell identifier to where apply batch wise scheme - attrib_name : str - Cell attribute type - operator : openmc.deplete.Operator - OpenMC operator object - axis : int {0,1,2} - Directional axis for geometrical parametrization, where 0, 1 and 2 stand - for 'x', 'y' and 'z', respectively. - 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. - density_treatment : str - Whether or not to keep constant volume or density after a depletion step - before the next one. - Default to 'constant-volume' - 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 - attrib_name : str {'translation', 'rotation'} - Cell attribute type - axis : int {0,1,2} - Directional axis for geometrical parametrization, where 0, 1 and 2 stand - for 'x', 'y' and 'z', respectively. - samples : int - Number of samples used to generate volume estimates for stochastic - volume calculations. - - """ - - def __init__( - self, - cell, - attrib_name, - operator, - bracket, - bracket_limit, - axis, - density_treatment="constant-volume", - bracketed_method="brentq", - tol=0.01, - target=1.0, - samples=1000000, - print_iterations=True, - search_for_keff_output=True, - ): - - super().__init__( - cell, - operator, - bracket, - bracket_limit, - axis, - density_treatment, - bracketed_method, - tol, - target, - print_iterations, - search_for_keff_output, - ) - - check_value("attrib_name", attrib_name, ("rotation", "translation")) - self.attrib_name = attrib_name - - if isinstance(self.cell.fill, openmc.Universe): - # check if universe contains 2 cells - check_length("universe cells", self.cell.fill.cells, 2) - self.universe_cells = [ - cell for cell in self.cell.fill.cells.values() if cell.fill.depletable - ] - elif isinstance(self.cell.fill, openmc.Lattice): - # check if lattice contains 2 universes - check_length("lattice universes", self.cell.fill.get_unique_universes(), 2) - self.universe_cells = [ - cell - for cell in self.cell.fill.get_all_cells().values() - if cell.fill.depletable - ] - else: - raise ValueError( - "{} s not a valid instance of " - " should be a {} or {} instance".format( - self.cell.fill, openmc.Universe, openmc.Lattice - ) - ) - - check_type("samples", samples, int) - self.samples = samples - - if self.universe_cells: - self._initialize_volume_calc() - - def _get_cell_attrib(self): - """Get cell attribute coefficient. - - Returns - ------- - coeff : float - cell coefficient - """ - if self.attrib_name == "translation": - return self.lib_cell.translation[self.axis] - elif self.attrib_name == "rotation": - return self.lib_cell.rotation[self.axis] - - def _set_cell_attrib(self, val): - """Set cell attribute to the cell instance. - - Attributes are only applied to a cell filled with a universe containing - two cells itself. - - Parameters - ---------- - val : float - Cell coefficient to set, in cm for translation and deg for rotation - """ - if self.attrib_name == "translation": - vector = self.lib_cell.translation - elif self.attrib_name == "rotation": - vector = self.lib_cell.rotation - - vector[self.axis] = val - setattr(self.lib_cell, self.attrib_name, vector) - - 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.universe_cells, self.samples, ll, ur) - self.model.settings.volume_calculations = mat_vol - - def _calculate_volumes(self): - """Perform stochastic volume calculation - - 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.universe_cells: - mat_id = cell.fill.id - volumes[str(mat_id)] = res.volumes[cell.id].n - volumes = comm.bcast(volumes) - return volumes - - -class TemperatureCellReactivityController(CellReactivityController): - """Reactivity control cell-based with temperature-attribute class. - - A user doesn't need to call this class directly. - Instead an instance of this class is automatically created by calling - :meth:`openmc.deplete.Integrator.add_reactivity_control` method from an - integrator class, such as :class:`openmc.deplete.CECMIntegrator`. - - .. versionadded:: 0.14.1 - - Parameters - ---------- - cell : openmc.Cell or int or str - OpenMC Cell identifier to where apply batch wise scheme - operator : openmc.deplete.Operator - OpenMC operator object - attrib_name : str - Cell attribute type - 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. - density_treatment : str - Whether or not to keep constant volume or density after a depletion step - before the next one. - Default to 'constant-volume' - 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 - attrib_name : str {'temperature'} - Cell attribute type - - """ - - def __init__( - self, - cell, - attrib_name, - operator, - bracket, - bracket_limit, - axis=None, - density_treatment="constant-volume", - bracketed_method="brentq", - tol=0.01, - target=1.0, - print_iterations=True, - search_for_keff_output=True, - ): - - super().__init__( - cell, - operator, - bracket, - bracket_limit, - axis, - density_treatment, - bracketed_method, - tol, - target, - print_iterations, - search_for_keff_output, - ) - - # Not needed but used for consistency with other classes - check_value("attrib_name", attrib_name, "temperature") - self.attrib_name = attrib_name - - # check if initial temperature has been set to right cell material - if isinstance(self.cell.fill, openmc.Universe): - cells = [ - cell for cell in self.cell.fill.cells.values() if cell.fill.temperature - ] - check_length("Only one cell with temperature", cells, 1) - self.cell = cells[0] - - check_type("temperature cell real", self.cell.fill.temperature, Real) - - def _get_cell_attrib(self): - """Get cell temperature. - - Returns - ------- - coeff : float - cell temperature, in Kelvin - """ - return self.lib_cell.get_temperature() - - def _set_cell_attrib(self, val): - """Set temperature value to the cell instance. - - Parameters - ---------- - val : float - Cell temperature to set, in Kelvin - """ - self.lib_cell.set_temperature(val) - - class MaterialReactivityController(ReactivityController): """Abstract class holding reactivity control material-based functions. - Specific classes for running reactivity control depletion calculations are - implemented as derived class of MaterialReactivityController. - Users should instantiate - :class:`openmc.deplete.reactivity_control.RefuelMaterialReactivityController` - rather than this class. - .. versionadded:: 0.14.1 Parameters ---------- material : openmc.Material or int or str - OpenMC Material identifier to where apply batch wise scheme + OpenMC Material identifier to where apply reactivity control operator : openmc.deplete.Operator OpenMC operator object - mat_vector : dict - Dictionary of material composition to parameterize, where a pair key value - represents a nuclide and its weight fraction, respectively. + 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. @@ -1051,10 +754,12 @@ class MaterialReactivityController(ReactivityController): 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. - density_treatment : str - Whether or not to keep constant volume or density after a depletion step - before the next one. - Default to 'constant-volume' + 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 @@ -1078,9 +783,14 @@ class MaterialReactivityController(ReactivityController): ---------- material : openmc.Material or int or str OpenMC Material identifier to where apply batch wise scheme - mat_vector : dict - Dictionary of material composition to parameterize, where a pair key value - represents a nuclide and its weight fraction, respectively. + 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 """ @@ -1088,10 +798,11 @@ class MaterialReactivityController(ReactivityController): self, material, operator, - mat_vector, + material_to_mix, bracket, bracket_limit, - density_treatment="constant-volume", + units='cc', + dilute=False, bracketed_method="brentq", tol=0.01, target=1.0, @@ -1103,7 +814,6 @@ class MaterialReactivityController(ReactivityController): operator, bracket, bracket_limit, - density_treatment, bracketed_method, tol, target, @@ -1113,16 +823,16 @@ class MaterialReactivityController(ReactivityController): self.material = self._get_material(material) - check_type("material vector", mat_vector, dict, str) - for nuc in mat_vector: + 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 - if not isclose(sum(mat_vector.values()), 1.0, abs_tol=0.01): - # Normalize material elements vector - sum_values = sum(mat_vector.values()) - for elm in mat_vector: - mat_vector[elm] /= sum_values - self.mat_vector = mat_vector + 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 @@ -1188,110 +898,36 @@ class MaterialReactivityController(ReactivityController): # right amount root = self._search_for_keff(0) - # Update concentration vector and volumes with new value - volumes = self._calculate_volumes(root) + # 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. -class RefuelMaterialReactivityController(MaterialReactivityController): - """Reactivity control material-based class for refuelling (addition or - removal) scheme. + Parameters + ---------- + param : float + grams or atoms or cc of material_to_mix to convert to cc - A user doesn't need to call this class directly. - Instead an instance of this class is automatically created by calling - :meth:`openmc.deplete.Integrator.add_reactivity_control` method from an - integrator class, such as :class:`openmc.deplete.CECMIntegrator`. + """ + 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 - .. versionadded:: 0.14.1 + self.material_to_mix.volume = param * multiplier - Parameters - ---------- - material : openmc.Material or int or str - OpenMC Material identifier to where apply batch wise scheme - operator : openmc.deplete.Operator - OpenMC operator object - attrib_name : str - Material attribute name - mat_vector : dict - Dictionary of material composition to parameterize, where a pair key value - represents a nuclide and its weight fraction, respectively. - 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. - density_treatment : str - Whether or not to keep constant volume or density after a depletion step - before the next one. - Default to 'constant-volume' - 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 - mat_vector : dict - Dictionary of material composition to parameterize, where a pair key value - represents a nuclide and its weight fraction, respectively. - attrib_name : str - Material attribute name - - """ - - def __init__( - self, - material, - attrib_name, - operator, - mat_vector, - bracket, - bracket_limit, - density_treatment="constant-volume", - bracketed_method="brentq", - tol=0.01, - target=1.0, - print_iterations=True, - search_for_keff_output=True, - ): - - super().__init__( - material, - operator, - mat_vector, - bracket, - bracket_limit, - density_treatment, - bracketed_method, - tol, - target, - print_iterations, - search_for_keff_output, - ) - - # Not needed but used for consistency with other classes - check_value("attrib_name", attrib_name, "refuel") - self.attrib_name = attrib_name def _model_builder(self, param): """Callable function which builds a model according to a passed @@ -1300,14 +936,10 @@ class RefuelMaterialReactivityController(MaterialReactivityController): Builds the parametric model to be passed to the :meth:`openmc.search.search_for_keff` method. - The parametrization can either be at constant volume or constant - density, according to user input. - Default is constant volume. - Parameters ---------- param : float - Model function variable, total grams of material to add or remove + Model function variable, cc of material_to_mix Returns ------- @@ -1317,68 +949,68 @@ class RefuelMaterialReactivityController(MaterialReactivityController): for rank in range(comm.size): number_i = comm.bcast(self.operator.number, root=rank) - for mat in number_i.materials: + for mat_id in number_i.materials: nuclides = [] densities = [] - if int(mat) == self.material.id: + if int(mat_id) == self.material.id: + self._set_mix_material_volume(param) - if self.density_treatment == "constant-density": - vol = number_i.get_mat_volume(mat) + ( - param / self.material.get_mass_density() + 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) - elif self.density_treatment == "constant-volume": - vol = number_i.get_mat_volume(mat) + # 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: - if nuc in self.mat_vector: - # units [#atoms/cm-b] - val = 1.0e-24 * ( - number_i.get_atom_density(mat, nuc) - + param - / atomic_mass(nuc) - * AVOGADRO - * self.mat_vector[nuc] - / vol - ) + 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] - else: - # get normalized atoms density in [atoms/b-cm] - val = 1.0e-24 * number_i[mat, nuc] / vol + atoms_per_bcm = 1.0e-24 * atoms / vol - if val > 0.0: + if atoms_per_bcm > 0.0: nuclides.append(nuc) - densities.append(val) + 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] - val = 1.0e-24 * number_i.get_atom_density(mat, nuc) + atoms_per_bcm = 1.0e-24 * number_i.get_atom_density(mat_id, nuc) - if val > 0.0: + if atoms_per_bcm > 0.0: nuclides.append(nuc) - densities.append(val) + densities.append(atoms_per_bcm) - # set nuclides and densities to the in-memory model - openmc.lib.materials[int(mat)].set_densities(nuclides, densities) + # 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 _calculate_volumes(self, res): + def _adjust_volumes(self, res): """Uses :meth:`openmc.deplete.reactivity_control._search_for_keff` - solution as grams of material to add or remove to calculate new material - volume. + solution as cc to update depletable volume if `self.dilute` attribute + is set to True. Parameters ---------- res : float - Solution in grams of material, coming from + Solution in cc of material_to_mix, coming from :meth:`openmc.deplete.reactivity_control._search_for_keff` Returns @@ -1390,12 +1022,10 @@ class RefuelMaterialReactivityController(MaterialReactivityController): number_i = self.operator.number volumes = {} - for mat in self.local_mats: - if int(mat) == self.material.id: - if self.density_treatment == "constant-density": - volumes[mat] = number_i.get_mat_volume(mat) + ( - res / self.material.get_mass_density() - ) - elif self.density_treatment == "constant-volume": - volumes[mat] = number_i.get_mat_volume(mat) + 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 From fa13a813d92cde29707f4192a00fe35fe5535cb8 Mon Sep 17 00:00:00 2001 From: church89 Date: Thu, 6 Jun 2024 13:21:22 +0200 Subject: [PATCH 2/2] update tests --- .../ref_depletion_with_refuel.h5 | Bin 38432 -> 39352 bytes .../ref_depletion_with_rotation.h5 | Bin 38432 -> 39352 bytes .../ref_depletion_with_translation.h5 | Bin 38432 -> 39352 bytes .../deplete_with_reactivity_control/test.py | 51 ++++-- .../test_deplete_reactivity_control.py | 170 +++++++++--------- 5 files changed, 125 insertions(+), 96 deletions(-) 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 26bd4839da3a1c67c9127ed0fa5ee7126340f53f..83a19702b0ef95740453e39a79ae1c44bf364e2f 100644 GIT binary patch delta 1323 zcmZuveM}o=81KDx?Oi9OgfXgaK!rgkjJCp{d>q%7^0m2H$hI6VQ?$aGx(%3fVG{xg zh}ALa@eFNA*vc?M49>DsN7Q9>x-DQPH4HUOBTF;_(^yf2F~_{_TFM_@a_{>*&+qqq z-0wNp02*SEU5SRA%oA@ROTpl$B81gJk0_S$>ex78lMvBH#uEd>A{;68ersApV>ApJ zob06dpNBBhpl1P-&;tIcOJ$ykA>3)a*AFAYSP7jTz%pP|oX8f3Pe)(`Fp)*v`HeWX z*P#%CT46c6NeEL3IJgojdM!#}htk4&lf7t->6ItbvYQnPtPKEp?|Gv1CA^d zXJJYx&4`W&ZH=#L0m6_Zso))f*1&k)2A}RFq;2Pl^n$BQj+q#bim<%HQ%q1>&HY3& zWOu-%F%grnSiyxb%*A*ziJ<_;#dG9nD#+=nzDBTd1>QB_dqIb22Ux2U881#WIRlNM zRv67TnWxtx%L9Xi6lhRs%SLocd~QMQk<}@ zjz+kH^#TnsZtnZy#?QJN-(;4Q#tl^$$Mh4O#xE~)n43?}83QcAz?JJA)}KBngZ}N~ zy6-+#oC_>lAn1OvC-DaCNyH@;i!Mh9n;|EvxhWQPtzuoT(;-!Lfnu(ee z_c$4MqVnI8s&RkG=uU&Lxcs#G~g8aHZhYJ;CC*sT~Dyc-shWmD~ez9Vc% zmfP~jb|qjdvIzbfm`E@0A1BW&Z?+wdY$01OpE=nzbcGauWwu#Wv*d7dYIgcZtz^T} z3_8+QK@y*=Q@A-*;=Oe#gGq~|B5M5%c5wL*D(tDfmroZRr5@~0KCc7csdW$grD=KX L>+K6uuTuX3)en;< delta 1018 zcmZvYe`r%z6vyAa&!zW0Gj(++I^QGy5y{$nVj!C|*oUt0ac3-{c6 zzUOn!J@<@KoV<*~R&0((KKU6l9vF9jJ)Hg?fte>Tgasr3a_y#O_9TmPpkNcEbOfbT z(F9pt5a(PGB~E9G@=K6Bth_{gfYDaRI=I6cVbm@_&RNQZYcx3kw~H)PHe`z-Q$`?J z*4&r)h>co9y=zz)svyu-(cBlZvLI9uNZVvOT9L(Yah(mO+(8KI2-qrbMQ(OLn=J?d z2LZtm!m^Wo>sWzt&<;`OE=V;jp*su?jC0lC;=f`sLqPGDtRJ##J*2C(vs62nOrbYa z@JX$9v8sQ20cMu)Z!>}4(8wiS8W5o9!+MNdaA{>L|27L&)g3_`MCwec9j@~%5BfBI z{(*(f0g)1FyhpOA_O`TXdGho3P}Ev~2H{)_jJW5p?EoF|=n)PbgE8+MPW}L^8~WZ1 zrOi(u=Qs29S_y{vG@=+G)PYb-|`Ek1OeqbKDItIo9)g(FH0Cl|;t>}jmQ!39!joHlSXg_xjLvEXcs z-f7h><+AUnCY+DaSM7S9%NJClvL=qi3w^>twuA?`>XzPq_2sEH<;JbiH=_VZ*>cP zZ)f`;Gh4eDFaP)h{qK*Hi~iPo`wsnbY5yWW`>gF=EHSVcKRRiQemU`Xt;nMuCeL0>y5ZRd zh9`XWjOUv*_EYm8-!}bdWMd%iYd@dKZt0$d9-S$YJ`J(6B}?ba3s-8*akwk}4?i+B A4FCWD 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 e82060e5d037f3433336c4ac836a42fffad31c92..3831319390ece1aebcc4ee14e36bc11b913ad9e0 100644 GIT binary patch delta 1284 zcmZuveM}o=81KD>_O4StW=T(FvAjFmjUp}Y=W6vZuv1GgoUjk9G8rKQm1>r- zz*nX%ewU zUqoace0@WR z?`b%JKYnNZYFvLk?xMZ1WhCXLU!K=HD}d=S%`TIW6xz z_h8Ye+S9nJ_DQAd`-z3eb5WYoURHlR)~P8i9IO57?1riI==x9g;&N5#aJ%N``Meog z@?h=5_UntZFfrbdzoj;BsB5MNuMC)TzWcknZ_sNFy6~Cd;A^u&*p**tP1pF3VON&= zbe#)1C1G1GcE37_;|D^Y{ptKKht6#^glw5A`t2(&?_vmlq~%AJF1hh3v!U|2AA0bj zPL_~(&4rJ?FlFh>S-}sh&W{(*p1@YRV10V19;f$Be6%%nAJ4o*#}oU#IJYiwxjl#Q i7$YPjxtn7x#-}$|{$T>o>G)^-*7PQQuwYxMaQy>ragiPX delta 987 zcmZvYZAep57{~89H|DuF(Q9Q|WLspqX0)kOYc;2v%UeMW0~^wdmfC_qD~gID>4k+w zZ9Z+O5Ns(FVxgTZBdLvEd?8}?g4U-&5&KXyDoP4Qx9-$}4xICx=l}f9|2*gPQ|#-; zE;Uwp-J0{5F~if+J-wk52vj_Q9xNaMkYPJju}v(BLJ%t^N$02fcJW zVPXandxPH*Z)4fM{Vl=PHtA8w@98Q!zP`0^v*ozRc>5)ng4yklivLhlU2+!T=nm*BjbU8_?K3AK>^KaMEivpn2^r-be?FzrN+4=g@g=ea z!>o1(7ec^FT68MT5f{VYw<AgBt zuk*ZS7Nd~yYrUrJ(aY0D+vtIT+=l_T!FvCX!%L;TWE#-&Etn+9EDyX%K0UUKjIJkRg< zyubH;&pEH61rFJiXu-i2yp5Qu+0P{i>j#%4i}e|}Z1J!d(UmNp30?_~l*VwhA)zrE zW=#(6w)FlZm}Qp#1C!Ar;psof9?v8^Y9VvOk~m97ABtcaSeE8-wJ96BF#_1sn$Wr` z&6=_A7emp0W;$04>p3_$bDC!urTASX!}$t*XqTn8M6H9#gKCUs`J-TyNH7VgR7;oz zj+ye=n-=IQqQO&$LswxfN2h2JqSbI^UpoYnz7UlD^mIL>Rio z8QgH&Mgn3Uf`l!KaNTVy1UN^ES)Y{Xw+nnctbm*jZh1f=BZmJ~%|LDkw>^@K6L0>f z>Tz&SoXf+9B6^|%Pu<0)ALqzuzLQ%gQ+@nVWQrou@6E4WJKI_k!~0eg3A660Sf0rK zQr+-vTZ}*5aoTU(x)j47?|doo46{n=omcyg?i(d9d^41jkjm@G#evp4omGCaxA8I8 z7v@pYpc9=vhPRM|_1bQ6qo4fF=>O#2jcxMO<#_{l|1LFKlVP3zu!n-cdui}sTskDb zyhi0*+Qa9&vb~wqttN+}*le{=^q#&^Y__~r^%RF2&EN0PEfdeF%t^nl_mP`ln|@Y~ z;R7RUCiR(uS2e++$_ZwV+ch^?`ASRsGv{2}l~EUDTe$~LU5~r6W_%qz+7ovrYjbGG zn{JNVa=z)1{%_UA*k_N849Y+4joUJ?GZj28a>wwm%v5!L$tt;+^zxJchKY4iXi!0sB+vH?>y5h5>!MoJS gBl@u`Q!i5wZ?3p!fWY|P#IQA+)DInV{n_4s00;AeApigX delta 986 zcmZvYZAep57{_DrxE(1CNF^ZcLR`Jd;UUW&ul zu-}Lcp}>MOSmS_a_HEtrQ3M7ZfiBD;4v={RH84#KN`QoUh-wLnDk2YYtrsWE0VzZW z(vlO9*e!iWI*&#j`wHL*lLx)U93;y0tiN29JnYR4;ln7cuELxB|a>c(E}w@ zG4>WiOIZy>S4^PmGD|Q{l!94T$6yTs$(_5Xt61$pE>)vm4VqZ0HicF_^ic_} zq0uZ0NO=%-lNQ0i`ch&HB5oeI9^Ytz~IUn-B2 zLHc87hK^09RHCpvOG2qWVF!~-JqP0~s|l*oD3FV}5UoC*L0IySst?hVXXlZ_v81`{ z?9(1g%0$y}uaXgPNJs;aWym_3@uH7@GX782p zx_4&VI_E3d($tbuzI<-8WvFrG;F5)sW!aV!?jF}W-f~1(RqlM-jh4rg zOP_w6n6MT6th?9V3H7acmBSwfZFiy{o6bcI(Lxhn>(UoF*WaJsZ@F9KOgCAdrxCle zn@lR0v9{2e{rc98>z^Jv-ix!pUTz+6)choe%0>j&4f)&ENKd;fb7)iL>Btv|X!(GA T9D;4R>7(|I^G%M<4Tk71Er}@4 diff --git a/tests/regression_tests/deplete_with_reactivity_control/test.py b/tests/regression_tests/deplete_with_reactivity_control/test.py index 006187417a..c9d0023f20 100644 --- a/tests/regression_tests/deplete_with_reactivity_control/test.py +++ b/tests/regression_tests/deplete_with_reactivity_control/test.py @@ -10,6 +10,10 @@ 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 @@ -78,13 +82,21 @@ def model(): return openmc.Model(geometry, materials, settings) -@pytest.mark.parametrize("obj, attribute, bracket_limit, axis, vec, ref_result", [ - ('trans_cell', 'translation', [-40,40], 2, None, 'depletion_with_translation'), - ('rot_cell', 'rotation', [-90,90], 2, None, 'depletion_with_rotation'), - ('f', 'refuel', [-100,100], None, {'U235':0.9, 'U238':0.1}, 'depletion_with_refuel') +@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_limit, - axis, vec, ref_result): +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) @@ -92,16 +104,23 @@ def test_reactivity_control(run_in_tmpdir, model, obj, attribute, bracket_limit, integrator = openmc.deplete.PredictorIntegrator( op, [1], 174., timestep_units = 'd') - kwargs = {'bracket': [-4,4], 'bracket_limit':bracket_limit, - 'tol': 0.1,} - - if vec is not None: - kwargs['mat_vector']=vec - - if axis is not None: - kwargs['axis'] = axis - - integrator.add_reactivity_control(obj, attribute, **kwargs) + 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 diff --git a/tests/unit_tests/test_deplete_reactivity_control.py b/tests/unit_tests/test_deplete_reactivity_control.py index 0f087be6ea..a202226db8 100644 --- a/tests/unit_tests/test_deplete_reactivity_control.py +++ b/tests/unit_tests/test_deplete_reactivity_control.py @@ -9,9 +9,8 @@ import openmc import openmc.lib from openmc.deplete import CoupledOperator from openmc.deplete import ( - GeometricalCellReactivityController, - TemperatureCellReactivityController, - RefuelMaterialReactivityController + CellReactivityController, + MaterialReactivityController ) CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml" @@ -74,62 +73,52 @@ def integrator(operator): return openmc.deplete.PredictorIntegrator( operator, [1,1], 0.0, timestep_units = 'd') -@pytest.mark.parametrize("case_name, obj, attribute, bracket, limit, axis, vec", [ - ('cell translation','universe_cell', 'translation', [-1,1], [-10,10], 2, None), - ('cell rotation', 'universe_cell', 'rotation', [-1,1], [-10,10], 2, None), - ('single cell temperature', 'fuel_cell', 'temperature', [-1,1], [-10,10], 2, None ), - ('multi cell', 'universe_cell', 'temperature', [-1,1], [-10,10], 2, None ), - ('material refuel', 'fuel', 'refuel', [-1,1], [-10,10], None, {'U235':0.1, 'U238':0.9}), - ('invalid_1', 'universe_cell', 'refuel', [-1,1], [-10,10], None, {'U235':0.1, 'U238':0.9}), - ('invalid_2', 'fuel', 'temperature', [-1,1], [-10,10], 2, None ), - ('invalid_3', 'fuel', 'translation', [-1,1], [-10,10], 2, None), +@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, vec): + 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 - kwargs = {'bracket': bracket, 'bracket_limit':limit, 'tol': 0.1,} - - if vec is not None: - kwargs['mat_vector']=vec - - if axis is not None: - kwargs['axis'] = axis - - if case_name == "invalid_1": - with pytest.raises(ValueError) as e: - integrator.add_reactivity_control(obj, attribute, **kwargs) - assert str(e.value) == 'Unable to set "Material name" to "universe_cell" '\ - 'since it is not in "[\'fuel\', \'water\']"' - elif case_name == "invalid_2": - with pytest.raises(ValueError) as e: - integrator.add_reactivity_control(obj, attribute, **kwargs) - assert str(e.value) == 'Unable to set "Cell name exists" to "fuel" since '\ - 'it is not in "[\'fuel_cell\', \'universe_cell\', \'\', \'\']"' - - elif case_name == "invalid_3": - with pytest.raises(ValueError) as e: - integrator.add_reactivity_control(obj, attribute, **kwargs) - assert str(e.value) == 'Unable to set "Cell name exists" to "fuel" since '\ - 'it is not in "[\'fuel_cell\', \'universe_cell\', \'\', \'\']"' else: - integrator.add_reactivity_control(obj, attribute, **kwargs) - if attribute in ('translation','rotation'): - assert integrator.reactivity_control.universe_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 + 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 - elif attribute == 'refuel': - assert integrator.reactivity_control.mat_vector == vec - - assert integrator.reactivity_control.attrib_name == attribute - 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 + 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), @@ -140,9 +129,14 @@ def test_cell_methods(run_in_tmpdir, model, operator, integrator, obj, attribute """ Test cell base class internal method """ - kwargs = {'bracket':[-1,1], 'bracket_limit':[-10,10], 'axis':2, 'tol':0.1} - - integrator.add_reactivity_control(obj, attribute, **kwargs) + 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() @@ -150,63 +144,79 @@ def test_cell_methods(run_in_tmpdir, model, operator, integrator, obj, attribute integrator.reactivity_control._set_cell_attrib(value_to_set) assert integrator.reactivity_control._get_cell_attrib() == value_to_set - vol = integrator.reactivity_control._calculate_volumes() + vol = integrator.reactivity_control._adjust_volumes() integrator.reactivity_control._update_x_and_set_volumes(operator.number.number, vol) - for cell in integrator.reactivity_control.universe_cells: + 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("nuclide, atoms_to_add", [ - ('U238', 1.0e22), - ('Xe135', 1.0e21) +@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, nuclide, - atoms_to_add): +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. """ - kwargs = {'bracket':[-1,1], 'bracket_limit':[-10,10], 'mat_vector':{}} - - integrator.add_reactivity_control('fuel', 'refuel', **kwargs) + 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() - #Increase number of atoms of U238 in fuel by fix amount and check the - # volume increase at constant-density - #extract fuel material from model materials + # 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[nuclide] - vol = operator.number.get_mat_volume(str(mat.id)) - operator.number.number[mat_index][nuc_index] += atoms_to_add - integrator.reactivity_control._update_volumes() + nuc_index = operator.number.index_nuc['U235'] + vol_before_mix = operator.number.get_mat_volume(str(mat.id)) - vol_to_compare = vol + (atoms_to_add * openmc.data.atomic_mass(nuclide) /\ - openmc.data.AVOGADRO / mat.density) + # 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) - assert operator.number.get_mat_volume(str(mat.id)) == pytest.approx(vol_to_compare) + 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(nuclide) - dens_to_compare = 1.0e-24 * operator.number.get_atom_density(str(mat.id), nuclide) + 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) - volumes = {str(mat.id): vol + 1} - new_x = integrator.reactivity_control._update_x_and_set_volumes(x, volumes) - dens_to_compare = 1.0e24 * volumes[str(mat.id)] *\ + # 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)) == volumes[str(mat.id)] + assert operator.number.get_mat_volume(str(mat.id)) == vol_after_mix[str(mat.id)] openmc.lib.finalize()