diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index fbf43a808e..3055561be7 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -581,3 +581,42 @@ class SILEQIIntegrator(SIIntegrator): proc_time += time1 + time2 return proc_time, [eos_conc, inter_conc], [res_bar] + + +def integrator_factory(method): + """This method is a factor for the integrator sub-classes + + Params + ------ + method : str + The type of integrator method to use. Valid values are: 'cecm', + 'predictor', 'cf4', 'epc_rk4', 'si_celi', 'si_leqi', 'celi', and 'leqi' + + Returns + ------- + integrator : Integrator + The type of integrator + + """ + + if method == 'cecm': + integrator = CECMIntegrator + elif method == 'predictor': + integrator = PredictorIntegrator + elif method == 'cf4': + integrator = CF4Integrator + elif method == 'epc_rk4': + integrator = EPCRK4Integrator + elif method == 'si_celi': + integrator = SICELIIntegrator + elif method == 'si_leqi': + integrator = SILEQIIntegrator + elif method == 'celi': + integrator = CELIIntegrator + elif method == 'leqi': + integrator = LEQIIntegrator + else: + msg = "Invalid integrator method: {}!".format(method) + raise ValueError(msg) + + return integrator diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 0e5cc9efa3..9c67d3bec0 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -536,7 +536,8 @@ class Operator(TransportOperator): # Initialize OpenMC library comm.barrier() - openmc.lib.init(intracomm=comm) + if not openmc.lib.LIB_INIT: + openmc.lib.init(intracomm=comm) # Generate tallies in memory materials = [openmc.lib.materials[int(i)] diff --git a/openmc/lib/__init__.py b/openmc/lib/__init__.py index 82dc92ba48..7e766a93b2 100644 --- a/openmc/lib/__init__.py +++ b/openmc/lib/__init__.py @@ -59,3 +59,5 @@ from .tally import * from .settings import settings from .math import * from .plot import * + +LIB_INIT = False diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 0da24aab88..e9601c3320 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -109,6 +109,7 @@ def global_bounding_box(): return llc, urc + def calculate_volumes(): """Run stochastic volume calculation""" _dll.openmc_calculate_volumes() @@ -147,6 +148,7 @@ def export_properties(filename=None): def finalize(): """Finalize simulation and free memory""" _dll.openmc_finalize() + openmc.lib.LIB_INIT = False def find_cell(xyz): @@ -253,6 +255,7 @@ def init(args=None, intracomm=None): intracomm = c_void_p(address) _dll.openmc_init(argc, argv, intracomm) + openmc.lib.LIB_INIT = True def is_statepoint_batch(): diff --git a/openmc/lib/material.py b/openmc/lib/material.py index 6770721ab6..fde197d3d8 100644 --- a/openmc/lib/material.py +++ b/openmc/lib/material.py @@ -172,7 +172,6 @@ class Material(_FortranObjectWithID): @property def nuclides(self): return self._get_densities()[0] - return nuclides @property def densities(self): diff --git a/openmc/model/model.py b/openmc/model/model.py index d6c91e73e1..7f58e9589c 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1,11 +1,17 @@ from collections.abc import Iterable from pathlib import Path import time +import warnings import h5py +import numpy as np import openmc -from openmc.checkvalue import check_type, check_value +from openmc.checkvalue import check_type, check_value, check_iterable_type, \ + check_length +import openmc.deplete as dep +from openmc.data.library import DataLibrary +from openmc.exceptions import DataError, InvalidIDError, SetupError class Model: @@ -31,6 +37,14 @@ class Model: Tallies information plots : openmc.Plots, optional Plot information + chain_file : str or Path, optional + Path to the depletion chain XML file. Defaults to the chain + found under the ``depletion_chain`` in the + :envvar:`OPENMC_CROSS_SECTIONS` environment variable if it exists. If a + str is provided it will be converted to a Path object. + fission_q : dict, optional + Dictionary of nuclides and their fission Q values [eV]. + If not given, values will be pulled from the ``chain_file``. Attributes ---------- @@ -44,11 +58,19 @@ class Model: Tallies information plots : openmc.Plots Plot information + chain_file : str or Path + Path to the depletion chain XML file. Defaults to the chain + found under the ``depletion_chain`` in the + :envvar:`OPENMC_CROSS_SECTIONS` environment variable if it exists. If a + str is provided it will be converted to a Path object. + fission_q : dict + Dictionary of nuclides and their fission Q values [eV]. + If not given, values will be pulled from the ``chain_file``. """ def __init__(self, geometry=None, materials=None, settings=None, - tallies=None, plots=None): + tallies=None, plots=None, chain_file=None, fission_q=None): self.geometry = openmc.Geometry() self.materials = openmc.Materials() self.settings = openmc.Settings() @@ -66,6 +88,19 @@ class Model: if plots is not None: self.plots = plots + self.chain_file = chain_file + self.fission_q = fission_q + self.depletion_operator = None + + if self.materials is None: + mats = self.geometry.get_all_materials().values() + else: + mats = self.materials + self._materials_by_id = {mat.id: mat for mat in mats} + cells = self.geometry.get_all_cells() + self._cells_by_id = {cell.id: cell for cell in cells.values()} + self._cells_by_name = {cell.name: cell for cell in cells.values()} + @property def geometry(self): return self._geometry @@ -86,6 +121,22 @@ class Model: def plots(self): return self._plots + @property + def chain_file(self): + return self._chain_file + + @property + def fission_q(self): + return self._fission_q + + @property + def depletion_operator(self): + return self._depletion_operator + + @property + def C_init(self): + return openmc.lib.LIB_INIT + @geometry.setter def geometry(self, geometry): check_type('geometry', geometry, openmc.Geometry) @@ -126,6 +177,25 @@ class Model: for plot in plots: self._plots.append(plot) + @chain_file.setter + def chain_file(self, chain_file): + check_type('chain_file', chain_file, (type(None), str, Path)) + if isinstance(chain_file, str): + self._chain_file = Path(chain_file) + else: + self._chain_file = chain_file + + @fission_q.setter + def fission_q(self, fission_q): + check_type('fission_q', fission_q, (type(None), dict)) + self._fission_q = fission_q + + @depletion_operator.setter + def depletion_operator(self, depletion_operator): + check_type('depletion_operator', depletion_operator, + (type(None), dep.Operator)) + self._depletion_operator = depletion_operator + @classmethod def from_xml(cls, geometry='geometry.xml', materials='materials.xml', settings='settings.xml'): @@ -151,8 +221,39 @@ class Model: settings = openmc.Settings.from_xml(settings) return cls(geometry, materials, settings) + def init_C_api(self, use_depletion_operator=False): + """Initializes the model in memory via the C-API + + Parameters + ---------- + directory : str + Directory to write XML files to. If it doesn't exist already, it + will be created. + use_depletion_operator : bool, optional + If True, the model will be loaded using the depletion operator + including all isotopes necessary from fission. This parameter will + use the :attr:`Model.chain_file` and :attr:`Model.fission_q` + attributes. + """ + + if use_depletion_operator: + # Create OpenMC transport operator + self.depletion_operator = \ + dep.Operator(self.geometry, self.settings, + str(self.chain_file), fission_q=self.fission_q) + else: + openmc.lib.hard_reset() + if dep.comm.rank == 0: + self.export_to_xml() + dep.comm.barrier() + openmc.lib.init(intracomm=dep.comm) + + def clear_C_api(self): + """Finalize simulation and free memory allocated for the C-API""" + openmc.lib.finalize() + def deplete(self, timesteps, chain_file=None, method='cecm', - fission_q=None, **kwargs): + fission_q=None, final_step=True, **kwargs): """Deplete model using specified timesteps/power Parameters @@ -169,26 +270,75 @@ class Model: fission_q : dict, optional Dictionary of nuclides and their fission Q values [eV]. If not given, values will be pulled from the ``chain_file``. + final_step : bool, optional + Indicate whether or not a transport solve should be run at the end + of the last timestep. + + .. versionadded:: 0.12.3 **kwargs Keyword arguments passed to integration function (e.g., :func:`openmc.deplete.integrator.cecm`) """ - # Import the depletion module. This is done here rather than the module - # header to delay importing openmc.lib (through openmc.deplete) which - # can be tough to install properly. - import openmc.deplete as dep - # Create OpenMC transport operator - op = dep.Operator( - self.geometry, self.settings, chain_file, - fission_q=fission_q, - ) + if self.C_init and self.depletion_operator is not None: + # Then the user has properly initialized the information and we can + # just carry forward + pass + elif self.C_init and self.depletion_operator is None: + # Then the user has initialzed the C-API but without the depletion + # isotopes loaded. We would have to reset and reload data, but + # doing so could lose user information. Therefore let us just + # provide an error and quit. + msg = "Model.deplete(...) cannot be called after " \ + "Model.init_C_api(...) if the use_depletion_operator " \ + "argument to Model.init_C_api(...) is False." + raise SetupError(msg) + else: + # To get here, the C-API is not initialized. So we can do that now + # To keep Model.deplete(...) API compatibility, we will allow the + # chain_file and fission_q params to be set since we havent loaded + # the API anyways + if chain_file is not None: + self.chain_file = chain_file + warnings.warn("The chain_file argument of Model.deplete(...) " + "has been deprecated and may be removed in a " + "future version. The Model.chain_file should be" + "used instead.", DeprecationWarning) + if fission_q is not None: + warnings.warn("The fission_q argument of Model.deplete(...) " + "has been deprecated and may be removed in a " + "future version. The Model.fission_q should be" + "used instead.", DeprecationWarning) + self.fission_q = fission_q + self.init_C_api(use_depletion_operator=True) - # Perform depletion - check_value('method', method, ('cecm', 'predictor', 'cf4', 'epc_rk4', - 'si_celi', 'si_leqi', 'celi', 'leqi')) - getattr(dep.integrator, method)(op, timesteps, **kwargs) + # Set up the integrator + integrator_class = dep.integrators.integrator_factory(method) + integrator = integrator_class(self.depletion_operator, + timesteps, **kwargs) + + # Now perform the depletion + integrator.integrate(final_step) + + # If we did not perform a transport calculation on the final step, then + # make the code update the C-API material inventory + if not final_step: + self.depletion_operator._update_materials() + + # Now make the python Materials match the C-API material data + for mat_id, mat in self._materials_by_id.items(): + if mat.depletable: + # Get the C data + c_mat = openmc.lib.materials[mat_id] + nuclides, densities = c_mat._get_densities() + # And now we can remove isotopes and add these ones in + atom_density = 0. + for nuc, density in zip(nuclides, densities): + mat.remove_nuclide(nuc) # Replace if it's there + mat.add_nuclide(nuc, density) + atom_density += density + mat.set_density('atom/b-cm', atom_density) def export_to_xml(self, directory='.'): """Export model to XML files. @@ -269,13 +419,18 @@ class Model: materials[mat_id].set_density('atom/b-cm', atom_density) def run(self, **kwargs): - """Creates the XML files, runs OpenMC, and returns the path to the last + """Runs OpenMC. If the C-API has been initialized, then the C-API is + used, otherwise, this method creates the XML files and runs OpenMC via + a system cal. In both cases this method returns the path to the last statepoint file generated. .. versionchanged:: 0.12 Instead of returning the final k-effective value, this function now returns the path to the final statepoint written. + .. versionchanged:: 0.12.3 + This method can utilize the C-API for execution + Parameters ---------- **kwargs @@ -289,16 +444,20 @@ class Model: """ - self.export_to_xml() - - # Setting tstart here ensures we don't pick up any pre-existing statepoint - # files in the output directory + # Setting tstart here ensures we don't pick up any pre-existing + # statepoint files in the output directory tstart = time.time() last_statepoint = None - openmc.run(**kwargs) + if self.C_init: + # Then run using the C-API + openmc.lib.run() + else: + # Then run via the command line + self.export_to_xml() + openmc.run(**kwargs) - # Get output directory and return the last statepoint written by this run + # Get output directory and return the last statepoint written this run if self.settings.output and 'path' in self.settings.output: output_dir = Path(self.settings.output['path']) else: @@ -309,3 +468,78 @@ class Model: tstart = mtime last_statepoint = sp return last_statepoint + + def _move_cell(self, cell_names_or_ids, vector, attrib_name): + # Method to do the same work whether it is a rotation or translation + check_type('cell_names_or_ids', cell_names_or_ids, Iterable, + (np.int, int, str)) + check_type('vector', vector, Iterable, (np.float, float)) + check_length('vector', vector, 3) + check_value('attrib_name', attrib_name, ('rotation', 'translation')) + + # Get the list of cell ids to use y converting from names and accepting + # only values that have actual ids + cell_ids = [None] * len(cell_names_or_ids) + for c, cell_name_or_id in enumerate(cell_names_or_ids): + if isinstance(cell_name_or_id, (int, np.int)): + if cell_name_or_id in self._cells_by_id: + cell_ids[c] = int(cell_name_or_id) + msg = 'Cell ID {} is not present in the model!'.format( + cell_name_or_id) + raise InvalidIDError(msg) + elif isinstance(cell_name_or_id, str): + if cell_name_or_id in self._cells_by_name: + cell_ids[c] = self._cells_by_name[cell_name_or_id] + else: + msg = 'Cell {} is not present in the model!'.format( + cell_name_or_id) + raise InvalidIDError(msg) + + # Now perform the motion + for cell_id in cell_ids: + cell = self._cells_by_id[cell_id] + if attrib_name == 'rotation': + cell.rotation = vector + elif attrib_name == 'translation': + cell.translation = vector + # Next lets keep what is in C-API memory up to date as well + if self.C_init: + C_cell = openmc.lib.cells[cell_id] + if attrib_name == 'rotation': + C_cell.rotation = vector + elif attrib_name == 'translation': + C_cell.translation = vector + + def rotate_cells(self, cell_names_or_ids, vector): + """Rotate the identified cell(s) by the specified rotation vector. + The rotation is only applied to cells filled with a universe. + + Parameters + ---------- + cell_names_or_ids : Iterable of str or int + The cell names (if str) or id (if int) that are to be translated + or rotated. This parameter can include a mix of names and ids. + vector : Iterable of float + The rotation vector of length 3 to apply. This array specifies the + angles in degrees about the x, y, and z axes, respectively. + + """ + + self._move_cell(cell_names_or_ids, vector, 'rotation') + + def translate_cells(self, cell_names_or_ids, vector): + """Translate the identified cell(s) by the specified translation vector. + The translation is only applied to cells filled with a universe. + + Parameters + ---------- + cell_names_or_ids : Iterable of str or int + The cell names (if str) or id (if int) that are to be translated + or rotated. This parameter can include a mix of names and ids. + vector : Iterable of float + The translation vector of length 3 to apply. This array specifies + the x, y, and z dimensions of the translation. + + """ + + self._move_cell(cell_names_or_ids, vector, 'translation')