diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index cf62debc0..0c870c9fb 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -126,6 +126,14 @@ std::string distribcell_path( int maximum_levels(int32_t univ); +//============================================================================== +//! Check whether or not a universe is the root universe using its ID. +//! \param univ_id The ID of the universe to check. +//! \return Whether or not it is the root universe. +//============================================================================== + +bool is_root_universe(int32_t univ_id); + //============================================================================== //! Deallocates global vectors and maps for cells, universes, and lattices. //============================================================================== diff --git a/man/man1/openmc.1 b/man/man1/openmc.1 index 4aa288a8c..331e42ea0 100644 --- a/man/man1/openmc.1 +++ b/man/man1/openmc.1 @@ -8,7 +8,10 @@ sections are available. .SH SYNOPSIS \fBopenmc\fR [\fIoptions\fR] [\fIpath\fR] .PP -It is assumed that if no +.I path +specifies either the path to a single model XML file containing the full model +or a directory containing either a model.xml file or a set of individual XML +files (settings.xml, materials.xml, geometry.xml). It is assumed that if no .I path is specified, the XML input files are present in the current directory. .SH OPTIONS @@ -46,6 +49,9 @@ The behavior of .B openmc is affected by the following environment variables. .TP +.B OPENMC_CHAIN_FILE +Indicates the path to a depletion chain XML file. +.TP .B OPENMC_CROSS_SECTIONS Indicates the default path to the cross_sections.xml summary file that is used to locate HDF5 format cross section libraries if the user has not specified the diff --git a/openmc/data/data.py b/openmc/data/data.py index 6793dff2e..16c7ecf54 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -207,8 +207,8 @@ _LOG_TWO = log(2.0) def atomic_mass(isotope): """Return atomic mass of isotope in atomic mass units. - Atomic mass data comes from the `Atomic Mass Evaluation 2016 - `_. + Atomic mass data comes from the `Atomic Mass Evaluation 2020 + `_. Parameters ---------- diff --git a/openmc/geometry.py b/openmc/geometry.py index ad86cc56b..ee4a58a3b 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -1,11 +1,14 @@ +from __future__ import annotations import os import typing from collections import OrderedDict, defaultdict -from collections.abc import Iterable from copy import deepcopy +from collections.abc import Iterable from pathlib import Path -from xml.etree import ElementTree as ET import warnings +from xml.etree import ElementTree as ET + +import numpy as np import openmc import openmc._xml as xml @@ -51,19 +54,19 @@ class Geometry: self._root_universe = univ @property - def root_universe(self): + def root_universe(self) -> openmc.UniverseBase: return self._root_universe @property - def bounding_box(self): + def bounding_box(self) -> np.ndarray: return self.root_universe.bounding_box @property - def merge_surfaces(self): + def merge_surfaces(self) -> bool: return self._merge_surfaces @property - def surface_precision(self): + def surface_precision(self) -> int: return self._surface_precision @root_universe.setter @@ -105,7 +108,7 @@ class Geometry: if universe.id in volume_calc.volumes: universe.add_volume_information(volume_calc) - def to_xml_element(self, remove_surfs=False): + def to_xml_element(self, remove_surfs=False) -> ET.Element: """Creates a 'geometry' element to be written to an XML file. Parameters @@ -164,7 +167,7 @@ class Geometry: tree.write(str(p), xml_declaration=True, encoding='utf-8') @classmethod - def from_xml_element(cls, elem, materials=None): + def from_xml_element(cls, elem, materials=None) -> Geometry: """Generate geometry from an XML element Parameters @@ -260,7 +263,7 @@ class Geometry: cls, path: PathLike = 'geometry.xml', materials: typing.Optional[typing.Union[PathLike, 'openmc.Materials']] = 'materials.xml' - ): + ) -> Geometry: """Generate geometry from XML file Parameters @@ -290,7 +293,7 @@ class Geometry: return cls.from_xml_element(root, materials) - def find(self, point): + def find(self, point) -> list: """Find cells/universes/lattices which contain a given point Parameters @@ -307,7 +310,7 @@ class Geometry: """ return self.root_universe.find(point) - def get_instances(self, paths): + def get_instances(self, paths) -> typing.Union[int, typing.List[int]]: """Return the instance number(s) for a cell/material in a geometry path. The instance numbers are used as indices into distributed @@ -354,7 +357,7 @@ class Geometry: return indices if return_list else indices[0] - def get_all_cells(self): + def get_all_cells(self) -> typing.OrderedDict[int, openmc.Cell]: """Return all cells in the geometry. Returns @@ -368,7 +371,7 @@ class Geometry: else: return OrderedDict() - def get_all_universes(self): + def get_all_universes(self) -> typing.OrderedDict[int, openmc.Universe]: """Return all universes in the geometry. Returns @@ -383,7 +386,7 @@ class Geometry: universes.update(self.root_universe.get_all_universes()) return universes - def get_all_materials(self): + def get_all_materials(self) -> typing.OrderedDict[int, openmc.Material]: """Return all materials within the geometry. Returns @@ -398,7 +401,7 @@ class Geometry: else: return OrderedDict() - def get_all_material_cells(self): + def get_all_material_cells(self) -> typing.OrderedDict[int, openmc.Cell]: """Return all cells filled by a material Returns @@ -417,7 +420,7 @@ class Geometry: return material_cells - def get_all_material_universes(self): + def get_all_material_universes(self) -> typing.OrderedDict[int, openmc.Universe]: """Return all universes having at least one material-filled cell. This method can be used to find universes that have at least one cell @@ -440,7 +443,7 @@ class Geometry: return material_universes - def get_all_lattices(self): + def get_all_lattices(self) -> typing.OrderedDict[int, openmc.Lattice]: """Return all lattices defined Returns @@ -458,7 +461,7 @@ class Geometry: return lattices - def get_all_surfaces(self): + def get_all_surfaces(self) -> typing.OrderedDict[int, openmc.Surface]: """ Return all surfaces used in the geometry @@ -475,7 +478,7 @@ class Geometry: surfaces = cell.region.get_surfaces(surfaces) return surfaces - def _get_domains_by_name(self, name, case_sensitive, matching, domain_type): + def _get_domains_by_name(self, name, case_sensitive, matching, domain_type) -> list: if not case_sensitive: name = name.lower() @@ -492,7 +495,9 @@ class Geometry: domains.sort(key=lambda x: x.id) return domains - def get_materials_by_name(self, name, case_sensitive=False, matching=False): + def get_materials_by_name( + self, name, case_sensitive=False, matching=False + ) -> typing.List[openmc.Material]: """Return a list of materials with matching names. Parameters @@ -513,7 +518,9 @@ class Geometry: """ return self._get_domains_by_name(name, case_sensitive, matching, 'material') - def get_cells_by_name(self, name, case_sensitive=False, matching=False): + def get_cells_by_name( + self, name, case_sensitive=False, matching=False + ) -> typing.List[openmc.Cell]: """Return a list of cells with matching names. Parameters @@ -534,7 +541,9 @@ class Geometry: """ return self._get_domains_by_name(name, case_sensitive, matching, 'cell') - def get_surfaces_by_name(self, name, case_sensitive=False, matching=False): + def get_surfaces_by_name( + self, name, case_sensitive=False, matching=False + ) -> typing.List[openmc.Surface]: """Return a list of surfaces with matching names. .. versionadded:: 0.13.3 @@ -557,7 +566,9 @@ class Geometry: """ return self._get_domains_by_name(name, case_sensitive, matching, 'surface') - def get_cells_by_fill_name(self, name, case_sensitive=False, matching=False): + def get_cells_by_fill_name( + self, name, case_sensitive=False, matching=False + ) -> typing.List[openmc.Cell]: """Return a list of cells with fills with matching names. Parameters @@ -602,7 +613,9 @@ class Geometry: return sorted(cells, key=lambda x: x.id) - def get_universes_by_name(self, name, case_sensitive=False, matching=False): + def get_universes_by_name( + self, name, case_sensitive=False, matching=False + ) -> typing.List[openmc.Universe]: """Return a list of universes with matching names. Parameters @@ -623,7 +636,9 @@ class Geometry: """ return self._get_domains_by_name(name, case_sensitive, matching, 'universe') - def get_lattices_by_name(self, name, case_sensitive=False, matching=False): + def get_lattices_by_name( + self, name, case_sensitive=False, matching=False + ) -> typing.List[openmc.Lattice]: """Return a list of lattices with matching names. Parameters @@ -644,7 +659,7 @@ class Geometry: """ return self._get_domains_by_name(name, case_sensitive, matching, 'lattice') - def remove_redundant_surfaces(self): + def remove_redundant_surfaces(self) -> typing.Dict[int, openmc.Surface]: """Remove and return all of the redundant surfaces. Uses surface_precision attribute of Geometry instance for rounding and @@ -707,7 +722,7 @@ class Geometry: # Recursively traverse the CSG tree to count all cell instances self.root_universe._determine_paths(instances_only=instances_only) - def clone(self): + def clone(self) -> Geometry: """Create a copy of this geometry with new unique IDs for all of its enclosed materials, surfaces, cells, universes and lattices.""" diff --git a/openmc/lattice.py b/openmc/lattice.py index 464daaf7e..b067b47a3 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -34,16 +34,16 @@ class Lattice(IDManagerMixin, ABC): Name of the lattice pitch : Iterable of float Pitch of the lattice in each direction in cm - outer : openmc.Universe + outer : openmc.UniverseBase A universe to fill all space outside the lattice - universes : Iterable of Iterable of openmc.Universe + universes : Iterable of Iterable of openmc.UniverseBase A two-or three-dimensional list/array of universes filling each element of the lattice """ next_id = 1 - used_ids = openmc.Universe.used_ids + used_ids = openmc.UniverseBase.used_ids def __init__(self, lattice_id=None, name=''): # Initialize Lattice class attributes @@ -79,7 +79,7 @@ class Lattice(IDManagerMixin, ABC): @outer.setter def outer(self, outer): - cv.check_type('outer universe', outer, openmc.Universe) + cv.check_type('outer universe', outer, openmc.UniverseBase) self._outer = outer @staticmethod @@ -92,7 +92,7 @@ class Lattice(IDManagerMixin, ABC): Group in HDF5 file universes : dict Dictionary mapping universe IDs to instances of - :class:`openmc.Universe`. + :class:`openmc.UniverseBase`. Returns ------- @@ -115,20 +115,20 @@ class Lattice(IDManagerMixin, ABC): ------- universes : collections.OrderedDict Dictionary whose keys are universe IDs and values are - :class:`openmc.Universe` instances + :class:`openmc.UniverseBase` instances """ univs = OrderedDict() for k in range(len(self._universes)): for j in range(len(self._universes[k])): - if isinstance(self._universes[k][j], openmc.Universe): + if isinstance(self._universes[k][j], openmc.UniverseBase): u = self._universes[k][j] univs[u._id] = u else: for i in range(len(self._universes[k][j])): u = self._universes[k][j][i] - assert isinstance(u, openmc.Universe) + assert isinstance(u, openmc.UniverseBase) univs[u._id] = u if self.outer is not None: @@ -246,7 +246,7 @@ class Lattice(IDManagerMixin, ABC): Returns ------- - openmc.Universe + openmc.UniverseBase Universe with given indices """ @@ -370,9 +370,9 @@ class RectLattice(Lattice): pitch : Iterable of float Pitch of the lattice in the x, y, and (if applicable) z directions in cm. - outer : openmc.Universe + outer : openmc.UniverseBase A universe to fill all space outside the lattice - universes : Iterable of Iterable of openmc.Universe + universes : Iterable of Iterable of openmc.UniverseBase A two- or three-dimensional list/array of universes filling each element of the lattice. The first dimension corresponds to the z-direction (if applicable), the second dimension corresponds to the y-direction, and @@ -633,11 +633,11 @@ class RectLattice(Lattice): cv.check_value('strategy', strategy, ('degenerate', 'lns')) cv.check_type('universes_to_ignore', universes_to_ignore, Iterable, - openmc.Universe) + openmc.UniverseBase) cv.check_type('materials_to_clone', materials_to_clone, Iterable, openmc.Material) cv.check_type('lattice_neighbors', lattice_neighbors, Iterable, - openmc.Universe) + openmc.UniverseBase) cv.check_value('number of lattice_neighbors', len(lattice_neighbors), (0, 8)) cv.check_type('key', key, types.FunctionType) @@ -973,7 +973,7 @@ class RectLattice(Lattice): Group in HDF5 file universes : dict Dictionary mapping universe IDs to instances of - :class:`openmc.Universe`. + :class:`openmc.UniverseBase`. Returns ------- @@ -999,7 +999,7 @@ class RectLattice(Lattice): lattice.outer = universes[outer] # Build array of Universe pointers for the Lattice - uarray = np.empty(universe_ids.shape, dtype=openmc.Universe) + uarray = np.empty(universe_ids.shape, dtype=openmc.UniverseBase) for z in range(universe_ids.shape[0]): for y in range(universe_ids.shape[1]): @@ -1054,9 +1054,9 @@ class HexLattice(Lattice): Pitch of the lattice in cm. The first item in the iterable specifies the pitch in the radial direction and, if the lattice is 3D, the second item in the iterable specifies the pitch in the axial direction. - outer : openmc.Universe + outer : openmc.UniverseBase A universe to fill all space outside the lattice - universes : Nested Iterable of openmc.Universe + universes : Nested Iterable of openmc.UniverseBase A two- or three-dimensional list/array of universes filling each element of the lattice. Each sub-list corresponds to one ring of universes and should be ordered from outermost ring to innermost ring. The universes @@ -1173,7 +1173,7 @@ class HexLattice(Lattice): @property def ndim(self): - return 2 if isinstance(self.universes[0][0], openmc.Universe) else 3 + return 2 if isinstance(self.universes[0][0], openmc.UniverseBase) else 3 @center.setter def center(self, center): @@ -1196,7 +1196,7 @@ class HexLattice(Lattice): @Lattice.universes.setter def universes(self, universes): - cv.check_iterable_type('lattice universes', universes, openmc.Universe, + cv.check_iterable_type('lattice universes', universes, openmc.UniverseBase, min_depth=2, max_depth=3) self._universes = universes @@ -2052,7 +2052,7 @@ class HexLattice(Lattice): Group in HDF5 file universes : dict Dictionary mapping universe IDs to instances of - :class:`openmc.Universe`. + :class:`openmc.UniverseBase`. Returns ------- diff --git a/openmc/material.py b/openmc/material.py index a2c6138d2..c84ff8c27 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -1,3 +1,4 @@ +from __future__ import annotations from collections import OrderedDict, defaultdict, namedtuple, Counter from collections.abc import Iterable from copy import deepcopy @@ -6,11 +7,11 @@ from pathlib import Path import re import typing # imported separately as py3.8 requires typing.Iterable import warnings -from typing import Optional +from typing import Optional, List, Union, Dict from xml.etree import ElementTree as ET -import h5py import numpy as np +import h5py import openmc import openmc.data @@ -134,7 +135,7 @@ class Material(IDManagerMixin): # If specified, a list of table names self._sab = [] - def __repr__(self): + def __repr__(self) -> str: string = 'Material\n' string += '{: <16}=\t{}\n'.format('\tID', self._id) string += '{: <16}=\t{}\n'.format('\tName', self._name) @@ -166,34 +167,34 @@ class Material(IDManagerMixin): return string @property - def name(self): + def name(self) -> Optional[str]: return self._name @property - def temperature(self): + def temperature(self) -> Optional[float]: return self._temperature @property - def density(self): + def density(self) -> Optional[float]: return self._density @property - def density_units(self): + def density_units(self) -> str: return self._density_units @property - def depletable(self): + def depletable(self) -> bool: return self._depletable @property - def paths(self): + def paths(self) -> List[str]: if self._paths is None: raise ValueError('Material instance paths have not been determined. ' 'Call the Geometry.determine_paths() method.') return self._paths @property - def num_instances(self): + def num_instances(self) -> int: if self._num_instances is None: raise ValueError( 'Number of material instances have not been determined. Call ' @@ -201,15 +202,15 @@ class Material(IDManagerMixin): return self._num_instances @property - def nuclides(self): + def nuclides(self) -> List[namedtuple]: return self._nuclides @property - def isotropic(self): + def isotropic(self) -> List[str]: return self._isotropic @property - def average_molar_mass(self): + def average_molar_mass(self) -> float: # Using the sum of specified atomic or weight amounts as a basis, sum # the mass and moles of the material mass = 0. @@ -226,11 +227,11 @@ class Material(IDManagerMixin): return mass / moles @property - def volume(self): + def volume(self) -> Optional[float]: return self._volume @property - def ncrystal_cfg(self): + def ncrystal_cfg(self) -> Optional[str]: return self._ncrystal_cfg @name.setter @@ -267,7 +268,7 @@ class Material(IDManagerMixin): self._isotropic = list(isotropic) @property - def fissionable_mass(self): + def fissionable_mass(self) -> float: if self.volume is None: raise ValueError("Volume must be set in order to determine mass.") density = 0.0 @@ -291,7 +292,7 @@ class Material(IDManagerMixin): return openmc.data.combine_distributions(dists, probs) if dists else None @classmethod - def from_hdf5(cls, group: h5py.Group): + def from_hdf5(cls, group: h5py.Group) -> Material: """Create material from HDF5 group Parameters @@ -346,7 +347,7 @@ class Material(IDManagerMixin): return material @classmethod - def from_ncrystal(cls, cfg, **kwargs): + def from_ncrystal(cls, cfg, **kwargs) -> Material: """Create material from NCrystal configuration string. Density, temperature, and material composition, and (ultimately) thermal @@ -881,7 +882,7 @@ class Material(IDManagerMixin): def make_isotropic_in_lab(self): self.isotropic = [x.name for x in self._nuclides] - def get_elements(self): + def get_elements(self) -> List[str]: """Returns all elements in the material .. versionadded:: 0.12 @@ -895,7 +896,7 @@ class Material(IDManagerMixin): return sorted({re.split(r'(\d+)', i)[0] for i in self.get_nuclides()}) - def get_nuclides(self, element: Optional[str] = None): + def get_nuclides(self, element: Optional[str] = None) -> List[str]: """Returns a list of all nuclides in the material, if the element argument is specified then just nuclides of that element are returned. @@ -925,7 +926,7 @@ class Material(IDManagerMixin): return matching_nuclides - def get_nuclide_densities(self): + def get_nuclide_densities(self) -> Dict[str, tuple]: """Returns all nuclides in the material and their densities Returns @@ -944,7 +945,7 @@ class Material(IDManagerMixin): return nuclides - def get_nuclide_atom_densities(self, nuclide: Optional[str] = None): + def get_nuclide_atom_densities(self, nuclide: Optional[str] = None) -> Dict[str, float]: """Returns one or all nuclides in the material and their atomic densities in units of atom/b-cm @@ -1030,7 +1031,7 @@ class Material(IDManagerMixin): return nuclides def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False, - volume: Optional[float] = None): + volume: Optional[float] = None) -> Union[Dict[str, float], float]: """Returns the activity of the material or for each nuclide in the material in units of [Bq], [Bq/g] or [Bq/cm3]. @@ -1077,7 +1078,7 @@ class Material(IDManagerMixin): return activity if by_nuclide else sum(activity.values()) def get_decay_heat(self, units: str = 'W', by_nuclide: bool = False, - volume: Optional[float] = None): + volume: Optional[float] = None) -> Union[Dict[str, float], float]: """Returns the decay heat of the material or for each nuclide in the material in units of [W], [W/g] or [W/cm3]. @@ -1125,7 +1126,7 @@ class Material(IDManagerMixin): return decayheat if by_nuclide else sum(decayheat.values()) - def get_nuclide_atoms(self, volume: Optional[float] = None): + def get_nuclide_atoms(self, volume: Optional[float] = None) -> Dict[str, float]: """Return number of atoms of each nuclide in the material .. versionadded:: 0.13.1 @@ -1154,7 +1155,7 @@ class Material(IDManagerMixin): atoms[nuclide] = 1.0e24 * atom_per_bcm * volume return atoms - def get_mass_density(self, nuclide: Optional[str] = None): + def get_mass_density(self, nuclide: Optional[str] = None) -> float: """Return mass density of one or all nuclides Parameters @@ -1176,7 +1177,7 @@ class Material(IDManagerMixin): mass_density += density_i return mass_density - def get_mass(self, nuclide: Optional[str] = None, volume: Optional[float] = None): + def get_mass(self, nuclide: Optional[str] = None, volume: Optional[float] = None) -> float: """Return mass of one or all nuclides. Note that this method requires that the :attr:`Material.volume` has @@ -1206,7 +1207,7 @@ class Material(IDManagerMixin): raise ValueError("Volume must be set in order to determine mass.") return volume*self.get_mass_density(nuclide) - def clone(self, memo: Optional[dict] = None): + def clone(self, memo: Optional[dict] = None) -> Material: """Create a copy of this material with a new unique ID. Parameters @@ -1245,7 +1246,7 @@ class Material(IDManagerMixin): return memo[self] - def _get_nuclide_xml(self, nuclide: str): + def _get_nuclide_xml(self, nuclide: str) -> ET.Element: xml_element = ET.Element("nuclide") xml_element.set("name", nuclide.name) @@ -1256,19 +1257,19 @@ class Material(IDManagerMixin): return xml_element - def _get_macroscopic_xml(self, macroscopic: str): + def _get_macroscopic_xml(self, macroscopic: str) -> ET.Element: xml_element = ET.Element("macroscopic") xml_element.set("name", macroscopic) return xml_element - def _get_nuclides_xml(self, nuclides: typing.Iterable[str]): + def _get_nuclides_xml(self, nuclides: typing.Iterable[str]) -> List[ET.Element]: xml_elements = [] for nuclide in nuclides: xml_elements.append(self._get_nuclide_xml(nuclide)) return xml_elements - def to_xml_element(self): + def to_xml_element(self) -> ET.Element: """Return XML representation of the material Returns @@ -1338,7 +1339,7 @@ class Material(IDManagerMixin): @classmethod def mix_materials(cls, materials, fracs: typing.Iterable[float], - percent_type: str = 'ao', name: Optional[str] = None): + percent_type: str = 'ao', name: Optional[str] = None) -> Material: """Mix materials together based on atom, weight, or volume fractions .. versionadded:: 0.12 @@ -1436,7 +1437,7 @@ class Material(IDManagerMixin): return new_mat @classmethod - def from_xml_element(cls, elem: ET.Element): + def from_xml_element(cls, elem: ET.Element) -> Material: """Generate material from an XML element Parameters @@ -1530,7 +1531,7 @@ class Materials(cv.CheckedList): self += materials @property - def cross_sections(self): + def cross_sections(self) -> Optional[Path]: return self._cross_sections @cross_sections.setter @@ -1638,7 +1639,7 @@ class Materials(cv.CheckedList): self._write_xml(fh) @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem) -> Material: """Generate materials collection from XML file Parameters @@ -1665,7 +1666,7 @@ class Materials(cv.CheckedList): return materials @classmethod - def from_xml(cls, path: PathLike = 'materials.xml'): + def from_xml(cls, path: PathLike = 'materials.xml') -> Material: """Generate materials collection from XML file Parameters diff --git a/openmc/model/model.py b/openmc/model/model.py index 4202cf60b..933e1172e 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1,3 +1,4 @@ +from __future__ import annotations from collections.abc import Iterable from contextlib import contextmanager from functools import lru_cache @@ -7,6 +8,7 @@ from numbers import Integral from tempfile import NamedTemporaryFile import warnings from xml.etree import ElementTree as ET +from typing import Optional, Dict import h5py @@ -93,27 +95,27 @@ class Model: self.plots = plots @property - def geometry(self): + def geometry(self) -> Optional[openmc.Geometry]: return self._geometry @property - def materials(self): + def materials(self) -> Optional[openmc.Materials]: return self._materials @property - def settings(self): + def settings(self) -> Optional[openmc.Settings]: return self._settings @property - def tallies(self): + def tallies(self) -> Optional[openmc.Tallies]: return self._tallies @property - def plots(self): + def plots(self) -> Optional[openmc.Plots]: return self._plots @property - def is_initialized(self): + def is_initialized(self) -> bool: try: import openmc.lib return openmc.lib.is_initialized @@ -122,7 +124,7 @@ class Model: @property @lru_cache(maxsize=None) - def _materials_by_id(self): + def _materials_by_id(self) -> dict: """Dictionary mapping material ID --> material""" if self.materials: mats = self.materials @@ -132,14 +134,14 @@ class Model: @property @lru_cache(maxsize=None) - def _cells_by_id(self): + def _cells_by_id(self) -> dict: """Dictionary mapping cell ID --> cell""" cells = self.geometry.get_all_cells() return {cell.id: cell for cell in cells.values()} @property @lru_cache(maxsize=None) - def _cells_by_name(self): + def _cells_by_name(self) -> Dict[int, openmc.Cell]: # Get the names maps, but since names are not unique, store a set for # each name key. In this way when the user requests a change by a name, # the change will be applied to all of the same name. @@ -152,7 +154,7 @@ class Model: @property @lru_cache(maxsize=None) - def _materials_by_name(self): + def _materials_by_name(self) -> Dict[int, openmc.Material]: if self.materials is None: mats = self.geometry.get_all_materials().values() else: @@ -207,7 +209,7 @@ class Model: @classmethod def from_xml(cls, geometry='geometry.xml', materials='materials.xml', settings='settings.xml', tallies='tallies.xml', - plots='plots.xml'): + plots='plots.xml') -> Model: """Create model from existing XML files Parameters @@ -597,12 +599,13 @@ class Model: def run(self, particles=None, threads=None, geometry_debug=False, restart_file=None, tracks=False, output=True, cwd='.', openmc_exec='openmc', mpi_args=None, event_based=None, - export_model_xml=True): - """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 call. In both cases this method returns the path to the last - statepoint file generated. + export_model_xml=True, **export_kwargs): + """Run 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 call. 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. @@ -630,27 +633,30 @@ class Model: output : bool, optional Capture OpenMC output from standard out cwd : str, optional - Path to working directory to run in. Defaults to the current - working directory. + Path to working directory to run in. Defaults to the current working + directory. openmc_exec : str, optional Path to OpenMC executable. Defaults to 'openmc'. mpi_args : list of str, optional - MPI execute command and any additional MPI arguments to pass, - e.g. ['mpiexec', '-n', '8']. + MPI execute command and any additional MPI arguments to pass, e.g. + ['mpiexec', '-n', '8']. event_based : None or bool, optional - Turns on event-based parallelism if True. If None, the value in - the Settings will be used. + Turns on event-based parallelism if True. If None, the value in the + Settings will be used. export_model_xml : bool, optional - Exports a single model.xml file rather than separate files. - Defaults to True. + Exports a single model.xml file rather than separate files. Defaults + to True. .. versionadded:: 0.13.3 + **export_kwargs + Keyword arguments passed to either :meth:`Model.export_to_model_xml` + or :meth:`Model.export_to_xml`. Returns ------- Path - Path to the last statepoint written by this run - (None if no statepoint was written) + Path to the last statepoint written by this run (None if no + statepoint was written) """ @@ -695,9 +701,9 @@ class Model: else: # Then run via the command line if export_model_xml: - self.export_to_model_xml() + self.export_to_model_xml(**export_kwargs) else: - self.export_to_xml() + self.export_to_xml(**export_kwargs) openmc.run(particles, threads, geometry_debug, restart_file, tracks, output, Path('.'), openmc_exec, mpi_args, event_based) diff --git a/openmc/plotter.py b/openmc/plotter.py index a37b6f042..0362a8d17 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -215,8 +215,9 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, return fig + def calculate_cexs(this, types, temperature=294., sab_name=None, - cross_sections=None, enrichment=None): + cross_sections=None, enrichment=None, ncrystal_cfg=None): """Calculates continuous-energy cross sections of a requested type. Parameters @@ -239,6 +240,8 @@ def calculate_cexs(this, types, temperature=294., sab_name=None, Enrichment for U235 in weight percent. For example, input 4.95 for 4.95 weight percent enriched U. Default is None (natural composition). + ncrystal_cfg : str, optional + Configuration string for NCrystal material. Returns ------- @@ -262,10 +265,10 @@ def calculate_cexs(this, types, temperature=294., sab_name=None, energy_grid, data = _calculate_cexs_elem_mat( this, types, temperature, cross_sections, sab_name, enrichment ) - else: energy_grid, xs = _calculate_cexs_nuclide( - this, types, temperature, sab_name, cross_sections + this, types, temperature, sab_name, cross_sections, + ncrystal_cfg ) # Convert xs (Iterable of Callable) to a grid of cross section values @@ -282,7 +285,7 @@ def calculate_cexs(this, types, temperature=294., sab_name=None, def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None, - cross_sections=None): + cross_sections=None, ncrystal_cfg=None): """Calculates continuous-energy cross sections of a requested type. Parameters @@ -303,6 +306,8 @@ def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None, Name of S(a,b) library to apply to MT=2 data when applicable. cross_sections : str, optional Location of cross_sections.xml file. Default is None. + ncrystal_cfg : str, optional + Configuration string for NCrystal material. Returns ------- @@ -420,6 +425,19 @@ def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None, [sab_sum, nuc[mt].xs[nucT]], [sab_Emax]) funcs.append(pw_funcs) + elif ncrystal_cfg: + import NCrystal + nc_scatter = NCrystal.createScatter(ncrystal_cfg) + nc_func = nc_scatter.crossSectionNonOriented + nc_emax = 5 # eV # this should be obtained from NCRYSTAL_MAX_ENERGY + energy_grid = np.union1d(np.geomspace(min(energy_grid), + 1.1*nc_emax, + 1000),energy_grid) # NCrystal does not have + # an intrinsic energy grid + pw_funcs = openmc.data.Regions1D( + [nc_func, nuc[mt].xs[nucT]], + [nc_emax]) + funcs.append(pw_funcs) else: funcs.append(nuc[mt].xs[nucT]) elif mt in nuc: @@ -522,12 +540,15 @@ def _calculate_cexs_elem_mat(this, types, temperature=294., # Load the library library = openmc.data.DataLibrary.from_xml(cross_sections) + ncrystal_cfg = None if isinstance(this, openmc.Material): # Expand elements in to nuclides with atomic densities nuc_fractions = this.get_nuclide_atom_densities() # Create a dict of [nuclide name] = nuclide object to carry forward # with a common nuclides format between openmc.Material and Elements nuclides = {nuclide: nuclide for nuclide in nuc_fractions} + # Add NCrystal cfg string if it exists + ncrystal_cfg = this.ncrystal_cfg else: # Expand elements in to nuclides with atomic densities nuclides = openmc.Element(this).expand(1., 'ao', enrichment=enrichment, @@ -563,7 +584,9 @@ def _calculate_cexs_elem_mat(this, types, temperature=294., name = nuclide[0] nuc = nuclide[1] sab_tab = sabs[name] - temp_E, temp_xs = calculate_cexs(nuc, types, T, sab_tab, cross_sections) + temp_E, temp_xs = calculate_cexs(nuc, types, T, sab_tab, cross_sections, + ncrystal_cfg=ncrystal_cfg + ) E.append(temp_E) # Since the energy grids are different, store the cross sections as # a tabulated function so they can be calculated on any grid needed. diff --git a/openmc/universe.py b/openmc/universe.py index 2ed96f4bc..67296dda8 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -1,3 +1,4 @@ +import typing from abc import ABC, abstractmethod from collections import OrderedDict from collections.abc import Iterable @@ -6,6 +7,7 @@ from numbers import Integral, Real from pathlib import Path from tempfile import TemporaryDirectory from xml.etree import ElementTree as ET +from warnings import warn import h5py import numpy as np @@ -87,7 +89,8 @@ class UniverseBase(ABC, IDManagerMixin): self._volume = volume_calc.volumes[self.id].n self._atoms = volume_calc.atoms[self.id] else: - raise ValueError('No volume information found for this universe.') + raise ValueError( + 'No volume information found for this universe.') else: raise ValueError('No volume information found for this universe.') @@ -168,7 +171,7 @@ class UniverseBase(ABC, IDManagerMixin): clone._cells = OrderedDict() for cell in self._cells.values(): clone.add_cell(cell.clone(clone_materials, clone_regions, - memo)) + memo)) # Memoize the clone memo[self] = clone @@ -293,9 +296,15 @@ class Universe(UniverseBase): return [self, cell] + cell.fill.find(p) return [] + # default kwargs that are passed to plt.legend in the plot method below. + _default_legend_kwargs = {'bbox_to_anchor': ( + 1.05, 1), 'loc': 2, 'borderaxespad': 0.0} + def plot(self, origin=(0., 0., 0.), width=(1., 1.), pixels=(200, 200), basis='xy', color_by='cell', colors=None, seed=None, - openmc_exec='openmc', axes=None, **kwargs): + openmc_exec='openmc', axes=None, legend=False, + legend_kwargs=_default_legend_kwargs, outline=False, + **kwargs): """Display a slice plot of the universe. Parameters @@ -329,6 +338,18 @@ class Universe(UniverseBase): Axes to draw to .. versionadded:: 0.13.1 + legend : bool + Whether a legend showing material or cell names should be drawn + + .. versionadded:: 0.13.4 + legend_kwargs : dict + Keyword arguments passed to :func:`matplotlib.pyplot.legend`. + + .. versionadded:: 0.13.4 + outline : bool + Whether outlines between color boundaries should be drawn + + .. versionadded:: 0.13.4 **kwargs Keyword arguments passed to :func:`matplotlib.pyplot.imshow` @@ -339,15 +360,19 @@ class Universe(UniverseBase): """ import matplotlib.image as mpimg + import matplotlib.patches as mpatches import matplotlib.pyplot as plt # Determine extents of plot if basis == 'xy': x, y = 0, 1 + xlabel, ylabel = 'x [cm]', 'y [cm]' elif basis == 'yz': x, y = 1, 2 + xlabel, ylabel = 'y [cm]', 'z [cm]' elif basis == 'xz': x, y = 0, 2 + xlabel, ylabel = 'x [cm]', 'z [cm]' x_min = origin[x] - 0.5*width[0] x_max = origin[x] + 0.5*width[0] y_min = origin[y] - 0.5*width[1] @@ -391,11 +416,57 @@ class Universe(UniverseBase): if axes is None: px = 1/plt.rcParams['figure.dpi'] fig, axes = plt.subplots() + axes.set_xlabel(xlabel) + axes.set_ylabel(ylabel) params = fig.subplotpars width = pixels[0]*px/(params.right - params.left) height = pixels[0]*px/(params.top - params.bottom) fig.set_size_inches(width, height) + if outline: + # Combine R, G, B values into a single int + rgb = (img * 256).astype(int) + image_value = (rgb[..., 0] << 16) + (rgb[..., 1] << 8) + (rgb[..., 2]) + + axes.contour( + image_value, + origin="upper", + colors="k", + linestyles="solid", + linewidths=1, + levels=np.unique(image_value), + extent=(x_min, x_max, y_min, y_max), + ) + + # add legend showing which colors represent which material + # or cell if that was requested + if legend: + if plot.colors is None: + raise ValueError("Must pass 'colors' dictionary if you " + "are adding a legend via legend=True.") + + if color_by == "cell": + expected_key_type = openmc.Cell + else: + expected_key_type = openmc.Material + + patches = [] + for key, color in plot.colors.items(): + + if isinstance(key, int): + raise TypeError( + "Cannot use IDs in colors dict for auto legend.") + elif not isinstance(key, expected_key_type): + raise TypeError( + "Color dict key type does not match color_by") + + # this works whether we're doing cells or materials + label = key.name if key.name != '' else key.id + key_patch = mpatches.Patch(color=color, label=label) + patches.append(key_patch) + + axes.legend(handles=patches, **legend_kwargs) + # Plot image and return the axes return axes.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs) @@ -738,8 +809,9 @@ class DAGMCUniverse(UniverseBase): @property def material_names(self): dagmc_file_contents = h5py.File(self.filename) - material_tags_hex=dagmc_file_contents['/tstt/tags/NAME'].get('values') - material_tags_ascii=[] + material_tags_hex = dagmc_file_contents['/tstt/tags/NAME'].get( + 'values') + material_tags_ascii = [] for tag in material_tags_hex: candidate_tag = tag.tobytes().decode().replace('\x00', '') # tags might be for temperature or reflective surfaces @@ -905,7 +977,8 @@ class DAGMCUniverse(UniverseBase): openmc.Universe Universe instance """ - bounding_cell = openmc.Cell(fill=self, cell_id=bounding_cell_id, region=self.bounding_region(**kwargs)) + bounding_cell = openmc.Cell( + fill=self, cell_id=bounding_cell_id, region=self.bounding_region(**kwargs)) return openmc.Universe(cells=[bounding_cell]) @classmethod @@ -972,4 +1045,4 @@ class DAGMCUniverse(UniverseBase): clone.volume = self.volume clone.auto_geom_ids = self.auto_geom_ids clone.auto_mat_ids = self.auto_mat_ids - return clone \ No newline at end of file + return clone diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 65b9b7bb0..2993db089 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -1,14 +1,18 @@ from __future__ import annotations from collections.abc import Iterable from numbers import Real, Integral -from typing import Iterable, List +import pathlib +import typing +from typing import Iterable, List, Optional, Union, Dict from xml.etree import ElementTree as ET import numpy as np +import h5py from openmc.filter import _PARTICLES from openmc.mesh import MeshBase, RectilinearMesh, UnstructuredMesh import openmc.checkvalue as cv +from openmc.checkvalue import PathLike from ._xml import get_text from .mixin import IDManagerMixin @@ -103,16 +107,20 @@ class WeightWindows(IDManagerMixin): next_id = 1 used_ids = set() - def __init__(self, mesh, lower_ww_bounds, - upper_ww_bounds=None, - upper_bound_ratio=None, - energy_bounds=None, - particle_type='neutron', - survival_ratio=3, - max_lower_bound_ratio=None, - max_split=10, - weight_cutoff=1.e-38, - id=None): + def __init__( + self, + mesh: MeshBase, + lower_ww_bounds: Iterable[float], + upper_ww_bounds: Optional[Iterable[float]] = None, + upper_bound_ratio: Optional[float] = None, + energy_bounds: Optional[Iterable[Real]] = None, + particle_type: str = 'neutron', + survival_ratio: float = 3, + max_lower_bound_ratio: Optional[float] = None, + max_split: int = 10, + weight_cutoff: float = 1.e-38, + id: Optional[int] = None + ): self.mesh = mesh self.id = id self.particle_type = particle_type @@ -161,7 +169,7 @@ class WeightWindows(IDManagerMixin): string += '{: <16}=\t{}\n'.format('\tWeight Cutoff', self._weight_cutoff) return string - def __eq__(self, other) -> bool: + def __eq__(self, other: WeightWindows) -> bool: # ensure that `other` is a WeightWindows object if not isinstance(other, WeightWindows): return False @@ -195,7 +203,7 @@ class WeightWindows(IDManagerMixin): return self._mesh @mesh.setter - def mesh(self, mesh): + def mesh(self, mesh: MeshBase): cv.check_type('Weight window mesh', mesh, MeshBase) self._mesh = mesh @@ -204,7 +212,7 @@ class WeightWindows(IDManagerMixin): return self._particle_type @particle_type.setter - def particle_type(self, pt): + def particle_type(self, pt: str): cv.check_value('Particle type', pt, _PARTICLES) self._particle_type = pt @@ -213,7 +221,7 @@ class WeightWindows(IDManagerMixin): return self._energy_bounds @energy_bounds.setter - def energy_bounds(self, bounds): + def energy_bounds(self, bounds: Iterable[float]): cv.check_type('Energy bounds', bounds, Iterable, Real) self._energy_bounds = np.asarray(bounds) @@ -228,7 +236,7 @@ class WeightWindows(IDManagerMixin): return self._lower_ww_bounds @lower_ww_bounds.setter - def lower_ww_bounds(self, bounds): + def lower_ww_bounds(self, bounds: Iterable[float]): cv.check_iterable_type('Lower WW bounds', bounds, Real, @@ -247,7 +255,7 @@ class WeightWindows(IDManagerMixin): return self._upper_ww_bounds @upper_ww_bounds.setter - def upper_ww_bounds(self, bounds): + def upper_ww_bounds(self, bounds: Iterable[float]): cv.check_iterable_type('Upper WW bounds', bounds, Real, @@ -266,7 +274,7 @@ class WeightWindows(IDManagerMixin): return self._survival_ratio @survival_ratio.setter - def survival_ratio(self, val): + def survival_ratio(self, val: float): cv.check_type('Survival ratio', val, Real) cv.check_greater_than('Survival ratio', val, 1.0, True) self._survival_ratio = val @@ -276,7 +284,7 @@ class WeightWindows(IDManagerMixin): return self._max_lower_bound_ratio @max_lower_bound_ratio.setter - def max_lower_bound_ratio(self, val): + def max_lower_bound_ratio(self, val: float): cv.check_type('Maximum lower bound ratio', val, Real) cv.check_greater_than('Maximum lower bound ratio', val, 1.0) self._max_lower_bound_ratio = val @@ -286,7 +294,7 @@ class WeightWindows(IDManagerMixin): return self._max_split @max_split.setter - def max_split(self, val): + def max_split(self, val: int): cv.check_type('Max split', val, Integral) self._max_split = val @@ -295,7 +303,7 @@ class WeightWindows(IDManagerMixin): return self._weight_cutoff @weight_cutoff.setter - def weight_cutoff(self, cutoff): + def weight_cutoff(self, cutoff: float): cv.check_type('Weight cutoff', cutoff, Real) cv.check_greater_than('Weight cutoff', cutoff, 0.0, True) self._weight_cutoff = cutoff @@ -343,7 +351,7 @@ class WeightWindows(IDManagerMixin): return element @classmethod - def from_xml_element(cls, elem, root) -> WeightWindows: + def from_xml_element(cls, elem: ET.Element, root: ET.Element) -> WeightWindows: """Generate weight window settings from an XML element Parameters @@ -398,7 +406,7 @@ class WeightWindows(IDManagerMixin): ) @classmethod - def from_hdf5(cls, group, meshes) -> WeightWindows: + def from_hdf5(cls, group: h5py.Group, meshes: Dict[int, MeshBase]) -> WeightWindows: """Create weight windows from HDF5 group Parameters @@ -443,7 +451,7 @@ class WeightWindows(IDManagerMixin): ) -def wwinp_to_wws(path) -> List[WeightWindows]: +def wwinp_to_wws(path: PathLike) -> List[WeightWindows]: """Create WeightWindows instances from a wwinp file .. versionadded:: 0.13.1 diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 04c8424aa..4d63f0892 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -587,32 +587,42 @@ std::pair DAGCell::distance( if (!dag_univ) fatal_error("DAGMC call made for particle in a non-DAGMC universe"); - moab::ErrorCode rval; + // initialize to lost particle conditions + int surf_idx = -1; + double dist = INFINITY; + moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_); moab::EntityHandle hit_surf; - double dist; + + // create the ray double pnt[3] = {r.x, r.y, r.z}; double dir[3] = {u.x, u.y, u.z}; - rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history()); - MB_CHK_ERR_CONT(rval); - int surf_idx; + MB_CHK_ERR_CONT( + dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history())); if (hit_surf != 0) { surf_idx = dag_univ->surf_idx_offset_ + dagmc_ptr_->index_by_handle(hit_surf); - } else { - // indicate that particle is lost - surf_idx = -1; - dist = INFINITY; - if (settings::run_mode == RunMode::PLOTTING) return {dist, surf_idx}; - if (!dagmc_ptr_->is_implicit_complement(vol) || - model::universe_map[dag_univ->id_] == model::root_universe) { - std::string material_id = p->material() == MATERIAL_VOID - ? "-1 (VOID)" - : std::to_string(model::materials[p->material()]->id()); - p->mark_as_lost( - fmt::format("No intersection found with DAGMC cell {}, material {}", - id_, material_id)); - } + } else if (!dagmc_ptr_->is_implicit_complement(vol) || + is_root_universe(dag_univ->id_)) { + // surface boundary conditions are ignored for projection plotting, meaning + // that the particle may move through the graveyard (bounding) volume and + // into the implicit complement on the other side where no intersection will + // be found. Treating this as a lost particle is problematic when plotting. + // Instead, the infinite distance and invalid surface index are returned. + if (settings::run_mode == RunMode::PLOTTING) return {INFTY, -1}; + + // the particle should be marked as lost immediately if an intersection + // isn't found in a volume that is not the implicit complement. In the case + // that the DAGMC model is the root universe of the geometry, even a missing + // intersection in the implicit complement should trigger this condition. + std::string material_id = + p->material() == MATERIAL_VOID + ? "-1 (VOID)" + : std::to_string(model::materials[p->material()]->id()); + auto lost_particle_msg = fmt::format( + "No intersection found with DAGMC cell {}, filled with material {}", id_, + material_id); + p->mark_as_lost(lost_particle_msg); } return {dist, surf_idx}; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index eabdc9f70..3ab5fcc2d 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -610,6 +610,11 @@ int maximum_levels(int32_t univ) return levels_below; } +bool is_root_universe(int32_t univ_id) +{ + return model::universe_map[univ_id] == model::root_universe; +} + //============================================================================== void free_memory_geometry() diff --git a/src/initialize.cpp b/src/initialize.cpp index 6383af64e..9c0a05d37 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -358,7 +358,7 @@ bool read_model_xml() for (const auto& input : other_inputs) { if (file_exists(settings::path_input + input)) { warning((fmt::format("Other XML file input(s) are present. These files " - "will be ignored in favor of the {} file.", + "may be ignored in favor of the {} file.", model_filename))); break; } @@ -392,8 +392,16 @@ bool read_model_xml() // Initialize distribcell_filters prepare_distribcell(); - if (check_for_node(root, "plots")) + if (check_for_node(root, "plots")) { read_plots_xml(root.child("plots")); + } else { + // When no element is present in the model.xml file, check for a + // regular plots.xml file + std::string filename = settings::path_input + "plots.xml"; + if (file_exists(filename)) { + read_plots_xml(); + } + } return true; } diff --git a/src/output.cpp b/src/output.cpp index 666ae15b3..2b0a4024c 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -268,7 +268,7 @@ void print_usage() { if (mpi::master) { fmt::print( - "Usage: openmc [options] [directory]\n\n" + "Usage: openmc [options] [path]\n\n" "Options:\n" " -c, --volume Run in stochastic volume calculation mode\n" " -g, --geometry-debug Run with geometry debugging on\n" diff --git a/src/settings.cpp b/src/settings.cpp index 2f8fbd1d0..6fa01708e 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -223,16 +223,11 @@ void read_settings_xml() std::string filename = settings::path_input + "settings.xml"; if (!file_exists(filename)) { if (run_mode != RunMode::PLOTTING) { - fatal_error( - fmt::format("Settings XML file '{}' does not exist! In order " - "to run OpenMC, you first need a set of input files; at a " - "minimum, this " - "includes settings.xml, geometry.xml, and materials.xml " - "or a single XML file containing all of these files. " - "Please consult " - "the user's guide at https://docs.openmc.org for further " - "information.", - filename)); + fatal_error("Could not find any XML input files! In order to run OpenMC, " + "you first need a set of input files; at a minimum, this " + "includes settings.xml, geometry.xml, and materials.xml or a " + "single model XML file. Please consult the user's guide at " + "https://docs.openmc.org for further information."); } else { // The settings.xml file is optional if we just want to make a plot. return; diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp index 79af65879..87f298baa 100644 --- a/src/tallies/trigger.cpp +++ b/src/tallies/trigger.cpp @@ -39,7 +39,8 @@ std::pair get_tally_uncertainty( auto mean = sum / n; // if the result has no contributions, return an invalid pair - if (mean == 0) return {-1 , -1}; + if (mean == 0) + return {-1, -1}; double std_dev = std::sqrt((sum_sq / n - mean * mean) / (n - 1)); double rel_err = (mean != 0.) ? std_dev / std::abs(mean) : 0.; @@ -193,8 +194,10 @@ void check_triggers() keff_ratio); } else { if (tally_ratio == INFINITY) { - msg = fmt::format("Triggers unsatisfied, no result tallied for score {} in tally {}", reaction_name(score), tally_id); - } else{ + msg = fmt::format( + "Triggers unsatisfied, no result tallied for score {} in tally {}", + reaction_name(score), tally_id); + } else { msg = fmt::format( "Triggers unsatisfied, max unc./thresh. is {} for {} in tally {}", tally_ratio, reaction_name(score), tally_id); @@ -211,13 +214,20 @@ void check_triggers() auto n_pred_batches = static_cast(n_active * max_ratio * max_ratio) + settings::n_inactive + 1; - std::string msg = - fmt::format("The estimated number of batches is {}", n_pred_batches); - if (n_pred_batches > settings::n_max_batches) { - msg.append(" --- greater than max batches"); - warning(msg); - } else { + if (max_ratio == INFINITY) { + std::string msg = + fmt::format("One or more tallies with triggers have no scores. Unable " + "to estimate the number of remaining batches."); write_message(msg, 7); + } else { + std::string msg = + fmt::format("The estimated number of batches is {}", n_pred_batches); + if (n_pred_batches > settings::n_max_batches) { + msg.append(" --- greater than max batches"); + warning(msg); + } else { + write_message(msg, 7); + } } } } diff --git a/tests/unit_tests/dagmc/broken_model.h5m b/tests/unit_tests/dagmc/broken_model.h5m new file mode 100644 index 000000000..4bc6c7d42 Binary files /dev/null and b/tests/unit_tests/dagmc/broken_model.h5m differ diff --git a/tests/unit_tests/dagmc/test_lost_particles.py b/tests/unit_tests/dagmc/test_lost_particles.py new file mode 100644 index 000000000..502bd795e --- /dev/null +++ b/tests/unit_tests/dagmc/test_lost_particles.py @@ -0,0 +1,81 @@ +import numpy as np +from pathlib import Path + +import openmc +import openmc.lib + +import pytest + +pytestmark = pytest.mark.skipif( + not openmc.lib._dagmc_enabled(), + reason="DAGMC CAD geometry is not enabled.") + + +@pytest.fixture +def broken_dagmc_model(request): + openmc.reset_auto_ids() + model = openmc.Model() + + ### MATERIALS ### + fuel = openmc.Material(name='no-void fuel') + fuel.set_density('g/cc', 10.29769) + fuel.add_nuclide('U233', 1.0) + + cladding = openmc.Material(name='clad') + cladding.set_density('g/cc', 6.55) + cladding.add_nuclide('Zr90', 1.0) + + h1 = openmc.Material(name='water') + h1.set_density('g/cc', 0.75) + h1.add_nuclide('H1', 1.0) + + model.materials = openmc.Materials([fuel, cladding, h1]) + + ### GEOMETRY ### + # create the DAGMC universe using a model that has many triangles + # removed + dagmc_file = Path(request.fspath).parent / "broken_model.h5m" + pincell_univ = openmc.DAGMCUniverse(filename=dagmc_file, auto_geom_ids=True) + + # create a 2 x 2 lattice using the DAGMC pincell + pitch = np.asarray((24.0, 24.0)) + lattice = openmc.RectLattice() + lattice.pitch = pitch + lattice.universes = [[pincell_univ] * 2] * 2 + lattice.lower_left = -pitch + + # clip the DAGMC geometry at +/- 10 cm w/ CSG planes + rpp = openmc.model.RectangularParallelepiped( + -pitch[0], pitch[0], -pitch[1], pitch[1], -10.0, 10.0, boundary_type='reflective') + bounding_cell = openmc.Cell(fill=lattice, region=-rpp) + + model.geometry = openmc.Geometry(root=[bounding_cell]) + + # settings + model.settings.particles = 100 + model.settings.batches = 10 + model.settings.inactive = 2 + model.settings.output = {'summary': False} + + model.export_to_xml() + + return model + + +def test_lost_particles(run_in_tmpdir, broken_dagmc_model): + broken_dagmc_model.export_to_xml() + # ensure that particles will be lost when cell intersections can't be found + # due to the removed triangles in this model + with pytest.raises(RuntimeError, match='Maximum number of lost particles has been reached.'): + openmc.run() + + # run this again, but with the dagmc universe as the root unvierse + for univ in broken_dagmc_model.geometry.get_all_universes().values(): + if isinstance(univ, openmc.DAGMCUniverse): + broken_dagmc_model.geometry.root_unvierse = univ + break + + broken_dagmc_model.export_to_xml() + with pytest.raises(RuntimeError, match='Maximum number of lost particles has been reached.'): + openmc.run() + diff --git a/tests/unit_tests/test_plots.py b/tests/unit_tests/test_plots.py index 15d158dd6..71ccd4a23 100644 --- a/tests/unit_tests/test_plots.py +++ b/tests/unit_tests/test_plots.py @@ -187,9 +187,9 @@ def test_plots(run_in_tmpdir): plots.export_to_xml() # from_xml - openmc.Plots.from_xml() - assert len(plots) - assert plots[0].origin == p1.origin - assert plots[0].colors == p1.colors - assert plots[0].mask_components == p1.mask_components - assert plots[1].origin == p2.origin + new_plots = openmc.Plots.from_xml() + assert len(new_plots) + assert new_plots[0].origin == p1.origin + assert new_plots[0].colors == p1.colors + assert new_plots[0].mask_components == p1.mask_components + assert new_plots[1].origin == p2.origin diff --git a/tests/unit_tests/test_universe.py b/tests/unit_tests/test_universe.py index a70bd4667..1bfbbf7ff 100644 --- a/tests/unit_tests/test_universe.py +++ b/tests/unit_tests/test_universe.py @@ -59,6 +59,7 @@ def test_plot(run_in_tmpdir, sphere_model): pixels=(10, 10), color_by='material', colors=colors, + outline=True, ) diff --git a/tools/ci/gha-install-vectfit.sh b/tools/ci/gha-install-vectfit.sh index 5884c6fac..8444c3036 100755 --- a/tools/ci/gha-install-vectfit.sh +++ b/tools/ci/gha-install-vectfit.sh @@ -39,7 +39,10 @@ cd $HOME git clone -b $XTENSOR_BLAS_BRANCH $XTENSOR_BLAS_REPO cd xtensor-blas && mkdir build && cd build && cmake .. && sudo make install +# Install wheel (remove when vectfit supports installation with build isolation) +pip install wheel + # Install vectfit cd $HOME git clone https://github.com/liangjg/vectfit.git -pip install ./vectfit +pip install --no-build-isolation ./vectfit