Merge branch 'develop' into adding_outlines_to_plot

This commit is contained in:
Jonathan Shimwell 2023-04-20 08:07:43 +01:00 committed by GitHub
commit cb7a443776
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
6 changed files with 204 additions and 117 deletions

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

@ -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
@ -604,7 +606,6 @@ class Model:
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.

View file

@ -7,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
@ -88,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.')
@ -169,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
@ -294,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, outline=False, **kwargs):
openmc_exec='openmc', axes=None, legend=False,
legend_kwargs=_default_legend_kwargs, outline=False,
**kwargs):
"""Display a slice plot of the universe.
Parameters
@ -330,13 +338,20 @@ class Universe(UniverseBase):
Axes to draw to
.. versionadded:: 0.13.1
**kwargs
Keyword arguments passed to :func:`matplotlib.pyplot.imshow`
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`
Returns
-------
@ -345,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]
@ -397,6 +416,8 @@ 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)
@ -417,6 +438,35 @@ class Universe(UniverseBase):
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)
@ -759,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
@ -926,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

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

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