Merge branch 'develop' into moving_voxel_plot_code

This commit is contained in:
Jonathan Shimwell 2023-04-20 18:19:28 +01:00 committed by GitHub
commit 63603dbabd
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
21 changed files with 452 additions and 199 deletions

View file

@ -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.
//==============================================================================

View file

@ -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

View file

@ -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
<https://www-nds.iaea.org/amdc/ame2016/AME2016-a.pdf>`_.
Atomic mass data comes from the `Atomic Mass Evaluation 2020
<https://doi.org/10.1088/1674-1137/abddaf>`_.
Parameters
----------

View file

@ -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."""

View file

@ -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
-------

View file

@ -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

View file

@ -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)

View file

@ -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.

View file

@ -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
return clone

View file

@ -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

View file

@ -587,32 +587,42 @@ std::pair<double, int32_t> 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};

View file

@ -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()

View file

@ -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 <plots> 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;
}

View file

@ -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"

View file

@ -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;

View file

@ -39,7 +39,8 @@ std::pair<double, double> 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<int>(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);
}
}
}
}

Binary file not shown.

View file

@ -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()

View file

@ -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

View file

@ -59,6 +59,7 @@ def test_plot(run_in_tmpdir, sphere_model):
pixels=(10, 10),
color_by='material',
colors=colors,
outline=True,
)

View file

@ -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