diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index 8a9509e900..c86a93c8a1 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -17,6 +17,7 @@ from .stepresult import * from .results import * from .integrators import * from .transfer_rates import * +from .batchwise import * from . import abc from . import cram from . import helpers diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 60774871c0..7f06fd332f 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -18,14 +18,16 @@ from warnings import warn from numpy import nonzero, empty, asarray from uncertainties import ufloat -from openmc.checkvalue import check_type, check_greater_than, PathLike +from openmc.checkvalue import checkvalue, check_type, check_greater_than, PathLike from openmc.mpi import comm from .stepresult import StepResult from .chain import Chain from .results import Results from .pool import deplete from .transfer_rates import TransferRates - +from openmc import Material, Cell, Universe +from .batchwise import (BatchwiseCellGeometrical, BatchwiseCellTemperature, + BatchwiseMaterialRefuel) __all__ = [ "OperatorResult", "TransportOperator", @@ -636,6 +638,7 @@ class Integrator(ABC): self.source_rates = asarray(source_rates) self.transfer_rates = None + self.batchwise = None if isinstance(solver, str): # Delay importing of cram module, which requires this file @@ -763,6 +766,14 @@ class Integrator(ABC): return (self.operator.prev_res[-1].time[-1], len(self.operator.prev_res) - 1) + def _get_bos_from_batchwise(self, step_index, bos_conc): + """Get BOS from criticality batch-wise control + """ + x = deepcopy(bos_conc) + # Get new vector after keff criticality control + x, root = self.batchwise.search_for_keff(x, step_index) + return x, root + def integrate( self, final_step: bool = True, @@ -797,6 +808,11 @@ class Integrator(ABC): # Solve transport equation (or obtain result from restart) if i > 0 or self.operator.prev_res is None: + # Update geometry/material according to batchwise definition + if self.batchwise and source_rate != 0.0: + n, root = self._get_bos_from_batchwise(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(i, source_rate, n) @@ -810,9 +826,8 @@ class Integrator(ABC): # 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 @@ -822,9 +837,13 @@ class Integrator(ABC): # solve) if output and final_step and comm.rank == 0: print(f"[openmc.deplete] t={t} (final operator evaluation)") + if self.batchwise and source_rate != 0.0: + n, root = self._get_bos_from_batchwise(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) + source_rate, self._i_res + len(self), proc_time, root) self.operator.write_bos_data(len(self) + self._i_res) self.operator.finalize() @@ -857,6 +876,30 @@ class Integrator(ABC): self.transfer_rates.set_transfer_rate(material, components, transfer_rate, transfer_rate_units, destination_material) + def add_batchwise(self, obj, attr, **kwargs): + """Add batchwise operation to integrator scheme. + + Parameters + ---------- + attr : str + Type of batchwise operation to add. `Trans` stands for geometrical + translation, `refuel` for material refueling and `dilute` for material + dilute. + **kwargs + keyword arguments that are passed to the batchwise class. + + """ + check_value('attribute', attr, ('translation', 'rotation', 'temperature', 'refuel')) + if attr in ('translation', 'rotation'): + batchwise = BatchwiseCellGeometrical + elif attr == 'temperature': + batchwise = BatchwiseCellTemperature + elif attr == 'refuel': + batchwise = BatchwiseMaterialRefuel + + self.batchwise = batchwise.from_params(obj, attr, self.operator, + self.operator.model, **kwargs) + @add_params class SIIntegrator(Integrator): r"""Abstract class for the Stochastic Implicit Euler integrators diff --git a/openmc/deplete/batchwise.py b/openmc/deplete/batchwise.py index 03dfcf9f04..7136f7e6b2 100644 --- a/openmc/deplete/batchwise.py +++ b/openmc/deplete/batchwise.py @@ -1,18 +1,77 @@ from abc import ABC, abstractmethod -from collections.abc import Iterable -from openmc.search import _SCALAR_BRACKETED_METHODS, search_for_keff -from openmc import Materials, Material, Cell -from openmc.data import atomic_mass, AVOGADRO, ELEMENT_SYMBOL - +from copy import deepcopy +from numbers import Real +from warnings import warn import numpy as np -import os -import h5py + +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_iterable_type, check_length) class Batchwise(ABC): + """Abstract class defining a generalized batchwise scheme. + In between transport and depletion steps batchwise operations can be added + with the aim of maintain a system critical. These operations act on OpenMC + Cell or Material, parametrizing one or more coefficients while running a + search_for_kef algorithm. + + Specific classes for running batchwise depletion calculations are + implemented as derived class of Batchwise + . + .. versionadded:: 0.13.4 + + Parameters + ---------- + operator : openmc.deplete.Operator + OpenMC operator object + model : openmc.model.Model + OpenMC model object + bracket : list of float + Bracketing interval to search for the solution, always 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 + Wether ot not to keep contant 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 + Wheter or not to print `search_for_keff` iterations. + Default to True + search_for_keff_output : Bool, Optional + Wheter or not to print transport iterations during `search_for_keff`. + 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, model, bracket, bracket_limit, - bracketed_method='brentq', tol=0.01, target=1.0, - print_iterations=True, search_for_keff_output=True): + density_treatment = 'constant-volume', 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 @@ -20,6 +79,10 @@ class Batchwise(ABC): self.model = model self.geometry = model.geometry + check_value('density_treatment', density_treatment, + ('constant-density','constant-volume')) + self.density_treatment = density_treatment + check_iterable_type('bracket', bracket, Real) check_length('bracket', bracket, 2) check_less_than('bracket values', bracket[0], bracket[1]) @@ -29,12 +92,18 @@ class Batchwise(ABC): check_length('bracket_limit', bracket_limit, 2) check_less_than('bracket limit values', bracket_limit[0], bracket_limit[1]) - self.bracket_limit = bracket_limit + check_value('bracketed_method', bracketed_method, + _SCALAR_BRACKETED_METHODS) self.bracketed_method = bracketed_method + + check_type('tol', tol, Real) self.tol = tol + + check_type('target', target, Real) self.target = target + self.print_iterations = print_iterations self.search_for_keff_output = search_for_keff_output @@ -46,7 +115,7 @@ class Batchwise(ABC): def bracketed_method(self, value): check_value('bracketed_method', value, _SCALAR_BRACKETED_METHODS) if value != 'brentq': - warn('brentq bracketed method is recomended here') + warn('WARNING: brentq bracketed method is recommended') self._bracketed_method = value @property @@ -69,48 +138,33 @@ class Batchwise(ABC): @abstractmethod def _model_builder(self, param): - - def _get_cell_id(self, val): - """Helper method for getting cell id from cell instance or cell name. + """ + Builds the parametric model to be passed to the `search_for_keff` + algorithm. + Callable function which builds a model according to a passed + parameter. This function must return an openmc.model.Model object. Parameters ---------- - val : Openmc.Cell or str or int representing Cell + param : parameter + model function variable Returns ------- - id : str - Cell id + _model : openmc.model.Model + OpenMC parametric model """ - if isinstance(val, Cell): - check_value('Cell id', val.id, [cell.id for cell in \ - self.geometry.get_all_cells().values()]) - val = val.id - - elif isinstance(val, str): - if val.isnumeric(): - check_value('Cell id', val, [str(cell.id) for cell in \ - self.geometry.get_all_cells().values()]) - val = int(val) - else: - check_value('Cell name', val, [cell.name for cell in \ - self.geometry.get_all_cells().values()]) - - val = [cell.id for cell in \ - self.geometry.get_all_cells().values() \ - if cell.name == val][0] - - elif isinstance(val, int): - check_value('Cell id', val, [cell.id for cell in \ - self.geometry.get_all_cells().values()]) - - else: - ValueError(f'Cell: {val} is not recognized') - - return val def _search_for_keff(self, val): """ Perform the criticality search for a given parametric model. - It supports geometrical or material based `search_for_keff`. + It supports cell or material based `search_for_keff`. + If the solution lies off the initial bracket, the method will iteratively + adapt the bracket to be able to run `search_for_keff` effectively. + The ratio between the bracket and the returned keffs values will be + the scaling factor to iteratively adapt the brackt unitl a valid solution + is found. + If the adapted bracket is moved too far off, i.e. above or below user + inputted bracket limits interval, the algorithm will stop and the closest + limit value will be used. Parameters ---------- val : float @@ -133,14 +187,15 @@ class Batchwise(ABC): # Run until a search_for_keff root is found or ouf ot limits root = None - while res == None: + while root == None: search = search_for_keff(self._model_builder, - bracket = [bracket[0] + val, bracket[1] + val], - tol = tol, - bracketed_method = self.bracketed_method, - target = self.target, - print_iterations = self.print_iterations, - run_args = {'output': self.search_for_keff_output}) + bracket = np.array(bracket) + val, + tol = 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: @@ -153,8 +208,8 @@ class Batchwise(ABC): else: # Set res with the closest limit and continue arg_min = abs(np.array(self.bracket_limit) - res).argmin() - warn('WARNING: Search_for_keff returned root out of '\ - 'bracket limit. Set root to {:.2f} and continue.' + warn("WARNING: 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] @@ -165,8 +220,8 @@ class Batchwise(ABC): if all(self.bracket_limit[0] < guess < self.bracket_limit[1] \ for guess in guesses): #Simple method to iteratively adapt the bracket - print('INFO: Function returned values below or above ' \ - 'target. Adapt bracket...') + print("Search_for_keff returned values below or above " + "target. Trying to iteratively adapt bracket...") # if the bracket ends up being smaller than the std of the # keff's closer value to target, no need to continue- @@ -181,7 +236,7 @@ class Batchwise(ABC): # 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 guess[np.argmax(keffs)] > guess[np.argmin(keffs)]: + if guesses[np.argmax(keffs)] > guesses[np.argmin(keffs)]: dir = 1 else: dir = -1 @@ -189,7 +244,7 @@ class Batchwise(ABC): bracket[np.argmax(keffs)] += grad * (self.target - \ max(keffs).n) * dir else: - if guess[np.argmax(keffs)] > guess[np.argmin(keffs)]: + if guesses[np.argmax(keffs)] > guesses[np.argmin(keffs)]: dir = -1 else: dir = 1 @@ -200,51 +255,50 @@ class Batchwise(ABC): else: # Set res with closest limit and continue arg_min = abs(np.array(self.bracket_limit) - guesses).argmin() - warn('WARNING: Adaptive iterative bracket went off '\ - 'bracket limits. Set root to {:.2f} and continue.' + warn("WARNING: Adaptive iterative bracket went off " + "bracket limits. Set root to {:.2f} and continue." .format(self.bracket_limit[arg_min])) root = self.bracket_limit[arg_min] else: - raise ValueError(f'ERROR: Search_for_keff output is not valid') + raise ValueError('ERROR: Search_for_keff output is not valid') return root - def _save_res(self, type, step_index, root): + def _update_volumes(self): """ - Save results to msr_results.h5 file. - Parameters - ---------- - type : str - String to characterize geometry and material results - step_index : int - depletion time step index - root : float or dict - Root of the search_for_keff function - """ - filename = 'msr_results.h5' - kwargs = {'mode': "a" if os.path.isfile(filename) else "w"} - - if comm.rank == 0: - with h5py.File(filename, **kwargs) as h5: - name = '_'.join([type, str(step_index)]) - if name in list(h5.keys()): - last = sorted([int(re.split('_',i)[1]) for i in h5.keys()])[-1] - step_index = last + 1 - h5.create_dataset('_'.join([type, str(step_index)]), data=root) - - def _update_volumes_after_depletion(self, x): - """ - After a depletion step, both materials volume and density change, due to + Update number volume. + After a depletion step, both material volume and density change, due to decay, transmutation reactions and transfer rates, if set. At present we lack an implementation to calculate density and volume - changes due to the different molecules speciation. Therefore, the assumption - is to consider the density constant and let the material volume - vary with the change in nuclides concentrations. - The method uses the nuclides concentrations coming from the previous Bateman - solution and calculates a new volume, keeping the mass density of the material - constant. It will then assign the volumes to the AtomNumber instance. + 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 updates the materials volume in AtomNumber dataset. + Parameters + ---------- + x : list of numpy.ndarray + Total atom concentrations + """ + 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', "_update_volumes" will + update material volumes, keeping the material total density constant, in + AtomNumber to renormalize the atom densities before assign them to the + model in memory. Parameters ---------- x : list of numpy.ndarray @@ -252,93 +306,8 @@ class Batchwise(ABC): """ self.operator.number.set_density(x) - for rank in range(comm.size): - number_i = comm.bcast(self.operator.number, root=rank) - - for i, mat in enumerate(number_i.materials): - # Total nuclides density - density = 0 - for nuc in number_i.nuclides: - # total number of atoms - val = number_i[mat, nuc] - # obtain nuclide density in atoms-g/mol - density += val * atomic_mass(nuc) - # Get mass dens from beginning, intended to be held constant - rho = openmc.lib.materials[int(mat)].get_density('g/cm3') - number_i.volume[i] = density / AVOGADRO / rho - -class BatchwiseCell(Batchwise): - - def __init__(self, operator, model, cell, attrib_name, axis, bracket, - bracket_limit, bracketed_method='brentq', tol=0.01, target=1.0, - print_iterations=True, search_for_keff_output=True): - - super().__init__(operator, model, bracket, bracket_limit, - bracketed_method, tol, target, print_iterations, - search_for_keff_output) - - self.cell_id = super()._get_cell_id(cell_id_or_name) - check_value('attrib_name', attrib_name, - ('rotation', 'translation')) - self.attrib_name = attrib_name - - #index of cell directionnal axis - check_value('axis', axis, (0,1,2)) - self.axis = axis - - # Initialize vector - self.vector = np.zeros(3) - - # materials that fill the attribute cell, if depletables - self.cell_material_ids = [cell.fill.id for cell in \ - self.geometry.get_all_cells()[self.cell_id].fill.cells.values() \ - if cell.fill.depletable] - - def _get_cell_attrib(self): - """ - Get cell attribute coefficient. - Returns - ------- - coeff : float - cell coefficient - """ - for cell in openmc.lib.cells.values(): - if cell.id == self.cell_id: - if self.attrib_name == 'translation': - return cell.translation[self.axis] - elif self.attrib_name == 'rotation': - return cell.rotation[self.axis] - - def _set_cell_attrib(self, val): - """ - Set cell attribute to the cell instance. - Attributes are only applied to cells filled with a universe - Parameters - ---------- - var : float - Surface coefficient to set - geometry : openmc.model.geometry - OpenMC geometry model - attrib_name : str - Currently only translation is implemented - """ - self.vector[self.axis] = val - for cell in openmc.lib.cells.values(): - if cell.id == self.cell_id or cell.name == self.cell_id: - setattr(cell, self.attrib_name, self.vector) - - def _update_materials(self, x): - """ - Assign concentration vectors from Bateman solution at previous - timestep to the in-memory model materials, after having recalculated the - material volume. - - Parameters - ---------- - x : list of numpy.ndarray - Total atom concentrations - """ - super()._update_volumes_after_depletion(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) @@ -351,56 +320,210 @@ class BatchwiseCell(Batchwise): # 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 > self.atom_density_limit: + if val > 0.0: nuclides.append(nuc) densities.append(val) - #set nuclide densities to model in memory (C-API) + # Update densities on C API side openmc.lib.materials[int(mat)].set_densities(nuclides, densities) - def _update_volumes(self): - openmc.lib.calculate_volumes() - res = openmc.VolumeCalculation.from_hdf5('volume_1.h5') + def _update_x_and_set_volumes(self, x, volumes): + """ + Update x with volumes, before assign them to AtomNumber materials. + Parameters + ---------- + x : list of numpy.ndarray + Total atoms concentrations + volumes : dict + Returns + ------- + x : list of numpy.ndarray + Updated total atoms concentrations + """ number_i = self.operator.number - for mat_idx, mat_id in enumerate(self.local_mats): - if mat_id in self.cell_material_ids: - number_i.volume[mat_idx] = res.volumes[int(mat_id)].n - def _update_x(self): - number_i = self.operator.number - for mat_idx, mat_id in enumerate(self.local_mats): - if mat_id in self.cell_material_ids: + for mat_idx, mat in enumerate(self.local_mats): + + if mat in volumes: + res_vol = volumes[mat] + # Normalize burbable nuclides in x vector without cross section data for nuc_idx, nuc in enumerate(number_i.burnable_nuclides): - x[mat_idx][nuc_idx] = number_i.volume[mat_idx] * \ - number_i.get_atom_density(mat_idx, nuc) + if nuc not in self.operator.nuclides_with_data: + # normalzie 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 BatchwiseCell(Batchwise): + """Abstract class holding batchwise cell-based functions. + + Specific classes for running batchwise depletion calculations are + implemented as derived class of BatchwiseCell. + + .. versionadded:: 0.13.4 + + Parameters + ---------- + cell : openmc.Cell or int or str + OpenMC Cell identifier to where apply batchwise scheme + operator : openmc.deplete.Operator + OpenMC operator object + model : openmc.model.Model + OpenMC model object + bracket : list of float + Bracketing interval to search for the solution, always 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 + Wether ot not to keep contant 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 + Wheter or not to print `search_for_keff` iterations. + Default to True + search_for_keff_output : Bool, Optional + Wheter 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 batchwise scheme + cell_materials : list of openmc.Material + Depletable materials that fill the Cell Universe. Only valid for + translation or rotation atttributes + """ + def __init__(self, cell, operator, model, 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__(operator, model, bracket, bracket_limit, density_treatment, + bracketed_method, tol, target, print_iterations, + search_for_keff_output) + + self.cell = self._get_cell(cell) + + # list of material fill the attribute cell, if depletables + self.cell_materials = None + + @abstractmethod + def _get_cell_attrib(self): + """ + Get cell attribute coefficient. + Returns + ------- + coeff : float + cell coefficient + """ + + @abstractmethod + def _set_cell_attrib(self, val): + """ + Set cell attribute to the cell instance. + Parameters + ---------- + val : float + cell coefficient to set + """ + + 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 + """ + if isinstance(val, Cell): + check_value('Cell exists', val, [cell for cell in \ + self.geometry.get_all_cells().values()]) + + elif isinstance(val, str): + if val.isnumeric(): + check_value('Cell id exists', int(val), [cell.id for cell in \ + self.geometry.get_all_cells().values()]) + val = [cell for cell in \ + self.geometry.get_all_cells().values() \ + if cell.id == int(val)][0] + else: + check_value('Cell name exists', val, [cell.name for cell in \ + self.geometry.get_all_cells().values()]) + + val = [cell for cell in \ + self.geometry.get_all_cells().values() \ + if cell.name == val][0] + + elif isinstance(val, int): + check_value('Cell id exists', val, [cell.id for cell in \ + self.geometry.get_all_cells().values()]) + val = [cell for cell in \ + self.geometry.get_all_cells().values() \ + if cell.id == val][0] + else: + ValueError(f'Cell: {val} is not supported') + + return val + def _model_builder(self, param): """ - Builds the parametric model that is passed to the `msr_search_for_keff` - function by setting the parametric variable to the geoemetrical cell. + Builds the parametric model that is passed to the `search_for_keff` + function by setting the parametric variable to the cell. Parameters ---------- param : model parametricl variable for examlple: cell translation coefficient Returns ------- - _model : openmc.model.Model + self.model : openmc.model.Model OpenMC parametric model """ self._set_cell_attrib(param) - #Calulate new volume and update if materials filling the cell are - # depletable - if self.cell_material_ids: - self._update_volumes() + # At this stage it not important to reassign the new volume, as the + # nuclide densities remain constant. However, if at least one of the cell + # materials is set as depletable, the volume change needs to be accounted for + # as an increase or reduction of number of atoms, i.e. vector x, before + # solving the depletion equations return self.model def search_for_keff(self, x, step_index): """ - Perform the criticality search on the parametric geometrical model. - Will set the root of the `search_for_keff` function to the cell - attribute. + Perform the criticality search parametrizing the cell coefficient. + The `search_for_keff` solution is then set as the new cell attribute Parameters ---------- x : list of numpy.ndarray @@ -414,22 +537,587 @@ class BatchwiseCell(Batchwise): val = self._get_cell_attrib() check_type('Cell coeff', val, Real) - # Update volume and concentration vectors before performing the search_for_keff - self._update_materials(x) + # 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 + super()._update_materials(x) # Calculate new cell attribute - root = super().search_for_keff(val) + root = super()._search_for_keff(val) # set results value as attribute in the geometry self._set_cell_attrib(root) - print('UPDATE: old value: {:.2f} cm --> ' \ - 'new value: {:.2f} cm'.format(val, root)) - # x needs to be updated after the search with new volumes if materials cell - # depletable - if self.cell_material_ids: - self._update_x() + # if at least one of the cell materials is depletable, calculate new + # volume and update x and number accordingly + # new volume + if self.cell_materials: + volumes = self._calculate_volumes() + x = super()._update_x_and_set_volumes(x, volumes) - #Store results - super()._save_res('geometry', step_index, root) - return x + return x, root + + +class BatchwiseCellGeometrical(BatchwiseCell): + """ + Batchwise 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_batchwise` method from an integrator + class, such as :class:`openmc.deplete.CECMIntegrator`. + + .. versionadded:: 0.13.4 + + Parameters + ---------- + cell : openmc.Cell or int or str + OpenMC Cell identifier to where apply batchwise scheme + attrib_name : str {'translation', 'rotation'} + Cell attribute type + operator : openmc.deplete.Operator + OpenMC operator object + model : openmc.model.Model + OpenMC model 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 + Bracketing interval to search for the solution, always 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 + Wether ot not to keep contant 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 + Wheter or not to print `search_for_keff` iterations. + Default to True + search_for_keff_output : Bool, Optional + Wheter 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 batchwise 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. + """ + def __init__(self, cell, attrib_name, operator, model, axis, 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__(cell, operator, model, bracket, bracket_limit, + 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 + + # check if cell is filled with 2 cells + check_length('fill materials', self.cell.fill.cells, 2) + #index of cell directionnal axis + check_value('axis', axis, (0,1,2)) + self.axis = axis + + # Initialize vector + self.vector = np.zeros(3) + + self.cell_materials = [cell.fill for cell in \ + self.cell.fill.cells.values() if cell.fill.depletable] + + if self.cell_materials: + self._initialize_volume_calc() + + @classmethod + def from_params(cls, obj, attr, operator, model, **kwargs): + return cls(obj, attr, operator, model, **kwargs) + + def _get_cell_attrib(self): + """ + Get cell attribute coefficient. + Returns + ------- + coeff : float + cell coefficient + """ + for cell in openmc.lib.cells.values(): + if cell.id == self.cell.id: + if self.attrib_name == 'translation': + return cell.translation[self.axis] + elif self.attrib_name == 'rotation': + return cell.rotation[self.axis] + + def _set_cell_attrib(self, val): + """ + Set cell attribute to the cell instance. + Attributes are only applied to cells filled with a universe + Parameters + ---------- + val : float + Cell coefficient to set in cm for translation and def for rotation + """ + self.vector[self.axis] = val + for cell in openmc.lib.cells.values(): + if cell.id == self.cell.id: + setattr(cell, self.attrib_name, self.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.cell_materials, 100000, ll, ur) + self.model.settings.volume_calculations = mat_vol + + def _calculate_volumes(self): + """ + Perform stochastic volume calculation + Returns + ------- + volumes : dict + Dictionary of calculate volumes, where key is the mat id + """ + openmc.lib.calculate_volumes() + volumes = {} + res = openmc.VolumeCalculation.from_hdf5('volume_1.h5') + for mat_idx, mat in enumerate(self.local_mats): + if int(mat) in [mat.id for mat in self.cell_materials]: + volumes[mat] = res.volumes[int(mat)].n + return volumes + +class BatchwiseCellTemperature(BatchwiseCell): + """ + Batchwise 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_batchwise` method from an integrator + class, such as :class:`openmc.deplete.CECMIntegrator`. + + .. versionadded:: 0.13.4 + + Parameters + ---------- + cell : openmc.Cell or int or str + OpenMC Cell identifier to where apply batchwise scheme + operator : openmc.deplete.Operator + OpenMC operator object + model : openmc.model.Model + OpenMC model object + bracket : list of float + Bracketing interval to search for the solution, always 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 + Wether ot not to keep contant 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 + Wheter or not to print `search_for_keff` iterations. + Default to True + search_for_keff_output : Bool, Optional + Wheter 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 batchwise scheme + """ + def __init__(self, cell, operator, model, 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__(cell, operator, model, bracket, bracket_limit, + density_treatment, bracketed_method, tol, target, + print_iterations, search_for_keff_output) + + # 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) + + @classmethod + def from_params(cls, obj, attr, operator, model, **kwargs): + return cls(obj, operator, model, **kwargs) + + def _get_cell_attrib(self): + """ + Get cell attribute coefficient. + Returns + ------- + coeff : float + cell temperature in Kelvin + """ + for cell in openmc.lib.cells.values(): + if cell.id == self.cell.id: + return cell.get_temperature() + + def _set_cell_attrib(self, val): + """ + Set cell attribute to the cell instance. + Parameters + ---------- + val : float + Cell temperature to set in Kelvin + """ + for cell in openmc.lib.cells.values(): + if cell.id == self.cell.id: + cell.set_temperature(val) + +class BatchwiseMaterial(Batchwise): + """Abstract class holding batchwise material-based functions. + + Specific classes for running batchwise depletion calculations are + implemented as derived class of BatchwiseMaterial. + + .. versionadded:: 0.13.4 + + Parameters + ---------- + material : openmc.Material or int or str + OpenMC Material identifier to where apply batchwise scheme + operator : openmc.deplete.Operator + OpenMC operator object + model : openmc.model.Model + OpenMC model object + mat_vector : dict + Dictionary of material composition to paremetrize, where a pair key, value + represents a nuclide and its weight fraction. + bracket : list of float + Bracketing interval to search for the solution, always 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 + Wether ot not to keep contant 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 + Wheter or not to print `search_for_keff` iterations. + Default to True + search_for_keff_output : Bool, Optional + Wheter 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 batchwise scheme + mat_vector : dict + Dictionary of material composition to paremetrize, where a pair key, value + represents a nuclide and its weight fraction. + """ + + def __init__(self, material, operator, model, 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__(operator, model, bracket, bracket_limit, + density_treatment, bracketed_method, tol, target, + print_iterations, search_for_keff_output) + + self.material = self._get_material(material) + + check_type("material vector", mat_vector, dict, str) + for nuc in mat_vector.keys(): + check_value("check nuclide exists", nuc, self.operator.nuclides_with_data) + + if round(sum(mat_vector.values()), 2) != 1.0: + # 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 + + @classmethod + def from_params(cls, obj, attr, operator, model, **kwargs): + return cls(obj, operator, model, **kwargs) + + 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/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 + + @abstractmethod + def _model_builder(self, param): + """ + Builds the parametric model that is passed to the `msr_search_for_keff` + function by updating the material densities and setting the parametric + variable as function of the nuclides vector. Since this is a paramteric + material addition (or removal), we can parametrize the volume as well. + Parameters + ---------- + param : + Model material function variable + Returns + ------- + _model : openmc.model.Model + Openmc parametric model + """ + + def search_for_keff(self, x, step_index): + """ + Perform the criticality search on the parametric material model. + Will set the root of the `search_for_keff` function to the atoms + concentrations vector. + Parameters + ---------- + x : list of numpy.ndarray + Total atoms concentrations + Returns + ------- + x : list of numpy.ndarray + Updated total atoms concentrations + """ + # Update AtomNumber with new conc vectors x. Materials are also updated + # even though they are also re-calculated when running the search_for_kef + super()._update_materials(x) + + # Solve search_for_keff and find new value + root = super()._search_for_keff(0) + + #Update concentration vector and volumes with new value + volumes = self._calculate_volumes(root) + x = super()._update_x_and_set_volumes(x, volumes) + + return x, root + +class BatchwiseMaterialRefuel(BatchwiseMaterial): + """ + Batchwise material-based class for refueling (addition or removal) scheme. + + 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_batchwise` method from an integrator + class, such as :class:`openmc.deplete.CECMIntegrator`. + + .. versionadded:: 0.13.4 + + Parameters + ---------- + material : openmc.Material or int or str + OpenMC Material identifier to where apply batchwise scheme + operator : openmc.deplete.Operator + OpenMC operator object + model : openmc.model.Model + OpenMC model object + mat_vector : dict + Dictionary of material composition to paremetrize, where a pair key, value + represents a nuclide and its weight fraction. + bracket : list of float + Bracketing interval to search for the solution, always 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 + Wether ot not to keep contant 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 + Wheter or not to print `search_for_keff` iterations. + Default to True + search_for_keff_output : Bool, Optional + Wheter 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 batchwise scheme + mat_vector : dict + Dictionary of material composition to paremetrize, where a pair key, value + represents a nuclide and its weight fraction. + """ + + def __init__(self, material, operator, model, 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, model, mat_vector, bracket, + bracket_limit, density_treatment, bracketed_method, tol, + target, print_iterations, search_for_keff_output) + + @classmethod + def from_params(cls, obj, attr, operator, model, **kwargs): + return cls(obj, operator, model, **kwargs) + + def _model_builder(self, param): + """ + Builds the parametric model that is passed to the `msr_search_for_keff` + function by updating the material densities and setting the parametric + variable to the material nuclides to add. We can either fix the total + volume or the material density according to user input. + Parameters + ---------- + param : + Model function variable, total grams of material to add or remove + 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 in number_i.materials: + nuclides = [] + densities = [] + + if int(mat) == self.material.id: + + if self.density_treatment == 'constant-density': + vol = number_i.get_mat_volume(mat) + \ + (param / self.material.get_mass_density()) + + elif self.density_treatment == 'constant-volume': + vol = number_i.get_mat_volume(mat) + + 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) + + else: + # get normalized atoms density in [atoms/b-cm] + val = 1.0e-24 * number_i[mat, nuc] / vol + + if val > 0.0: + nuclides.append(nuc) + densities.append(val) + + 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) + + if val > 0.0: + nuclides.append(nuc) + densities.append(val) + + #set nuclides and densities to the in-memory model + openmc.lib.materials[int(mat)].set_densities(nuclides, densities) + + # alwyas need to return a model + return self.model + + def _calculate_volumes(self, res): + """ + """ + 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) + return volumes diff --git a/openmc/deplete/stepresult.py b/openmc/deplete/stepresult.py index 372f7cf25d..3aa8b45dc0 100644 --- a/openmc/deplete/stepresult.py +++ b/openmc/deplete/stepresult.py @@ -74,7 +74,7 @@ class StepResult: self.mat_to_hdf5_ind = None self.data = None - + self.batchwise = None def __repr__(self): t = self.time[0] dt = self.time[1] - self.time[0] @@ -349,6 +349,10 @@ class StepResult: "depletion time", (1,), maxshape=(None,), dtype="float64") + handle.create_dataset( + "batchwise_root", (1,), maxshape=(None,), + dtype="float64") + def _to_hdf5(self, handle, index, parallel=False): """Converts results object into an hdf5 object. @@ -379,6 +383,7 @@ class StepResult: time_dset = handle["/time"] source_rate_dset = handle["/source_rate"] proc_time_dset = handle["/depletion time"] + root_dset = handle["/batchwise_root"] # Get number of results stored number_shape = list(number_dset.shape) @@ -412,6 +417,10 @@ class StepResult: 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 +444,7 @@ class StepResult: proc_time_dset[index] = ( self.proc_time / (comm.size * self.n_hdf5_mats) ) + root_dset[index] = self.batchwise @classmethod def from_hdf5(cls, handle, step): @@ -458,11 +468,13 @@ class StepResult: else: # Older versions used "power" instead of "source_rate" source_rate_dset = handle["/power"] + root_dset = handle["/batchwise_root"] results.data = number_dset[step, :, :, :] results.k = eigenvalues_dset[step, :] results.time = time_dset[step, :] results.source_rate = source_rate_dset[step, 0] + results.batchwise = root_dset[step] if "depletion time" in handle: proc_time_dset = handle["/depletion time"] @@ -509,7 +521,7 @@ class StepResult: @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 @@ -564,6 +576,7 @@ class StepResult: results.proc_time = proc_time if results.proc_time is not None: results.proc_time = comm.reduce(proc_time, op=MPI.SUM) + results.batchwise = 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 7f6b3b5b9a..4818fbe10b 100644 --- a/openmc/search.py +++ b/openmc/search.py @@ -12,7 +12,7 @@ _SCALAR_BRACKETED_METHODS = ['brentq', 'brenth', 'ridder', 'bisect'] def _search_keff(guess, target, model_builder, model_args, print_iterations, - run_args, guesses, results): + run_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 @@ -51,7 +51,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 +74,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. @@ -193,12 +197,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 - 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: + print(f'WARNING: {root_res.flag}') + return guesses, results + # In case the root finder is not successful + except Exception as e: + print(f'WARNING: {e}') + return guesses, results diff --git a/tests/unit_tests/test_deplete_batchwise.py b/tests/unit_tests/test_deplete_batchwise.py new file mode 100644 index 0000000000..7c05786fa5 --- /dev/null +++ b/tests/unit_tests/test_deplete_batchwise.py @@ -0,0 +1,135 @@ +""" Tests for Batchwise class """ + +from pathlib import Path +from math import exp + +import pytest +import numpy as np + +import openmc +import openmc.lib +from openmc.deplete import CoupledOperator +from openmc.deplete import (BatchwiseCellGeometrical, BatchwiseCellTemperature, + BatchwiseMaterialRefuel) + +CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml" + +@pytest.fixture +def model(): + f = openmc.Material(name="f") + 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="w") + w.add_element("O", 1) + w.add_element("H", 2) + w.set_density("g/cc", 1.0) + w.depletable = True + h = openmc.Material(name='h') + h.set_density('g/cm3', 0.001598) + h.add_element('He', 2.4044e-4) + radii = [0.42, 0.45] + f.volume = np.pi * radii[0] ** 2 + w.volume = np.pi * (radii[1]**2 - radii[0]**2) + materials = openmc.Materials([f, w, h]) + surf_wh = openmc.ZPlane(z0=0) + surf_f = openmc.Sphere(r=radii[0]) + surf_w = openmc.Sphere(r=radii[1], boundary_type='vacuum') + cell_w = openmc.Cell(fill=w, region=-surf_wh) + cell_h = openmc.Cell(fill=h, region=+surf_wh) + uni_wh = openmc.Universe(cells=(cell_w, cell_h)) + cell_f = openmc.Cell(name='f',fill=f, region=-surf_f) + cell_wh = openmc.Cell(name='wh',fill=uni_wh, region=+surf_f & -surf_w) + geometry = openmc.Geometry([cell_f, cell_wh]) + settings = openmc.Settings() + settings.particles = 1000 + settings.inactive = 10 + settings.batches = 50 + return openmc.Model(geometry, materials, settings) + +@pytest.mark.parametrize("object_name, attribute", [ + ('wh', 'translation'), + ('wh', 'rotation'), + ('f', 'temperature'), + ('f', 'refuel'), + ]) +def test_attributes(model, object_name, attribute): + """ + Test classes attributes are set correctly + """ + op = CoupledOperator(model, CHAIN_PATH) + + bracket = [-1,1] + bracket_limit = [-10,10] + axis = 2 + mat_vector = {'U235':0.1, 'U238':0.9} + if attribute in ('translation','rotation'): + bw = BatchwiseCellGeometrical(object_name, attribute, op, model, axis, bracket, bracket_limit) + assert bw.cell_materials == [cell.fill for cell in \ + model.geometry.get_cells_by_name(object_name)[0].fill.cells.values() \ + if cell.fill.depletable] + assert bw.attrib_name == attribute + assert bw.axis == axis + + elif attribute == 'temperature': + bw = BatchwiseCellTemperature(object_name, op, model, bracket, bracket_limit) + elif attribute == 'refuel': + bw = BatchwiseMaterialRefuel(object_name, op, model, mat_vector, bracket, bracket_limit) + assert bw.mat_vector == mat_vector + + assert bw.bracket == bracket + assert bw.bracket_limit == bracket_limit + assert bw.burn_mats == op.burnable_mats + assert bw.local_mats == op.local_mats + +@pytest.mark.parametrize("attribute", [ + ('translation'), + ('rotation') + ]) +def test_cell_methods(run_in_tmpdir, model, attribute): + """ + Test cell base class internal method + """ + bracket = [-1,1] + bracket_limit = [-10,10] + axis = 2 + op = CoupledOperator(model, CHAIN_PATH) + integrator = openmc.deplete.PredictorIntegrator( + op, [1,1], 0.0, timestep_units = 'd') + integrator.add_batchwise('wh', attribute, axis=axis, bracket=bracket, + bracket_limit=bracket_limit) + + model.export_to_xml() + openmc.lib.init() + integrator.batchwise._set_cell_attrib(0) + assert integrator.batchwise._get_cell_attrib() == 0 + vol = integrator.batchwise._calculate_volumes() + for cell in integrator.batchwise.cell_materials: + assert vol[str(cell.id)] == pytest.approx([mat.volume for mat in model.materials if mat.id == cell.id]) + openmc.lib.finalize() + +def test_update_materials_methods(run_in_tmpdir): + """ + It's the same methods for all classe so one class is enough + """ + bracket = [-1,1] + bracket_limit = [-10,10] + axis = 2 + op = CoupledOperator(model, CHAIN_PATH) + integrator = openmc.deplete.PredictorIntegrator( + op, [1,1], 0.0, timestep_units = 'd') + integrator.add_batchwise('wh', 'translation', axis=axis, bracket=bracket, + bracket_limit=bracket_limit) + + model.export_to_xml() + openmc.lib.init() + x = op.number.number + #Double number of atoms of U238 in fuel + mat_index = op.number.index_mat['1'] + nuc_index = op.number.index_nuc['U238'] + vol = op.number.get_mat_volume('1') + op.number.number[mat_index][nuc_idx] *= 2 + integrator.batchwise._update_volumes(op.number.number) + + vol = integrator.batchwise._calculate_volumes()