diff --git a/docs/source/conf.py b/docs/source/conf.py index 8dbbfc6d444..dc5acd5f62d 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -120,10 +120,9 @@ # -- Options for HTML output --------------------------------------------------- # The theme to use for HTML and HTML Help pages -if not on_rtd: - import sphinx_rtd_theme - html_theme = 'sphinx_rtd_theme' - html_theme_path = [sphinx_rtd_theme.get_html_theme_path()] +import sphinx_rtd_theme +html_theme = 'sphinx_rtd_theme' +html_theme_path = [sphinx_rtd_theme.get_html_theme_path()] html_logo = '_images/openmc_logo.png' diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 1016fbaa3d7..a52de1b9048 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -60,13 +60,15 @@ fission. ```` Element -------------------- -The ```` element indicates two kinds of cutoffs. The first is the weight -cutoff used below which particles undergo Russian roulette. Surviving particles -are assigned a user-determined weight. Note that weight cutoffs and Russian -rouletting are not turned on by default. The second is the energy cutoff which -is used to kill particles under certain energy. The energy cutoff should not be -used unless you know particles under the energy are of no importance to results -you care. This element has the following attributes/sub-elements: +The ```` element indicates three kinds of cutoffs. The first is the +weight cutoff used below which particles undergo Russian roulette. Surviving +particles are assigned a user-determined weight. Note that weight cutoffs and +Russian rouletting are not turned on by default. The second is the energy cutoff +which is used to kill particles under certain energy. The energy cutoff should +not be used unless you know particles under the energy are of no importance to +results you care. The third is the time cutoff used to kill particles whose time +exceeds a specific cutoff. Particles will be killed exactly at the specified +time. :weight: The weight below which particles undergo Russian roulette. @@ -99,6 +101,26 @@ you care. This element has the following attributes/sub-elements: *Default*: 0.0 + :time_neutron + The time above which neutrons will be killed. + + *Default*: Infinity + + :time_photon + The time above which photons will be killed. + + *Default*: Infinity + + :time_electron + The time above which electrons will be killed. + + *Default*: Infinity + + :time_positron + The time above which positorns will be killed. + + *Default*: Infinity + ---------------------------- ```` ---------------------------- @@ -345,6 +367,15 @@ either "false" or "true". *Default*: false +----------------------- +```` Element +----------------------- + +The ```` element is used to set the seed for the pseudorandom number +generator during generation of colors in plots. + + *Default*: 1 + --------------------- ```` Element --------------------- @@ -1175,4 +1206,3 @@ mesh-based weight windows. The ``weight_windows_file`` element has no attributes and contains the path to a weight windows HDF5 file to load during simulation initialization. - diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 58dbe7f0ec0..69089d14d33 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -225,7 +225,7 @@ distributed over a cube centered at the origin with an edge length of 10 cm, and emitting a pulse of particles from 0 to 10 µs, one would run:: - source = openmc.Source() + source = openmc.IndependentSource() source.space = openmc.stats.Box((-5, -5, -5), (5, 5, 5)) source.angle = openmc.stats.Isotropic() source.energy = openmc.stats.Discrete([10.0e6], [1.0]) @@ -237,11 +237,11 @@ that indicates the relative strength of a source distribution if multiple are used. For example, to create two sources, one that should be sampled 70% of the time and another that should be sampled 30% of the time:: - src1 = openmc.Source() + src1 = openmc.IndependentSource() src1.strength = 0.7 ... - src2 = openmc.Source() + src2 = openmc.IndependentSource() src2.strength = 0.3 ... @@ -251,7 +251,7 @@ Finally, the :attr:`Source.particle` attribute can be used to indicate the source should be composed of particles other than neutrons. For example, the following would generate a photon source:: - source = openmc.Source() + source = openmc.IndependentSource() source.particle = 'photon' ... @@ -266,7 +266,7 @@ File-based Sources OpenMC can use a pregenerated HDF5 source file by specifying the ``filename`` argument to :class:`openmc.Source`:: - settings.source = openmc.Source(filename='source.h5') + settings.source = openmc.IndependentSource(filename='source.h5') Statepoint and source files are generated automatically when a simulation is run and can be used as the starting source in a new simulation. Alternatively, a diff --git a/docs/source/usersguide/troubleshoot.rst b/docs/source/usersguide/troubleshoot.rst index b1cd80f3e65..f86ed727405 100644 --- a/docs/source/usersguide/troubleshoot.rst +++ b/docs/source/usersguide/troubleshoot.rst @@ -28,7 +28,7 @@ commands: .. code-block:: sh mkdir build-debug && cd build-debug - cmake -Ddebug=on /path/to/openmc + cmake -DCMAKE_BUILD_TYPE=Debug /path/to/openmc make Now when you re-run your problem, it should report exactly where the program diff --git a/include/openmc/capi.h b/include/openmc/capi.h index ac69c652fce..a444814f585 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -191,6 +191,16 @@ int openmc_weight_windows_get_bounds(int32_t index, const double** lower_bounds, const double** upper_bounds, size_t* size); int openmc_weight_windows_set_bounds(int32_t index, const double* lower_bounds, const double* upper_bounds, size_t size); +int openmc_weight_windows_get_survival_ratio(int32_t index, double* ratio); +int openmc_weight_windows_set_survival_ratio(int32_t index, double ratio); +int openmc_weight_windows_get_max_lower_bound_ratio( + int32_t index, double* lb_ratio); +int openmc_weight_windows_set_max_lower_bound_ratio( + int32_t index, double lb_ratio); +int openmc_weight_windows_get_weight_cutoff(int32_t index, double* cutoff); +int openmc_weight_windows_set_weight_cutoff(int32_t index, double cutoff); +int openmc_weight_windows_get_max_split(int32_t index, int* max_split); +int openmc_weight_windows_set_max_split(int32_t index, int max_split); size_t openmc_weight_windows_size(); int openmc_weight_windows_export(const char* filename = nullptr); int openmc_weight_windows_import(const char* filename = nullptr); diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index 2da3c836398..d77c34bd5f5 100644 --- a/include/openmc/mgxs.h +++ b/include/openmc/mgxs.h @@ -16,20 +16,6 @@ namespace openmc { -//============================================================================== -// Cache contains the cached data for an MGXS object -//============================================================================== - -struct CacheData { - double sqrtkT; // last temperature corresponding to t - int t; // temperature index - int a; // angle index - // last angle that corresponds to a - double u; - double v; - double w; -}; - //============================================================================== // MGXS contains the mgxs data for a nuclide/material //============================================================================== @@ -96,10 +82,9 @@ class Mgxs { bool equiv(const Mgxs& that); public: - std::string name; // name of dataset, e.g., UO2 - double awr; // atomic weight ratio - bool fissionable; // Is this fissionable - vector cache; // index and data cache + std::string name; // name of dataset, e.g., UO2 + double awr; // atomic weight ratio + bool fissionable; // Is this fissionable Mgxs() = default; @@ -135,13 +120,15 @@ class Mgxs { //! @param mu Cosine of the change-in-angle, for scattering quantities; //! use nullptr if irrelevant. //! @param dg delayed group index; use nullptr if irrelevant. + //! @param t Temperature index. + //! @param a Angle index. //! @return Requested cross section value. - double get_xs( - MgxsType xstype, int gin, const int* gout, const double* mu, const int* dg); + double get_xs(MgxsType xstype, int gin, const int* gout, const double* mu, + const int* dg, int t, int a); - inline double get_xs(MgxsType xstype, int gin) + inline double get_xs(MgxsType xstype, int gin, int t, int a) { - return get_xs(xstype, gin, nullptr, nullptr, nullptr); + return get_xs(xstype, gin, nullptr, nullptr, nullptr, t, a); } //! \brief Samples the fission neutron energy and if prompt or delayed. @@ -150,7 +137,10 @@ class Mgxs { //! @param dg Sampled delayed group index. //! @param gout Sampled outgoing energy group. //! @param seed Pseudorandom seed pointer - void sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed); + //! @param t Temperature index. + //! @param a Angle index. + void sample_fission_energy( + int gin, int& dg, int& gout, uint64_t* seed, int t, int a); //! \brief Samples the outgoing energy and angle from a scatter event. //! @@ -159,23 +149,37 @@ class Mgxs { //! @param mu Sampled cosine of the change-in-angle. //! @param wgt Weight of the particle to be adjusted. //! @param seed Pseudorandom seed pointer. + //! @param t Temperature index. + //! @param a Angle index. void sample_scatter( - int gin, int& gout, double& mu, double& wgt, uint64_t* seed); + int gin, int& gout, double& mu, double& wgt, uint64_t* seed, int t, int a); //! \brief Calculates cross section quantities needed for tracking. //! //! @param p The particle whose attributes set which MGXS to get. void calculate_xs(Particle& p); - //! \brief Sets the temperature index in cache given a temperature + //! \brief Sets the temperature index in the particle's cache. + //! + //! @param p Particle. + void set_temperature_index(Particle& p); + + //! \brief Gets the temperature index given a temperature. //! //! @param sqrtkT Temperature of the material. - void set_temperature_index(double sqrtkT); + //! @return The temperature index corresponding to sqrtkT. + int get_temperature_index(double sqrtkT) const; + + //! \brief Sets the angle index in the particle's cache. + //! + //! @param p Particle. + void set_angle_index(Particle& p); - //! \brief Sets the angle index in cache given a direction + //! \brief Gets the angle index given a direction. //! //! @param u Incoming particle direction. - void set_angle_index(Direction u); + //! @return The angle index corresponding to u. + int get_angle_index(const Direction& u) const; //! \brief Provide const access to list of XsData held by this const vector& get_xsdata() const { return xs; } diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 376c9abf305..80a561f4d05 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -40,6 +40,10 @@ class Particle : public ParticleData { double speed() const; + //! moves the particle by the distance length to its next location + //! \param length the distance the particle is moved + void move_distance(double length); + //! create a secondary particle // //! stores the current phase space attributes of the particle in the diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 6d86d528e70..d0813f9c5df 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -170,6 +170,18 @@ struct MacroXS { double pair_production; //!< macroscopic pair production xs }; +//============================================================================== +// Cache contains the cached data for an MGXS object +//============================================================================== + +struct CacheDataMG { + int material {-1}; //!< material index + double sqrtkT; //!< last temperature corresponding to t + int t {0}; //!< temperature index + int a {0}; //!< angle index + Direction u; //!< angle that corresponds to a +}; + //============================================================================== // Information about nearest boundary crossing //============================================================================== @@ -231,6 +243,7 @@ class ParticleData { vector neutron_xs_; //!< Microscopic neutron cross sections vector photon_xs_; //!< Microscopic photon cross sections MacroXS macro_xs_; //!< Macroscopic cross sections + CacheDataMG mg_xs_cache_; //!< Multigroup XS cache int64_t id_; //!< Unique ID ParticleType type_ {ParticleType::neutron}; //!< Particle type (n, p, e, etc.) @@ -349,6 +362,8 @@ class ParticleData { ElementMicroXS& photon_xs(int i) { return photon_xs_[i]; } MacroXS& macro_xs() { return macro_xs_; } const MacroXS& macro_xs() const { return macro_xs_; } + CacheDataMG& mg_xs_cache() { return mg_xs_cache_; } + const CacheDataMG& mg_xs_cache() const { return mg_xs_cache_; } int64_t& id() { return id_; } const int64_t& id() const { return id_; } diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 3ce13c89f90..9f7882643ae 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -92,6 +92,8 @@ extern ElectronTreatment electron_treatment; //!< how to treat secondary electrons extern array energy_cutoff; //!< Energy cutoff in [eV] for each particle type +extern array + time_cutoff; //!< Time cutoff in [s] for each particle type extern int legendre_to_tabular_points; //!< number of points to convert Legendres extern int max_order; //!< Maximum Legendre order for multigroup data diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 992c8fdc155..9852fe679d0 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -155,6 +155,22 @@ class WeightWindows { void set_particle_type(ParticleType p_type); + double survival_ratio() const { return survival_ratio_; } + + double& survival_ratio() { return survival_ratio_; } + + double max_lower_bound_ratio() const { return max_lb_ratio_; } + + double& max_lower_bound_ratio() { return max_lb_ratio_; } + + int max_split() const { return max_split_; } + + int& max_split() { return max_split_; } + + double weight_cutoff() const { return weight_cutoff_; } + + double& weight_cutoff() { return weight_cutoff_; } + //---------------------------------------------------------------------------- // Accessors int32_t id() const { return id_; } diff --git a/openmc/cell.py b/openmc/cell.py index b0be1ea79ae..8572a640d09 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -1,4 +1,3 @@ -from collections import OrderedDict from collections.abc import Iterable from math import cos, sin, pi from numbers import Real @@ -86,7 +85,7 @@ class Cell(IDManagerMixin): calculated in a stochastic volume calculation and added via the :meth:`Cell.add_volume_information` method. For 'distribmat' cells it is the total volume of all instances. - atoms : collections.OrderedDict + atoms : dict Mapping of nuclides to the total number of atoms for each nuclide present in the cell, or in all of its instances for a 'distribmat' fill. For example, {'U235': 1.0e22, 'U238': 5.0e22, ...}. @@ -316,7 +315,7 @@ def atoms(self): # Assumes that volume is total volume of all instances # Also assumes that all instances have the same volume partial_volume = self.volume / len(self.fill) - self._atoms = OrderedDict() + self._atoms = {} for mat in self.fill: for key, atom_per_bcm in mat.get_nuclide_atom_densities().items(): # To account for overlap of nuclides between distribmat @@ -389,13 +388,13 @@ def get_nuclide_densities(self): Returns ------- - nuclides : collections.OrderedDict + nuclides : dict Dictionary whose keys are nuclide names and values are 2-tuples of (nuclide, density) """ - nuclides = OrderedDict() + nuclides = {} if self.fill_type == 'material': nuclides.update(self.fill.get_nuclide_densities()) @@ -423,13 +422,13 @@ def get_all_cells(self, memo=None): Returns ------- - cells : collections.orderedDict + cells : dict Dictionary whose keys are cell IDs and values are :class:`Cell` instances """ - cells = OrderedDict() + cells = {} if memo and self in memo: return cells @@ -447,12 +446,12 @@ def get_all_materials(self, memo=None): Returns ------- - materials : collections.OrderedDict + materials : dict Dictionary whose keys are material IDs and values are :class:`Material` instances """ - materials = OrderedDict() + materials = {} if self.fill_type == 'material': materials[self.fill.id] = self.fill elif self.fill_type == 'distribmat': @@ -473,13 +472,13 @@ def get_all_universes(self): Returns ------- - universes : collections.OrderedDict + universes : dict Dictionary whose keys are universe IDs and values are :class:`Universe` instances """ - universes = OrderedDict() + universes = {} if self.fill_type == 'universe': universes[self.fill.id] = self.fill @@ -560,6 +559,67 @@ def clone(self, clone_materials=True, clone_regions=True, memo=None): return memo[self] + def plot(self, *args, **kwargs): + """Display a slice plot of the cell. + + .. versionadded:: 0.13.4 + + Parameters + ---------- + origin : iterable of float + Coordinates at the origin of the plot. If left as None then the + bounding box center will be used to attempt to ascertain the origin. + Defaults to (0, 0, 0) if the bounding box is not finite + width : iterable of float + Width of the plot in each basis direction. If left as none then the + bounding box width will be used to attempt to ascertain the plot + width. Defaults to (10, 10) if the bounding box is not finite + pixels : Iterable of int or int + If iterable of ints provided, then this directly sets the number of + pixels to use in each basis direction. If int provided, then this + sets the total number of pixels in the plot and the number of pixels + in each basis direction is calculated from this total and the image + aspect ratio. + basis : {'xy', 'xz', 'yz'} + The basis directions for the plot + color_by : {'cell', 'material'} + Indicate whether the plot should be colored by cell or by material + colors : dict + Assigns colors to specific materials or cells. Keys are instances of + :class:`Cell` or :class:`Material` and values are RGB 3-tuples, RGBA + 4-tuples, or strings indicating SVG color names. Red, green, blue, + and alpha should all be floats in the range [0.0, 1.0], for example: + + .. code-block:: python + + # Make water blue + water = openmc.Cell(fill=h2o) + universe.plot(..., colors={water: (0., 0., 1.)) + seed : int + Seed for the random number generator + openmc_exec : str + Path to OpenMC executable. + axes : matplotlib.Axes + Axes to draw to + legend : bool + Whether a legend showing material or cell names should be drawn + legend_kwargs : dict + Keyword arguments passed to :func:`matplotlib.pyplot.legend`. + outline : bool + Whether outlines between color boundaries should be drawn + axis_units : {'km', 'm', 'cm', 'mm'} + Units used on the plot axis + **kwargs + Keyword arguments passed to :func:`matplotlib.pyplot.imshow` + + Returns + ------- + matplotlib.axes.Axes + Axes containing resulting image + + """ + return openmc.Universe(cells=[self]).plot(*args, **kwargs) + def create_xml_subelement(self, xml_element, memo=None): """Add the cell's xml representation to an incoming xml element diff --git a/openmc/data/ace.py b/openmc/data/ace.py index 91cbe41960d..7c348bd176c 100644 --- a/openmc/data/ace.py +++ b/openmc/data/ace.py @@ -15,7 +15,6 @@ """ -from collections import OrderedDict import enum from pathlib import Path import struct @@ -544,7 +543,7 @@ def get_libraries_from_xsdir(path): # Create list of ACE libraries -- we use an ordered dictionary while # building to get O(1) membership checks while retaining insertion order - libraries = OrderedDict() + libraries = {} for line in lines: words = line.split() if len(words) < 3: @@ -576,7 +575,7 @@ def get_libraries_from_xsdata(path): with open(xsdata, 'r') as xsdata_file: # As in get_libraries_from_xsdir, we use a dict for O(1) membership # check while retaining insertion order - libraries = OrderedDict() + libraries = {} for line in xsdata_file: words = line.split() if len(words) >= 9: diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 11d1456f853..e0c0032f072 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -1,4 +1,3 @@ -from collections import OrderedDict from collections.abc import Mapping, MutableMapping from io import StringIO from math import log10 @@ -76,7 +75,7 @@ class IncidentNeutron(EqualityMixin): state. name : str Name of the nuclide using the GNDS naming convention - reactions : collections.OrderedDict + reactions : dict Contains the cross sections, secondary angle and energy distributions, and other associated data for each reaction. The keys are the MT values and the values are Reaction objects. @@ -107,7 +106,7 @@ def __init__(self, name, atomic_number, mass_number, metastable, self.kTs = kTs self.energy = {} self._fission_energy = None - self.reactions = OrderedDict() + self.reactions = {} self._urr = {} self._resonances = None diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 1d2546f0e55..7c96e9fac0d 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -1,4 +1,3 @@ -from collections import OrderedDict from collections.abc import Mapping, Callable from copy import deepcopy from io import StringIO @@ -432,7 +431,7 @@ class IncidentPhoton(EqualityMixin): the projection of the electron momentum on the scattering vector, :math:`p_z` for each subshell). Note that subshell occupancies may not match the atomic relaxation data. - reactions : collections.OrderedDict + reactions : dict Contains the cross sections for each photon reaction. The keys are MT values and the values are instances of :class:`PhotonReaction`. @@ -441,7 +440,7 @@ class IncidentPhoton(EqualityMixin): def __init__(self, atomic_number): self.atomic_number = atomic_number self._atomic_relaxation = None - self.reactions = OrderedDict() + self.reactions = {} self.compton_profiles = {} self.bremsstrahlung = {} diff --git a/openmc/deplete/atom_number.py b/openmc/deplete/atom_number.py index 78ceecca66a..b24c048cc92 100644 --- a/openmc/deplete/atom_number.py +++ b/openmc/deplete/atom_number.py @@ -2,8 +2,6 @@ An ndarray to store atom densities with string, integer, or slice indexing. """ -from collections import OrderedDict - import numpy as np from openmc import Material @@ -47,8 +45,8 @@ class AtomNumber: """ def __init__(self, local_mats, nuclides, volume, n_nuc_burn): - self.index_mat = OrderedDict((mat, i) for i, mat in enumerate(local_mats)) - self.index_nuc = OrderedDict((nuc, i) for i, nuc in enumerate(nuclides)) + self.index_mat = {mat: i for i, mat in enumerate(local_mats)} + self.index_nuc = {nuc: i for i, nuc in enumerate(nuclides)} self.volume = np.ones(len(local_mats)) for mat, val in volume.items(): diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 51fb0ad9cf1..5c80b1842c4 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -9,7 +9,7 @@ import math import os import re -from collections import OrderedDict, defaultdict, namedtuple +from collections import defaultdict, namedtuple from collections.abc import Mapping, Iterable from numbers import Real, Integral from warnings import warn @@ -248,7 +248,7 @@ class Chain: Nuclides present in the chain. reactions : list of str Reactions that are tracked in the depletion chain - nuclide_dict : OrderedDict of str to int + nuclide_dict : dict of str to int Maps a nuclide name to an index in nuclides. fission_yields : None or iterable of dict List of effective fission yields for materials. Each dictionary @@ -264,7 +264,7 @@ class Chain: def __init__(self): self.nuclides = [] self.reactions = [] - self.nuclide_dict = OrderedDict() + self.nuclide_dict = {} self._fission_yields = None def __contains__(self, nuclide): diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index 5fc4ff81a13..daadafe2059 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -6,7 +6,6 @@ """ from abc import abstractmethod -from collections import OrderedDict from warnings import warn import numpy as np @@ -178,7 +177,7 @@ def _get_burnable_mats(self): ------- burnable_mats : list of str List of burnable material IDs - volume : OrderedDict of str to float + volume : dict of str to float Volume of each material in [cm^3] nuclides : list of str Nuclides in order of how they'll appear in the simulation. @@ -187,7 +186,7 @@ def _get_burnable_mats(self): burnable_mats = set() model_nuclides = set() - volume = OrderedDict() + volume = {} self.heavy_metal = 0.0 @@ -258,7 +257,7 @@ def _extract_number(self, local_mats, volume, all_nuclides, prev_res=None): ---------- local_mats : list of str Material IDs to be managed by this process - volume : OrderedDict of str to float + volume : dict of str to float Volumes for the above materials in [cm^3] all_nuclides : list of str Nuclides to be used in the simulation. diff --git a/openmc/deplete/reaction_rates.py b/openmc/deplete/reaction_rates.py index 299f1133f62..b806dfad84b 100644 --- a/openmc/deplete/reaction_rates.py +++ b/openmc/deplete/reaction_rates.py @@ -29,11 +29,11 @@ class ReactionRates(np.ndarray): Attributes ---------- - index_mat : OrderedDict of str to int + index_mat : dict of str to int A dictionary mapping material ID as string to index. - index_nuc : OrderedDict of str to int + index_nuc : dict of str to int A dictionary mapping nuclide name as string to index. - index_rx : OrderedDict of str to int + index_rx : dict of str to int A dictionary mapping reaction name as string to index. n_mat : int Number of materials. diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index f250ab36397..d205bf5718e 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -2,7 +2,7 @@ import bisect import math import typing # required to prevent typing.Union namespace overwriting Union -from typing import Iterable, Optional, Tuple +from typing import Iterable, Optional, Tuple, List from warnings import warn import h5py @@ -95,6 +95,59 @@ def from_hdf5(cls, filename: PathLike): ) return cls(filename) + def get_activity( + self, + mat: typing.Union[Material, str], + units: str = "Bq/cm3", + by_nuclide: bool = False, + volume: Optional[float] = None + ) -> Tuple[np.ndarray, typing.Union[np.ndarray, List[dict]]]: + """Get activity of material over time. + + Parameters + ---------- + mat : openmc.Material, str + Material object or material id to evaluate + units : {'Bq', 'Bq/g', 'Bq/cm3'} + Specifies the type of activity to return, options include total + activity [Bq], specific [Bq/g] or volumetric activity [Bq/cm3]. + by_nuclide : bool + Specifies if the activity should be returned for the material as a + whole or per nuclide. Default is False. + volume : float, optional + Volume of the material. If not passed, defaults to using the + :attr:`Material.volume` attribute. + + Returns + ------- + times : numpy.ndarray + Array of times in [s] + activities : numpy.ndarray or List[dict] + Array of total activities if by_nuclide = False (default) + or list of dictionaries of activities by nuclide if + by_nuclide = True. + + """ + if isinstance(mat, Material): + mat_id = str(mat.id) + elif isinstance(mat, str): + mat_id = mat + else: + raise TypeError('mat should be of type openmc.Material or str') + + times = np.empty_like(self, dtype=float) + if by_nuclide: + activities = [None] * len(self) + else: + activities = np.empty_like(self, dtype=float) + + # Evaluate activity for each depletion time + for i, result in enumerate(self): + times[i] = result.time[0] + activities[i] = result.get_material(mat_id).get_activity(units, by_nuclide, volume) + + return times, activities + def get_atoms( self, mat: typing.Union[Material, str], @@ -158,6 +211,60 @@ def get_atoms( return times, concentrations + def get_decay_heat( + self, + mat: typing.Union[Material, str], + units: str = "W", + by_nuclide: bool = False, + volume: Optional[float] = None + ) -> Tuple[np.ndarray, typing.Union[np.ndarray, List[dict]]]: + """Get decay heat of material over time. + + Parameters + ---------- + mat : openmc.Material, str + Material object or material id to evaluate. + units : {'W', 'W/g', 'W/cm3'} + Specifies the units of decay heat to return. Options include total + heat [W], specific [W/g] or volumetric heat [W/cm3]. + by_nuclide : bool + Specifies if the decay heat should be returned for the material as a + whole or per nuclide. Default is False. + volume : float, optional + Volume of the material. If not passed, defaults to using the + :attr:`Material.volume` attribute. + + Returns + ------- + times : numpy.ndarray + Array of times in [s] + decay_heat : numpy.ndarray or List[dict] + Array of total decay heat values if by_nuclide = False (default) + or list of dictionaries of decay heat values by nuclide if + by_nuclide = True. + """ + + if isinstance(mat, Material): + mat_id = str(mat.id) + elif isinstance(mat, str): + mat_id = mat + else: + raise TypeError('mat should be of type openmc.Material or str') + + times = np.empty_like(self, dtype=float) + if by_nuclide: + decay_heat = [None] * len(self) + else: + decay_heat = np.empty_like(self, dtype=float) + + # Evaluate decay heat for each depletion time + for i, result in enumerate(self): + times[i] = result.time[0] + decay_heat[i] = result.get_material(mat_id).get_decay_heat( + units, by_nuclide, volume) + + return times, decay_heat + def get_mass(self, mat: typing.Union[Material, str], nuc: str, @@ -405,10 +512,11 @@ def get_step_where( if math.isclose(time, times[ix], rel_tol=rtol, abs_tol=atol): return ix + closest = min(times, key=lambda t: abs(time - t)) raise ValueError( - "A value of {} {} was not found given absolute and " - "relative tolerances {} and {}.".format( - time, time_units, atol, rtol) + f"A value of {time} {time_units} was not found given absolute and " + f"relative tolerances {atol} and {rtol}. Closest time is {closest} " + f"{time_units}." ) def export_to_materials( diff --git a/openmc/deplete/stepresult.py b/openmc/deplete/stepresult.py index 689627bbfdf..9af6b5fee76 100644 --- a/openmc/deplete/stepresult.py +++ b/openmc/deplete/stepresult.py @@ -6,7 +6,6 @@ import copy import warnings -from collections import OrderedDict from pathlib import Path import h5py @@ -43,13 +42,13 @@ class StepResult: Number of nuclides. rates : list of ReactionRates The reaction rates for each substep. - volume : OrderedDict of str to float + volume : dict of str to float Dictionary mapping mat id to volume. - index_mat : OrderedDict of str to int + index_mat : dict of str to int A dictionary mapping mat ID as string to index. - index_nuc : OrderedDict of str to int + index_nuc : dict of str to int A dictionary mapping nuclide name as string to index. - mat_to_hdf5_ind : OrderedDict of str to int + mat_to_hdf5_ind : dict of str to int A dictionary mapping mat ID as string to global index. n_hdf5_mats : int Number of materials in entire geometry. @@ -462,11 +461,11 @@ def from_hdf5(cls, handle, step): results.proc_time = np.array([np.nan]) # Reconstruct dictionaries - results.volume = OrderedDict() - results.index_mat = OrderedDict() - results.index_nuc = OrderedDict() - rxn_nuc_to_ind = OrderedDict() - rxn_to_ind = OrderedDict() + results.volume = {} + results.index_mat = {} + results.index_nuc = {} + rxn_nuc_to_ind = {} + rxn_to_ind = {} for mat, mat_handle in handle["/materials"].items(): vol = mat_handle.attrs["volume"] diff --git a/openmc/element.py b/openmc/element.py index bc3195aea4b..01651bebb6f 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -1,4 +1,3 @@ -from collections import OrderedDict import re import lxml.etree as ET @@ -124,7 +123,7 @@ def expand(self, percent, percent_type, enrichment=None, natural_nuclides = {name for name, abundance in natural_isotopes(self)} # Create dict to store the expanded nuclides and abundances - abundances = OrderedDict() + abundances = {} # If cross_sections is None, get the cross sections from the global # configuration diff --git a/openmc/filter.py b/openmc/filter.py index 9b3dd516291..d8cf59c0073 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1,5 +1,4 @@ from abc import ABCMeta -from collections import OrderedDict from collections.abc import Iterable import hashlib from itertools import product @@ -1677,7 +1676,7 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs): level_key = f'level {i_level + 1}' # Create a dictionary for this level for Pandas Multi-index - level_dict = OrderedDict() + level_dict = {} # Use the first distribcell path to determine if level # is a universe/cell or lattice level diff --git a/openmc/geometry.py b/openmc/geometry.py index d4e7aa201d0..82a3fe61c5c 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -1,7 +1,7 @@ from __future__ import annotations import os import typing -from collections import OrderedDict, defaultdict +from collections import defaultdict from copy import deepcopy from collections.abc import Iterable from pathlib import Path @@ -39,11 +39,16 @@ class Geometry: """ - def __init__(self, root=None): + def __init__( + self, + root: typing.Optional[openmc.UniverseBase] = None, + merge_surfaces: bool = False, + surface_precision: int = 10 + ): self._root_universe = None self._offsets = {} - self.merge_surfaces = False - self.surface_precision = 10 + self.merge_surfaces = merge_surfaces + self.surface_precision = surface_precision if root is not None: if isinstance(root, openmc.UniverseBase): self.root_universe = root @@ -357,41 +362,41 @@ def get_instances(self, paths) -> typing.Union[int, typing.List[int]]: return indices if return_list else indices[0] - def get_all_cells(self) -> typing.OrderedDict[int, openmc.Cell]: + def get_all_cells(self) -> typing.Dict[int, openmc.Cell]: """Return all cells in the geometry. Returns ------- - collections.OrderedDict + dict Dictionary mapping cell IDs to :class:`openmc.Cell` instances """ if self.root_universe is not None: return self.root_universe.get_all_cells(memo=set()) else: - return OrderedDict() + return {} - def get_all_universes(self) -> typing.OrderedDict[int, openmc.Universe]: + def get_all_universes(self) -> typing.Dict[int, openmc.Universe]: """Return all universes in the geometry. Returns ------- - collections.OrderedDict + dict Dictionary mapping universe IDs to :class:`openmc.Universe` instances """ - universes = OrderedDict() + universes = {} universes[self.root_universe.id] = self.root_universe universes.update(self.root_universe.get_all_universes()) return universes - def get_all_materials(self) -> typing.OrderedDict[int, openmc.Material]: + def get_all_materials(self) -> typing.Dict[int, openmc.Material]: """Return all materials within the geometry. Returns ------- - collections.OrderedDict + dict Dictionary mapping material IDs to :class:`openmc.Material` instances @@ -399,19 +404,19 @@ def get_all_materials(self) -> typing.OrderedDict[int, openmc.Material]: if self.root_universe is not None: return self.root_universe.get_all_materials(memo=set()) else: - return OrderedDict() + return {} - def get_all_material_cells(self) -> typing.OrderedDict[int, openmc.Cell]: + def get_all_material_cells(self) -> typing.Dict[int, openmc.Cell]: """Return all cells filled by a material Returns ------- - collections.OrderedDict + dict Dictionary mapping cell IDs to :class:`openmc.Cell` instances that are filled with materials or distributed materials. """ - material_cells = OrderedDict() + material_cells = {} for cell in self.get_all_cells().values(): if cell.fill_type in ('material', 'distribmat'): @@ -420,7 +425,7 @@ def get_all_material_cells(self) -> typing.OrderedDict[int, openmc.Cell]: return material_cells - def get_all_material_universes(self) -> typing.OrderedDict[int, openmc.Universe]: + def get_all_material_universes(self) -> typing.Dict[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 @@ -428,12 +433,12 @@ def get_all_material_universes(self) -> typing.OrderedDict[int, openmc.Universe] Returns ------- - collections.OrderedDict + dict Dictionary mapping universe IDs to :class:`openmc.Universe` instances with at least one material-filled cell """ - material_universes = OrderedDict() + material_universes = {} for universe in self.get_all_universes().values(): for cell in universe.cells.values(): @@ -443,16 +448,16 @@ def get_all_material_universes(self) -> typing.OrderedDict[int, openmc.Universe] return material_universes - def get_all_lattices(self) -> typing.OrderedDict[int, openmc.Lattice]: + def get_all_lattices(self) -> typing.Dict[int, openmc.Lattice]: """Return all lattices defined Returns ------- - collections.OrderedDict + dict Dictionary mapping lattice IDs to :class:`openmc.Lattice` instances """ - lattices = OrderedDict() + lattices = {} for cell in self.get_all_cells().values(): if cell.fill_type == 'lattice': @@ -461,17 +466,17 @@ def get_all_lattices(self) -> typing.OrderedDict[int, openmc.Lattice]: return lattices - def get_all_surfaces(self) -> typing.OrderedDict[int, openmc.Surface]: + def get_all_surfaces(self) -> typing.Dict[int, openmc.Surface]: """ Return all surfaces used in the geometry Returns ------- - collections.OrderedDict + dict Dictionary mapping surface IDs to :class:`openmc.Surface` instances """ - surfaces = OrderedDict() + surfaces = {} for cell in self.get_all_cells().values(): if cell.region is not None: @@ -729,3 +734,61 @@ def clone(self) -> Geometry: clone = deepcopy(self) clone.root_universe = self.root_universe.clone() return clone + + def plot(self, *args, **kwargs): + """Display a slice plot of the geometry. + + .. versionadded:: 0.13.4 + + Parameters + ---------- + origin : iterable of float + Coordinates at the origin of the plot. If left as None then the + bounding box center will be used to attempt to ascertain the origin. + Defaults to (0, 0, 0) if the bounding box is not finite + width : iterable of float + Width of the plot in each basis direction. If left as none then the + bounding box width will be used to attempt to ascertain the plot + width. Defaults to (10, 10) if the bounding box is not finite + pixels : Iterable of int or int + If iterable of ints provided, then this directly sets the number of + pixels to use in each basis direction. If int provided, then this + sets the total number of pixels in the plot and the number of pixels + in each basis direction is calculated from this total and the image + aspect ratio. + basis : {'xy', 'xz', 'yz'} + The basis directions for the plot + color_by : {'cell', 'material'} + Indicate whether the plot should be colored by cell or by material + colors : dict + Assigns colors to specific materials or cells. Keys are instances of + :class:`Cell` or :class:`Material` and values are RGB 3-tuples, RGBA + 4-tuples, or strings indicating SVG color names. Red, green, blue, + and alpha should all be floats in the range [0.0, 1.0], for example: + .. code-block:: python + # Make water blue + water = openmc.Cell(fill=h2o) + universe.plot(..., colors={water: (0., 0., 1.)) + seed : int + Seed for the random number generator + openmc_exec : str + Path to OpenMC executable. + axes : matplotlib.Axes + Axes to draw to + legend : bool + Whether a legend showing material or cell names should be drawn + legend_kwargs : dict + Keyword arguments passed to :func:`matplotlib.pyplot.legend`. + outline : bool + Whether outlines between color boundaries should be drawn + axis_units : {'km', 'm', 'cm', 'mm'} + Units used on the plot axis + **kwargs + Keyword arguments passed to :func:`matplotlib.pyplot.imshow` + Returns + ------- + matplotlib.axes.Axes + Axes containing resulting image + """ + + return self.root_universe.plot(*args, **kwargs) diff --git a/openmc/lattice.py b/openmc/lattice.py index 0983f6d0c7a..926a016efe3 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -1,5 +1,4 @@ from abc import ABC -from collections import OrderedDict from collections.abc import Iterable from copy import deepcopy from math import sqrt, floor @@ -113,13 +112,13 @@ def get_unique_universes(self): Returns ------- - universes : collections.OrderedDict + universes : dict Dictionary whose keys are universe IDs and values are :class:`openmc.UniverseBase` instances """ - univs = OrderedDict() + univs = {} for k in range(len(self._universes)): for j in range(len(self._universes[k])): if isinstance(self._universes[k][j], openmc.UniverseBase): @@ -164,12 +163,12 @@ def get_all_cells(self, memo=None): Returns ------- - cells : collections.OrderedDict + cells : dict Dictionary whose keys are cell IDs and values are :class:`Cell` instances """ - cells = OrderedDict() + cells = {} if memo and self in memo: return cells @@ -189,13 +188,13 @@ def get_all_materials(self, memo=None): Returns ------- - materials : collections.OrderedDict + materials : dict Dictionary whose keys are material IDs and values are :class:`Material` instances """ - materials = OrderedDict() + materials = {} # Append all Cells in each Cell in the Universe to the dictionary cells = self.get_all_cells(memo) @@ -209,7 +208,7 @@ def get_all_universes(self): Returns ------- - universes : collections.OrderedDict + universes : dict Dictionary whose keys are universe IDs and values are :class:`Universe` instances @@ -217,7 +216,7 @@ def get_all_universes(self): # Initialize a dictionary of all Universes contained by the Lattice # in each nested Universe level - all_universes = OrderedDict() + all_universes = {} # Get all unique Universes contained in each of the lattice cells unique_universes = self.get_unique_universes() @@ -1127,7 +1126,7 @@ def num_rings(self): @property def orientation(self): return self._orientation - + @orientation.setter def orientation(self, orientation): cv.check_value('orientation', orientation.lower(), ('x', 'y')) diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py index d142d6cd208..9f294a42c04 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -229,7 +229,8 @@ def id_map(plot): ------- id_map : numpy.ndarray A NumPy array with shape (vertical pixels, horizontal pixels, 3) of - OpenMC property ids with dtype int32 + OpenMC property ids with dtype int32. The last dimension of the array + contains, in order, cell IDs, cell instances, and material IDs. """ img_data = np.zeros((plot.v_res, plot.h_res, 3), diff --git a/openmc/lib/weight_windows.py b/openmc/lib/weight_windows.py index e881f1d3120..2984eb23318 100644 --- a/openmc/lib/weight_windows.py +++ b/openmc/lib/weight_windows.py @@ -69,6 +69,38 @@ _dll.openmc_weight_windows_get_bounds.restype = c_int _dll.openmc_weight_windows_get_bounds.errcheck = _error_handler +_dll.openmc_weight_windows_get_survival_ratio.argtypes = [c_int32, POINTER(c_double)] +_dll.openmc_weight_windows_get_survival_ratio.restype = c_int +_dll.openmc_weight_windows_get_survival_ratio.errcheck = _error_handler + +_dll.openmc_weight_windows_set_survival_ratio.argtypes = [c_int32, c_double] +_dll.openmc_weight_windows_set_survival_ratio.restype = c_int +_dll.openmc_weight_windows_set_survival_ratio.errcheck = _error_handler + +_dll.openmc_weight_windows_get_max_lower_bound_ratio.argtypes = [c_int32, POINTER(c_double)] +_dll.openmc_weight_windows_get_max_lower_bound_ratio.restype = c_int +_dll.openmc_weight_windows_get_max_lower_bound_ratio.errcheck = _error_handler + +_dll.openmc_weight_windows_set_max_lower_bound_ratio.argtypes = [c_int32, c_double] +_dll.openmc_weight_windows_set_max_lower_bound_ratio.restype = c_int +_dll.openmc_weight_windows_set_max_lower_bound_ratio.errcheck = _error_handler + +_dll.openmc_weight_windows_get_weight_cutoff.argtypes = [c_int32, POINTER(c_double)] +_dll.openmc_weight_windows_get_weight_cutoff.restype = c_int +_dll.openmc_weight_windows_get_weight_cutoff.errcheck = _error_handler + +_dll.openmc_weight_windows_set_weight_cutoff.argtypes = [c_int32, c_double] +_dll.openmc_weight_windows_set_weight_cutoff.restype = c_int +_dll.openmc_weight_windows_set_weight_cutoff.errcheck = _error_handler + +_dll.openmc_weight_windows_get_max_split.argtypes = [c_int32, POINTER(c_int)] +_dll.openmc_weight_windows_get_max_split.restype = c_int +_dll.openmc_weight_windows_get_max_split.errcheck = _error_handler + +_dll.openmc_weight_windows_set_max_split.argtypes = [c_int32, c_int] +_dll.openmc_weight_windows_set_max_split.restype = c_int +_dll.openmc_weight_windows_set_max_split.errcheck = _error_handler + class WeightWindows(_FortranObjectWithID): """WeightWindows stored internally. @@ -201,6 +233,46 @@ def bounds(self, bounds): _dll.openmc_weight_windows_set_bounds(self._index, lower_p, upper_p, lower.size) + @property + def survival_ratio(self): + ratio = c_double() + _dll.openmc_weight_windows_get_survival_ratio(self._index, ratio) + return ratio.value + + @survival_ratio.setter + def survival_ratio(self, ratio): + _dll.openmc_weight_windows_set_survival_ratio(self._index, ratio) + + @property + def max_lower_bound_ratio(self): + lb_ratio = c_double() + _dll.openmc_weight_windows_get_max_lower_bound_ratio(self._index, lb_ratio) + return lb_ratio.value + + @max_lower_bound_ratio.setter + def max_lower_bound_ratio(self, lb_ratio): + _dll.openmc_weight_windows_set_max_lower_bound_ratio(self._index, lb_ratio) + + @property + def weight_cutoff(self): + cutoff = c_double() + _dll.openmc_weight_windows_get_weight_cutoff(self._index, cutoff) + return cutoff.value + + @weight_cutoff.setter + def weight_cutoff(self, cutoff): + _dll.openmc_weight_windows_set_weight_cutoff(self._index, cutoff) + + @property + def max_split(self): + max_split = c_int() + _dll.openmc_weight_windows_get_max_split(self._index, max_split) + return max_split.value + + @max_split.setter + def max_split(self, max_split): + _dll.openmc_weight_windows_set_max_split(self._index, max_split) + def update_magic(self, tally, value='mean', threshold=1.0, ratio=5.0): """Update weight window values using the MAGIC method diff --git a/openmc/material.py b/openmc/material.py index 4e3fdac8159..f0e9117e63d 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -1,5 +1,5 @@ from __future__ import annotations -from collections import OrderedDict, defaultdict, namedtuple, Counter +from collections import defaultdict, namedtuple, Counter from collections.abc import Iterable from copy import deepcopy from numbers import Real @@ -145,6 +145,7 @@ def __repr__(self) -> str: string += f' [{self._density_units}]\n' string += '{: <16}=\t{} [cm^3]\n'.format('\tVolume', self._volume) + string += '{: <16}=\t{}\n'.format('\tDepletable', self._depletable) string += '{: <16}\n'.format('\tS(a,b) Tables') @@ -641,7 +642,8 @@ def remove_macroscopic(self, macroscopic: str): def add_element(self, element: str, percent: float, percent_type: str = 'ao', enrichment: Optional[float] = None, enrichment_target: Optional[str] = None, - enrichment_type: Optional[str] = None): + enrichment_type: Optional[str] = None, + cross_sections: Optional[str] = None): """Add a natural element to the material Parameters @@ -668,6 +670,8 @@ def add_element(self, element: str, percent: float, percent_type: str = 'ao', Default is: 'ao' for two-isotope enrichment; 'wo' for U enrichment .. versionadded:: 0.12 + cross_sections : str, optional + Location of cross_sections.xml file. Notes ----- @@ -746,7 +750,8 @@ def add_element(self, element: str, percent: float, percent_type: str = 'ao', percent_type, enrichment, enrichment_target, - enrichment_type): + enrichment_type, + cross_sections): self.add_nuclide(*nuclide) def add_elements_from_formula(self, formula: str, percent_type: str = 'ao', @@ -937,8 +942,7 @@ def get_nuclide_densities(self) -> Dict[str, tuple]: """ - # keep ordered dictionary for testing purposes - nuclides = OrderedDict() + nuclides = {} for nuclide in self._nuclides: nuclides[nuclide.name] = nuclide @@ -1023,7 +1027,7 @@ def get_nuclide_atom_densities(self, nuclide: Optional[str] = None) -> Dict[str, nuc_densities = density * nuc_densities - nuclides = OrderedDict() + nuclides = {} for n, nuc in enumerate(nucs): if nuclide is None or nuclide == nuc: nuclides[nuc] = nuc_densities[n] @@ -1456,7 +1460,7 @@ def from_xml_element(cls, elem: ET.Element) -> Material: if "cfg" in elem.attrib: cfg = elem.get("cfg") return Material.from_ncrystal(cfg, material_id=mat_id) - + mat = cls(mat_id) mat.name = elem.get('name') diff --git a/openmc/mesh.py b/openmc/mesh.py index f9eaff6d5f8..5313ebb398b 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1,16 +1,19 @@ +import typing +import warnings from abc import ABC, abstractmethod from collections.abc import Iterable -from collections import OrderedDict from math import pi -from numbers import Real, Integral +from numbers import Integral, Real from pathlib import Path -import warnings -import lxml.etree as ET +from typing import Optional, Sequence, Tuple +import h5py +import lxml.etree as ET import numpy as np -import openmc.checkvalue as cv import openmc +import openmc.checkvalue as cv +from openmc.checkvalue import PathLike from ._xml import get_text from .mixin import IDManagerMixin from .surface import _BOUNDARY_TYPES @@ -38,7 +41,7 @@ class MeshBase(IDManagerMixin, ABC): next_id = 1 used_ids = set() - def __init__(self, mesh_id=None, name=''): + def __init__(self, mesh_id: Optional[int] = None, name: str = ''): # Initialize Mesh class attributes self.id = mesh_id self.name = name @@ -48,7 +51,7 @@ def name(self): return self._name @name.setter - def name(self, name): + def name(self, name: str): if name is not None: cv.check_type(f'name for mesh ID="{self._id}"', name, str) self._name = name @@ -68,7 +71,7 @@ def _volume_dim_check(self): 'Volumes cannot be provided.') @classmethod - def from_hdf5(cls, group): + def from_hdf5(cls, group: h5py.Group): """Create mesh from HDF5 group Parameters @@ -98,7 +101,7 @@ def from_hdf5(cls, group): raise ValueError('Unrecognized mesh type: "' + mesh_type + '"') @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generates a mesh from an XML element Parameters @@ -261,10 +264,10 @@ def num_mesh_cells(self): return np.prod(self.dimension) def write_data_to_vtk(self, - filename, - datasets=None, - volume_normalization=True, - curvilinear=False): + filename: PathLike, + datasets: Optional[dict] = None, + volume_normalization: bool = True, + curvilinear: bool = False): """Creates a VTK object of the mesh Parameters @@ -393,9 +396,9 @@ def _insert_point(pnt): else: return result - ### Add all points to the unstructured grid, maintaining a flat list of IDs as we go ### + # Add all points to the unstructured grid, maintaining a flat list of IDs as we go ### - # flat array storind point IDs for a given vertex + # flat array storing point IDs for a given vertex # in the grid point_ids = [] @@ -452,7 +455,7 @@ def _insert_point(pnt): return vtk_grid - def _check_vtk_datasets(self, datasets): + def _check_vtk_datasets(self, datasets: dict): """Perform some basic checks that the datasets are valid for this mesh Parameters @@ -493,7 +496,7 @@ class RegularMesh(StructuredMesh): name : str Name of the mesh dimension : Iterable of int - The number of mesh cells in each direction. + The number of mesh cells in each direction (x, y, z). n_dimension : int Number of mesh dimensions. lower_left : Iterable of float @@ -502,6 +505,9 @@ class RegularMesh(StructuredMesh): upper_right : Iterable of float The upper-right corner of the structured mesh. If only two coordinate are given, it is assumed that the mesh is an x-y mesh. + bounding_box: openmc.BoundingBox + Axis-aligned bounding box of the cell defined by the upper-right and lower- + left coordinates width : Iterable of float The width of mesh cells in each direction. indices : Iterable of tuple @@ -510,7 +516,7 @@ class RegularMesh(StructuredMesh): """ - def __init__(self, mesh_id=None, name=''): + def __init__(self, mesh_id: Optional[int] = None, name: str = ''): super().__init__(mesh_id, name) self._dimension = None @@ -523,7 +529,7 @@ def dimension(self): return tuple(self._dimension) @dimension.setter - def dimension(self, dimension): + def dimension(self, dimension: typing.Iterable[int]): cv.check_type('mesh dimension', dimension, Iterable, Integral) cv.check_length('mesh dimension', dimension, 1, 3) self._dimension = dimension @@ -540,7 +546,7 @@ def lower_left(self): return self._lower_left @lower_left.setter - def lower_left(self, lower_left): + def lower_left(self, lower_left: typing.Iterable[Real]): cv.check_type('mesh lower_left', lower_left, Iterable, Real) cv.check_length('mesh lower_left', lower_left, 1, 3) self._lower_left = lower_left @@ -560,7 +566,7 @@ def upper_right(self): return [l + w * d for l, w, d in zip(ls, ws, dims)] @upper_right.setter - def upper_right(self, upper_right): + def upper_right(self, upper_right: typing.Iterable[Real]): cv.check_type('mesh upper_right', upper_right, Iterable, Real) cv.check_length('mesh upper_right', upper_right, 1, 3) self._upper_right = upper_right @@ -584,7 +590,7 @@ def width(self): return [(u - l) / d for u, l, d in zip(us, ls, dims)] @width.setter - def width(self, width): + def width(self, width: typing.Iterable[Real]): cv.check_type('mesh width', width, Iterable, Real) cv.check_length('mesh width', width, 1, 3) self._width = width @@ -595,7 +601,7 @@ def width(self, width): @property def cartesian_vertices(self): - """Returns vertices in cartesian coordiantes. Identical to ``vertices`` for RegularMesh and RectilinearMesh + """Returns vertices in cartesian coordinates. Identical to ``vertices`` for RegularMesh and RectilinearMesh """ return self.vertices @@ -674,7 +680,7 @@ def __repr__(self): return string @classmethod - def from_hdf5(cls, group): + def from_hdf5(cls, group: h5py.Group): mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) # Read and assign mesh properties @@ -691,7 +697,13 @@ def from_hdf5(cls, group): return mesh @classmethod - def from_rect_lattice(cls, lattice, division=1, mesh_id=None, name=''): + def from_rect_lattice( + cls, + lattice: 'openmc.RectLattice', + division: int = 1, + mesh_id: Optional[int] = None, + name: str = '' + ): """Create mesh from an existing rectangular lattice Parameters @@ -717,7 +729,7 @@ def from_rect_lattice(cls, lattice, division=1, mesh_id=None, name=''): shape = np.array(lattice.shape) width = lattice.pitch*shape - mesh = cls(mesh_id, name) + mesh = cls(mesh_id=mesh_id, name=name) mesh.lower_left = lattice.lower_left mesh.upper_right = lattice.lower_left + width mesh.dimension = shape*division @@ -727,10 +739,10 @@ def from_rect_lattice(cls, lattice, division=1, mesh_id=None, name=''): @classmethod def from_domain( cls, - domain, - dimension=(10, 10, 10), - mesh_id=None, - name='' + domain: typing.Union['openmc.Cell', 'openmc.Region', 'openmc.Universe', 'openmc.Geometry'], + dimension: Sequence[int] = (10, 10, 10), + mesh_id: Optional[int] = None, + name: str = '' ): """Create mesh from an existing openmc cell, region, universe or geometry by making use of the objects bounding box property. @@ -740,9 +752,9 @@ def from_domain( domain : {openmc.Cell, openmc.Region, openmc.Universe, openmc.Geometry} The object passed in will be used as a template for this mesh. The bounding box of the property of the object passed will be used to - set the lower_left and upper_right of the mesh instance + set the lower_left and upper_right and of the mesh instance dimension : Iterable of int - The number of mesh cells in each direction. + The number of mesh cells in each direction (x, y, z). mesh_id : int Unique identifier for the mesh name : str @@ -760,7 +772,7 @@ def from_domain( (openmc.Cell, openmc.Region, openmc.Universe, openmc.Geometry), ) - mesh = cls(mesh_id, name) + mesh = cls(mesh_id=mesh_id, name=name) mesh.lower_left = domain.bounding_box[0] mesh.upper_right = domain.bounding_box[1] mesh.dimension = dimension @@ -797,7 +809,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate mesh from an XML element Parameters @@ -812,7 +824,7 @@ def from_xml_element(cls, elem): """ mesh_id = int(get_text(elem, 'id')) - mesh = cls(mesh_id) + mesh = cls(mesh_id=mesh_id) mesh_type = get_text(elem, 'type') if mesh_type is not None: @@ -836,7 +848,7 @@ def from_xml_element(cls, elem): return mesh - def build_cells(self, bc=None): + def build_cells(self, bc: Optional[str] = None): """Generates a lattice of universes with the same dimensionality as the mesh object. The individual cells/universes produced will not have material definitions applied and so downstream code @@ -961,6 +973,7 @@ def build_cells(self, bc=None): return root_cell, cells + def Mesh(*args, **kwargs): warnings.warn("Mesh has been renamed RegularMesh. Future versions of " "OpenMC will not accept the name Mesh.") @@ -984,7 +997,7 @@ class RectilinearMesh(StructuredMesh): name : str Name of the mesh dimension : Iterable of int - The number of mesh cells in each direction. + The number of mesh cells in each direction (x, y, z). n_dimension : int Number of mesh dimensions (always 3 for a RectilinearMesh). x_grid : numpy.ndarray @@ -999,7 +1012,7 @@ class RectilinearMesh(StructuredMesh): """ - def __init__(self, mesh_id=None, name=''): + def __init__(self, mesh_id: int = None, name: str = ''): super().__init__(mesh_id, name) self._x_grid = None @@ -1106,11 +1119,11 @@ def __repr__(self): return string @classmethod - def from_hdf5(cls, group): + def from_hdf5(cls, group: h5py.Group): mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) # Read and assign mesh properties - mesh = cls(mesh_id) + mesh = cls(mesh_id=mesh_id) mesh.x_grid = group['x_grid'][()] mesh.y_grid = group['y_grid'][()] mesh.z_grid = group['z_grid'][()] @@ -1118,7 +1131,7 @@ def from_hdf5(cls, group): return mesh @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate a rectilinear mesh from an XML element Parameters @@ -1132,8 +1145,8 @@ def from_xml_element(cls, elem): Rectilinear mesh object """ - id = int(get_text(elem, 'id')) - mesh = cls(id) + mesh_id = int(get_text(elem, 'id')) + mesh = cls(mesh_id=mesh_id) mesh.x_grid = [float(x) for x in get_text(elem, 'x_grid').split()] mesh.y_grid = [float(y) for y in get_text(elem, 'y_grid').split()] mesh.z_grid = [float(z) for z in get_text(elem, 'z_grid').split()] @@ -1171,6 +1184,17 @@ class CylindricalMesh(StructuredMesh): Parameters ---------- + r_grid : numpy.ndarray + 1-D array of mesh boundary points along the r-axis. + Requirement is r >= 0. + z_grid : numpy.ndarray + 1-D array of mesh boundary points along the z-axis. + phi_grid : numpy.ndarray + 1-D array of mesh boundary points along the phi-axis in radians. + The default value is [0, 2π], i.e. the full phi range. + origin : numpy.ndarray + 1-D array of length 3 the (x,y,z) origin of the mesh in + cartesian coordinates mesh_id : int Unique identifier for the mesh name : str @@ -1183,7 +1207,8 @@ class CylindricalMesh(StructuredMesh): name : str Name of the mesh dimension : Iterable of int - The number of mesh cells in each direction. + The number of mesh cells in each direction (r_grid, + phi_grid, z_grid). n_dimension : int Number of mesh dimensions (always 3 for a CylindricalMesh). r_grid : numpy.ndarray @@ -1200,16 +1225,33 @@ class CylindricalMesh(StructuredMesh): indices : Iterable of tuple An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, 1), ...] + lower_left : Iterable of float + The lower-left corner of the structured mesh. If only two coordinate + are given, it is assumed that the mesh is an x-y mesh. + upper_right : Iterable of float + The upper-right corner of the structured mesh. If only two coordinate + are given, it is assumed that the mesh is an x-y mesh. + bounding_box: openmc.OpenMC + Axis-aligned cartesian bounding box of cell defined by upper-right and lower- + left coordinates """ - def __init__(self, mesh_id=None, name=''): + def __init__( + self, + r_grid: Sequence[float], + z_grid: Sequence[float], + phi_grid: Sequence[float] = (0, 2*pi), + origin: Sequence[float] = (0., 0., 0.), + mesh_id: Optional[int] = None, + name: str = '', + ): super().__init__(mesh_id, name) - self._r_grid = None - self._phi_grid = [0.0, 2*pi] - self._z_grid = None - self.origin = (0., 0., 0.) + self._r_grid = r_grid + self._phi_grid = phi_grid + self._z_grid = z_grid + self.origin = origin @property def dimension(self): @@ -1252,7 +1294,7 @@ def phi_grid(self, grid): @property def z_grid(self): return self._z_grid - + @z_grid.setter def z_grid(self, grid): cv.check_type('mesh z_grid', grid, Iterable, Real) @@ -1272,6 +1314,27 @@ def indices(self): for p in range(1, np + 1) for r in range(1, nr + 1)) + + @property + def lower_left(self): + return np.array(( + self.origin[0] - self.r_grid[-1], + self.origin[1] - self.r_grid[-1], + self.origin[2] - self.z_grid[-1] + )) + + @property + def upper_right(self): + return np.array(( + self.origin[0] + self.r_grid[-1], + self.origin[1] + self.r_grid[-1], + self.origin[2] + self.z_grid[-1] + )) + + @property + def bounding_box(self): + return openmc.BoundingBox(self.lower_left, self.upper_right) + def __repr__(self): fmt = '{0: <16}{1}{2}\n' string = super().__repr__() @@ -1295,14 +1358,16 @@ def __repr__(self): return string @classmethod - def from_hdf5(cls, group): + def from_hdf5(cls, group: h5py.Group): mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) # Read and assign mesh properties - mesh = cls(mesh_id) - mesh.r_grid = group['r_grid'][()] - mesh.phi_grid = group['phi_grid'][()] - mesh.z_grid = group['z_grid'][()] + mesh = cls( + mesh_id=mesh_id, + r_grid = group['r_grid'][()], + phi_grid = group['phi_grid'][()], + z_grid = group['z_grid'][()], + ) if 'origin' in group: mesh.origin = group['origin'][()] @@ -1311,11 +1376,11 @@ def from_hdf5(cls, group): @classmethod def from_domain( cls, - domain, - dimension=(10, 10, 10), - mesh_id=None, - phi_grid_bounds=(0.0, 2*pi), - name='' + domain: typing.Union['openmc.Cell', 'openmc.Region', 'openmc.Universe', 'openmc.Geometry'], + dimension: Sequence[int] = (10, 10, 10), + mesh_id: Optional[int] = None, + phi_grid_bounds: Sequence[float] = (0.0, 2*pi), + name: str = '' ): """Creates a regular CylindricalMesh from an existing openmc domain. @@ -1338,8 +1403,8 @@ def from_domain( Returns ------- - openmc.RegularMesh - RegularMesh instance + openmc.CylindricalMesh + CylindricalMesh instance """ cv.check_type( @@ -1348,7 +1413,6 @@ def from_domain( (openmc.Cell, openmc.Region, openmc.Universe, openmc.Geometry), ) - mesh = cls(mesh_id, name) # loaded once to avoid reading h5m file repeatedly cached_bb = domain.bounding_box @@ -1360,21 +1424,24 @@ def from_domain( cached_bb[1][1], ] ) - mesh.r_grid = np.linspace( + r_grid = np.linspace( 0, max_bounding_box_radius, num=dimension[0]+1 ) - mesh.phi_grid = np.linspace( + phi_grid = np.linspace( phi_grid_bounds[0], phi_grid_bounds[1], num=dimension[1]+1 ) - mesh.z_grid = np.linspace( + z_grid = np.linspace( cached_bb[0][2], cached_bb[1][2], num=dimension[2]+1 ) + mesh = cls( + r_grid=r_grid, z_grid=z_grid, phi_grid=phi_grid, mesh_id=mesh_id, name=name + ) return mesh @@ -1407,7 +1474,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate a cylindrical mesh from an XML element Parameters @@ -1423,11 +1490,13 @@ def from_xml_element(cls, elem): """ mesh_id = int(get_text(elem, 'id')) - mesh = cls(mesh_id) - mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] - mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] - mesh.z_grid = [float(x) for x in get_text(elem, "z_grid").split()] - mesh.origin = [float(x) for x in get_text(elem, "origin", default=[0., 0., 0.]).split()] + mesh = cls( + r_grid = [float(x) for x in get_text(elem, "r_grid").split()], + phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()], + z_grid = [float(x) for x in get_text(elem, "z_grid").split()], + origin = [float(x) for x in get_text(elem, "origin", default=[0., 0., 0.]).split()], + mesh_id=mesh_id, + ) return mesh @@ -1453,7 +1522,7 @@ def cartesian_vertices(self): return self._convert_to_cartesian(self.vertices, self.origin) @staticmethod - def _convert_to_cartesian(arr, origin): + def _convert_to_cartesian(arr, origin: Sequence[float]): """Converts an array with xyz values in the first dimension (shape (3, ...)) to Cartesian coordinates. """ @@ -1470,6 +1539,18 @@ class SphericalMesh(StructuredMesh): Parameters ---------- + r_grid : numpy.ndarray + 1-D array of mesh boundary points along the r-axis. + Requirement is r >= 0. + phi_grid : numpy.ndarray + 1-D array of mesh boundary points along the phi-axis in radians. + The default value is [0, 2π], i.e. the full phi range. + theta_grid : numpy.ndarray + 1-D array of mesh boundary points along the theta-axis in radians. + The default value is [0, π], i.e. the full theta range. + origin : numpy.ndarray + 1-D array of length 3 the (x,y,z) origin of the mesh in + cartesian coordinates mesh_id : int Unique identifier for the mesh name : str @@ -1482,7 +1563,8 @@ class SphericalMesh(StructuredMesh): name : str Name of the mesh dimension : Iterable of int - The number of mesh cells in each direction. + The number of mesh cells in each direction (r_grid, + theta_grid, phi_grid). n_dimension : int Number of mesh dimensions (always 3 for a SphericalMesh). r_grid : numpy.ndarray @@ -1500,16 +1582,33 @@ class SphericalMesh(StructuredMesh): indices : Iterable of tuple An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, 1), ...] + lower_left : numpy.ndarray + The lower-left corner of the structured mesh. If only two coordinate + are given, it is assumed that the mesh is an x-y mesh. + upper_right : numpy.ndarray + The upper-right corner of the structured mesh. If only two coordinate + are given, it is assumed that the mesh is an x-y mesh. + bounding_box : openmc.BoundingBox + Axis-aligned bounding box of the cell defined by the upper-right and lower- + left coordinates """ - def __init__(self, mesh_id=None, name=''): + def __init__( + self, + r_grid: Sequence[float], + phi_grid: Sequence[float] = (0, 2*pi), + theta_grid: Sequence[float] = (0, pi), + origin: Sequence[float] = (0., 0., 0.), + mesh_id: Optional[int] = None, + name: str = '', + ): super().__init__(mesh_id, name) - self._r_grid = None - self._theta_grid = [0, pi] - self._phi_grid = [0, 2*pi] - self.origin = (0., 0., 0.) + self._r_grid = r_grid + self._theta_grid = theta_grid + self._phi_grid = phi_grid + self.origin = origin @property def dimension(self): @@ -1572,6 +1671,20 @@ def indices(self): for t in range(1, nt + 1) for r in range(1, nr + 1)) + @property + def lower_left(self): + r = self.r_grid[-1] + return np.array((self.origin[0] - r, self.origin[1] - r, self.origin[2] - r)) + + @property + def upper_right(self): + r = self.r_grid[-1] + return np.array((self.origin[0] + r, self.origin[1] + r, self.origin[2] + r)) + + @property + def bounding_box(self): + return openmc.BoundingBox(self.lower_left, self.upper_right) + def __repr__(self): fmt = '{0: <16}{1}{2}\n' string = super().__repr__() @@ -1595,14 +1708,16 @@ def __repr__(self): return string @classmethod - def from_hdf5(cls, group): + def from_hdf5(cls, group: h5py.Group): mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) # Read and assign mesh properties - mesh = cls(mesh_id) - mesh.r_grid = group['r_grid'][()] - mesh.theta_grid = group['theta_grid'][()] - mesh.phi_grid = group['phi_grid'][()] + mesh = cls( + r_grid = group['r_grid'][()], + theta_grid = group['theta_grid'][()], + phi_grid = group['phi_grid'][()], + mesh_id=mesh_id, + ) if 'origin' in group: mesh.origin = group['origin'][()] @@ -1637,7 +1752,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate a spherical mesh from an XML element Parameters @@ -1652,11 +1767,13 @@ def from_xml_element(cls, elem): """ mesh_id = int(get_text(elem, 'id')) - mesh = cls(mesh_id) - mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] - mesh.theta_grid = [float(x) for x in get_text(elem, "theta_grid").split()] - mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] - mesh.origin = [float(x) for x in get_text(elem, "origin", default=[0., 0., 0.]).split()] + mesh = cls( + mesh_id=mesh_id, + r_grid = [float(x) for x in get_text(elem, "r_grid").split()], + theta_grid = [float(x) for x in get_text(elem, "theta_grid").split()], + phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()], + origin = [float(x) for x in get_text(elem, "origin", default=[0., 0., 0.]).split()], + ) return mesh @@ -1682,7 +1799,7 @@ def cartesian_vertices(self): return self._convert_to_cartesian(self.vertices, self.origin) @staticmethod - def _convert_to_cartesian(arr, origin): + def _convert_to_cartesian(arr, origin: Sequence[float]): """Converts an array with xyz values in the first dimension (shape (3, ...)) to Cartesian coordinates. """ @@ -1757,8 +1874,8 @@ class UnstructuredMesh(MeshBase): _LINEAR_TET = 0 _LINEAR_HEX = 1 - def __init__(self, filename, library, mesh_id=None, name='', - length_multiplier=1.0): + def __init__(self, filename: PathLike, library: str, mesh_id: Optional[int] = None, + name: str = '', length_multiplier: float = 1.0): super().__init__(mesh_id, name) self.filename = filename self._volumes = None @@ -1783,7 +1900,7 @@ def library(self): return self._library @library.setter - def library(self, lib): + def library(self, lib: str): cv.check_value('Unstructured mesh library', lib, ('moab', 'libmesh')) self._library = lib @@ -1792,7 +1909,7 @@ def size(self): return self._size @size.setter - def size(self, size): + def size(self, size: int): cv.check_type("Unstructured mesh size", size, Integral) self._size = size @@ -1801,7 +1918,7 @@ def output(self): return self._output @output.setter - def output(self, val): + def output(self, val: bool): cv.check_type("Unstructured mesh output value", val, bool) self._output = val @@ -1819,7 +1936,7 @@ def volumes(self): return self._volumes @volumes.setter - def volumes(self, volumes): + def volumes(self, volumes: typing.Iterable[Real]): cv.check_type("Unstructured mesh volumes", volumes, Iterable, Real) self._volumes = volumes @@ -1851,7 +1968,7 @@ def n_elements(self): return self._n_elements @n_elements.setter - def n_elements(self, val): + def n_elements(self, val: int): cv.check_type('Number of elements', val, Integral) self._n_elements = val @@ -1883,7 +2000,7 @@ def __repr__(self): self.length_multiplier) return string - def centroid(self, bin): + def centroid(self, bin: int): """Return the vertex averaged centroid of an element Parameters @@ -1927,7 +2044,12 @@ def write_vtk_mesh(self, **kwargs): ) self.write_data_to_vtk(**kwargs) - def write_data_to_vtk(self, filename=None, datasets=None, volume_normalization=True): + def write_data_to_vtk( + self, + filename: Optional[PathLike] = None, + datasets: Optional[dict] = None, + volume_normalization: bool = True + ): """Map data to unstructured VTK mesh elements. Parameters @@ -2019,12 +2141,12 @@ def write_data_to_vtk(self, filename=None, datasets=None, volume_normalization=T writer.Write() @classmethod - def from_hdf5(cls, group): + def from_hdf5(cls, group: h5py.Group): mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) filename = group['filename'][()].decode() library = group['library'][()].decode() - mesh = cls(filename, library, mesh_id=mesh_id) + mesh = cls(filename=filename, library=library, mesh_id=mesh_id) vol_data = group['volumes'][()] mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],)) mesh.n_elements = mesh.volumes.size @@ -2063,7 +2185,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate unstructured mesh object from XML element Parameters diff --git a/openmc/mgxs/__init__.py b/openmc/mgxs/__init__.py index 459b6a874db..4f12a8fb42c 100644 --- a/openmc/mgxs/__init__.py +++ b/openmc/mgxs/__init__.py @@ -13,14 +13,20 @@ - "XMAS-172_" designed for LWR analysis ([SAR1990]_, [SAN2004]_) - "SHEM-361_" designed for LWR analysis to eliminate self-shielding calculations of thermal resonances ([HFA2005]_, [SAN2007]_, [HEB2008]_) -- "SCALE-44" designed for criticality analysis ([ZAL1999]_) +- "SCALE-X" (where X is 44 which is designed for criticality analysis + and 252 is designed for thermal reactors) for the SCALE code suite + ([ZAL1999]_ and [REARDEN2013]_) +- "MPACT-X" (where X is 51 (PWR), 60 (BWR), 69 (Magnox)) from the MPACT_ reactor + physics code ([KIM2019]_ and [KIM2020]_) - "ECCO-1968_" designed for fine group reactor cell calculations for fast, intermediate and thermal reactor applications ([SAR1990]_) - activation_ energy group structures "VITAMIN-J-42", "VITAMIN-J-175", "TRIPOLI-315", "CCFE-709_" and "UKAEA-1102_" .. _CASMO: https://www.studsvik.com/SharepointFiles/CASMO-5%20Development%20and%20Applications.pdf -.. _SCALE: https://www-nds.iaea.org/publications/indc/indc-czr-0001.pdf +.. _SCALE44: https://www-nds.iaea.org/publications/indc/indc-czr-0001.pdf +.. _SCALE252: https://oecd-nea.org/science/wpncs/amct/workingarea/meeting2013/EGAMCT2013_08.pdf +.. _MPACT: https://vera.ornl.gov/mpact/ .. _XMAS-172: https://www-nds.iaea.org/wimsd/energy.htm .. _SHEM-361: https://www.polymtl.ca/merlin/downloads/FP214.pdf .. _activation: https://fispact.ukaea.uk/wiki/Keyword:GETXS @@ -47,7 +53,18 @@ Santamarina-Hfaiedh energy mesh between 22.5 eV and 11.4 keV. International Conference on the Physics of Reactors 2008, PHYSOR 08. 2. 929-938. .. [ZAL1999] K. Záleský and L. Marková (1999), Assessment of Nuclear Data Needs - for Broad-Group SCALE Library Related to VVER Spent Fuel Applications, IAEA. + for Broad-Group SCALE Library Related to VVER Spent Fuel Applications, IAEA. SCALE44_. +.. [REARDEN2013] B. T. Rearden, M. E. Dunn, D. Wiarda, C. Celik, K. Bekar, + M. L. Williams, D. E. Peplow, M. A. Jessee, C. M. Perfetti, + I. C. Gauld, W. A. Wieselquist, J. P. Lefebvre, R. A. Lefebvre, + W. J. Marshall, A. B. Thompson, F. Havluj, S. E. Skutnik, + K. J. Dugan. (2013). Overview of SCALE 6.2. OECD. SCALE252_. +.. [KIM2019] Kim, K.S., Williams, M., Wiarda, D., & Clarno, K. (2019). Development + of the multigroup cross section library for the CASL neutronics simulator MPACT: + Method and procedure. Annals of Nuclear Energy, 133. pp. 46-58. +.. [KIM2020] Kim, K.S., Ade, B., & Luciano, N. (2020). Development + of the MPACT 69-group Library for Magnox Reactor Analysis using VERA. + Proceedings of International Conference on Physics of Reactors PHYSOR2020. """ GROUP_STRUCTURES['CASMO-2'] = np.array([ @@ -80,6 +97,30 @@ 3.25e-1, 3.5e-1, 3.75e-1, 4.e-1, 6.25e-1, 1., 1.77, 3., 4.75, 6., 8.1, 1.e1, 3.e1, 1.e2, 5.5e2, 3.e3, 1.7e4, 2.5e4, 1.e5, 4.e5, 9.e5, 1.4e6, 1.85e6, 2.354e6, 2.479e6, 3.e6, 4.8e6, 6.434e6, 8.1873e6, 2.e7]) +GROUP_STRUCTURES['MPACT-51'] = np.array([ + 0., 1.e-2, 3.e-2, 4.e-2, 6.e-2, 8.e-2, 1.e-1, 1.5e-1, 2.e-1, 2.75e-1, + 3.5e-1, 5.e-1, 6.25e-1, 7.5e-1, 9.25e-1, 9.75e-1, 1.010, 1.080, 1.130, + 1.175, 1.250, 1.450, 1.860, 2.470, 3.730, 4.700, 5.000, 5.400, 6.250, + 7.150, 8.100, 1.19e+1, 1.44e+1, 3.e+1, 4.83e+1, 7.6e+1, 1.43e+2, 3.05e+2, + 9.5e+2, 2.25e+3, 9.5e+3, 2.e+4, 5.e+4, 7.3e+4, 2.e+5, 4.92e+5, 8.2e+5, + 1.356e+6, 2.354e+6, 4.304e+6, 6.434e+6, 2.e+7]) +GROUP_STRUCTURES['MPACT-60'] = np.array([ + 0., 1.e-2, 3.e-2, 4.e-2, 6.e-2, 8.e-2, 1.e-1, 1.5e-1, 2.e-1, 2.75e-1, + 3.5e-1, 5.e-1, 6.25e-1, 7.5e-1, 9.25e-1, 9.75e-1, 1.01, 1.08, 1.13, + 1.175, 1.25, 1.45, 1.86, 2.47, 3.73, 4.7, 5., 5.4, 6.25, 7.15, 8.1, + 1.19e+1, 1.44e+1, 3.e+1, 4.83e+1, 7.6e+1, 1.43e+2, 2.095e+2, 3.05e+2, + 6.7e+2, 9.5e+2, 1.55e+3, 2.25e+3, 3.9e+3, 9.5e+3, 1.3e+4, 2.e+4, 3.e+4, + 5.e+4, 7.3e+4, 1.283e+5, 2.e+5, 3.3e+5, 4.92e+5, 6.7e+5, 8.2e+5, 1.356e+6, + 2.354e+6, 4.304e+6, 6.434e+6, 2.e+7]) +GROUP_STRUCTURES['MPACT-69'] = np.array([ + 0., 1.e-2, 3.e-2, 4.e-2, 6.e-2, 8.e-2, 9.e-2, 1.e-1, 1.25e-1, 1.5e-1, + 1.75e-1, 2.e-1, 2.25e-1, 2.5e-1, 2.75e-1, 3.e-1, 3.25e-1, 3.5e-1, 3.75e-1, + 4.e-1, 4.5e-1, 5.e-1, 5.5e-1, 6.e-1, 6.25e-1, 6.5e-1, 7.5e-1, 8.5e-1, + 9.25e-1, 9.75e-1, 1.01, 1.08, 1.13, 1.175, 1.25, 1.45, 1.86, 2.47, 3., + 3.73, 4.7, 5., 5.4, 6.25, 7.15, 8.1, 1.e+1, 1.19e+1, 1.44e+1, 3.e+1, + 4.83e+1, 7.6e+1, 1.43e+2, 3.05e+2, 5.5e+2, 9.5e+2, 2.25e+3, 3.9e+3, 9.5e+3, + 2.e+4, 5.e+4, 7.3e+4, 2.e+5, 4.92e+5, 8.2e+5, 1.356e+6, 2.354e+6, 4.304e+6, + 6.434e+6, 2.e+7]) GROUP_STRUCTURES['CASMO-70'] = np.array([ 0., 5.e-3, 1.e-2, 1.5e-2, 2.e-2, 2.5e-2, 3.e-2, 3.5e-2, 4.2e-2, 5.e-2, 5.8e-2, 6.7e-2, 8.e-2, 1.e-1, 1.4e-1, 1.8e-1, 2.2e-1, @@ -157,6 +198,34 @@ 1.1052e7, 1.1618e7, 1.2214e7, 1.2523e7, 1.2840e7, 1.3499e7, 1.3840e7, 1.4191e7, 1.4550e7, 1.4918e7, 1.5683e7, 1.6487e7, 1.6905e7, 1.7332e7, 1.9640e7]) +GROUP_STRUCTURES['SCALE-252'] = np.array([ + 0., 1.e-4, 5.e-4, 7.5e-4, 1.e-3, 1.2e-3, 1.5e-3, 2.e-3, 2.5e-3, 3.e-3, + 4.e-3, 5.e-3, 7.5e-3, 1.e-2, 2.53e-2, 3.e-2, 4.e-2, 5.e-2, 6.e-2, 7.e-2, + 8.e-2, 9.e-2, 1.e-1, 1.25e-1, 1.5e-1, 1.75e-1, 2.e-1, 2.25e-1, 2.5e-1, + 2.75e-1, 3.e-1, 3.25e-1, 3.5e-1, 3.75e-1, 4.e-1, 4.5e-1, 5.e-1, 5.5e-1, + 6.e-1, 6.25e-1, 6.5e-1, 7.e-1, 7.5e-1, 8.e-1, 8.5e-1, 9.e-1, 9.25e-1, + 9.5e-1, 9.75e-1, 1., 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09, + 1.1, 1.11, 1.12, 1.13, 1.14, 1.15, 1.175, 1.2, 1.225, 1.25, 1.3, 1.35, 1.4, + 1.45, 1.5, 1.59, 1.68, 1.77, 1.86, 1.94, 2., 2.12, 2.21, 2.3, 2.38, 2.47, + 2.57, 2.67, 2.77, 2.87, 2.97, 3., 3.1, 3.2, 3.5, 3.73, 4.1, 4.7, 5., 5.4, + 6., 6.25, 6.5, 6.75, 6.875, 7., 7.15, 8.1, 9.1, 1.e+1, 1.15e+1, 1.19e+1, + 1.29e+1, 1.44e+1, 1.6e+1, 1.7e+1, 1.85e+1, 1.94e+1, 2.e+1, 2.05e+1, + 2.12e+1, 2.175e+1, 2.25e+1, 2.5e+1, 2.75e+1, 3.e+1, 3.125e+1, 3.175e+1, + 3.325e+1, 3.375e+1, 3.5e+1, 3.55e+1, 3.6e+1, 3.7e+1, 3.713e+1, 3.727e+1, + 3.763e+1, 3.8e+1, 3.91e+1, 3.96e+1, 4.1e+1, 4.24e+1, 4.4e+1, 4.52e+1, + 4.83e+1, 5.06e+1, 5.34e+1, 5.8e+1, 6.1e+1, 6.3e+1, 6.5e+1, 6.75e+1, 7.2e+1, + 7.6e+1, 8.e+1, 8.17e+1, 9.e+1, 9.7e+1, 1.012e+2, 1.05e+2, 1.08e+2, 1.13e+2, + 1.16e+2, 1.175e+2, 1.19e+2, 1.22e+2, 1.43e+2, 1.7e+2, 1.8e+2, 1.877e+2, + 1.885e+2, 1.915e+2, 1.93e+2, 2.02e+2, 2.074e+2, 2.095e+2, 2.2e+2, 2.4e+2, + 2.85e+2, 3.05e+2, 5.5e+2, 6.7e+2, 6.83e+2, 9.5e+2, 1.15e+3, 1.5e+3, + 1.55e+3, 1.8e+3, 2.2e+3, 2.25e+3, 2.5e+3, 3.e+3, 3.74e+3, 3.9e+3, 5.7e+3, + 8.03e+3, 9.5e+3, 1.3e+4, 1.7e+4, 2.e+4, 3.e+4, 4.5e+4, 5.e+4, 5.2e+4, + 6.e+4, 7.3e+4, 7.5e+4, 8.2e+4, 8.5e+4, 1.e+5, 1.283e+5, 1.49e+5, 2.e+5, + 2.7e+5, 3.3e+5, 4.e+5, 4.2e+5, 4.4e+5, 4.7e+5, 4.92e+5, 5.5e+5, 5.73e+5, + 6.e+5, 6.7e+5, 6.79e+5, 7.5e+5, 8.2e+5, 8.611e+5, 8.75e+5, 9.e+5, 9.2e+5, + 1.01e+6, 1.1e+6, 1.2e+6, 1.25e+6, 1.317e+6, 1.356e+6, 1.4e+6, 1.5e+6, + 1.85e+6, 2.354e+6, 2.479e+6, 3.e+6, 4.304e+6, 4.8e+6, 6.434e+6, 8.187e+6, + 1.e+7, 1.284e+7, 1.384e+7, 1.455e+7, 1.568e+7, 1.733e+7, 2.e+7]) GROUP_STRUCTURES['TRIPOLI-315'] = np.array([ 1.0e-5, 1.1e-4, 3.000e-3, 5.500e-3, 1.000e-2, 1.500e-2, 2.000e-2, 3.000e-2, 3.200e-2, 3.238e-2, 4.300e-2, 5.900e-2, 7.700e-2, 9.500e-2, 1.000e-1, diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 988b7cc7566..2fa2e0f0567 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -1,4 +1,3 @@ -from collections import OrderedDict from collections.abc import Iterable import copy from numbers import Integral @@ -83,7 +82,7 @@ class Library: tally_trigger : openmc.Trigger An (optional) tally precision trigger given to each tally used to compute the cross section - all_mgxs : collections.OrderedDict + all_mgxs : dict MGXS objects keyed by domain ID and cross section type sp_filename : str The filename of the statepoint with tally data used to the @@ -119,7 +118,7 @@ def __init__(self, geometry, by_nuclide=False, self._legendre_order = 0 self._histogram_bins = 16 self._tally_trigger = None - self._all_mgxs = OrderedDict() + self._all_mgxs = {} self._sp_filename = None self._keff = None self._sparse = False @@ -159,9 +158,9 @@ def __deepcopy__(self, memo): clone._keff = self._keff clone._sparse = self.sparse - clone._all_mgxs = OrderedDict() + clone._all_mgxs = {} for domain in self.domains: - clone.all_mgxs[domain.id] = OrderedDict() + clone.all_mgxs[domain.id] = {} for mgxs_type in self.mgxs_types: mgxs = copy.deepcopy(self.all_mgxs[domain.id][mgxs_type]) clone.all_mgxs[domain.id][mgxs_type] = mgxs @@ -191,7 +190,7 @@ def name(self): def name(self, name): cv.check_type('name', name, str) self._name = name - + @property def mgxs_types(self): return self._mgxs_types @@ -306,7 +305,7 @@ def nuclides(self, nuclides): @property def energy_groups(self): return self._energy_groups - + @energy_groups.setter def energy_groups(self, energy_groups): cv.check_type('energy groups', energy_groups, openmc.mgxs.EnergyGroups) @@ -324,7 +323,7 @@ def num_delayed_groups(self, num_delayed_groups): cv.check_greater_than('num delayed groups', num_delayed_groups, 0, equality=True) self._num_delayed_groups = num_delayed_groups - + @property def num_polar(self): return self._num_polar @@ -382,7 +381,7 @@ def scatter_format(self, scatter_format): self.correction = None self._scatter_format = scatter_format - + @property def legendre_order(self): return self._legendre_order @@ -503,7 +502,7 @@ def build_library(self): # Initialize MGXS for each domain and mgxs type and store in dictionary for domain in self.domains: - self.all_mgxs[domain.id] = OrderedDict() + self.all_mgxs[domain.id] = {} for mgxs_type in self.mgxs_types: if mgxs_type in openmc.mgxs.MDGXS_TYPES: mgxs = openmc.mgxs.MDGXS.get_mgxs( diff --git a/openmc/mgxs/mdgxs.py b/openmc/mgxs/mdgxs.py index 7643757d469..3a6732748bb 100644 --- a/openmc/mgxs/mdgxs.py +++ b/openmc/mgxs/mdgxs.py @@ -1,4 +1,3 @@ -from collections import OrderedDict import copy import itertools from numbers import Integral @@ -90,7 +89,7 @@ class MDGXS(MGXS): the multi-group cross section estimator : {'tracklength', 'analog'} The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section rxn_rate_tally : openmc.Tally Derived tally for the reaction rate tally used in the numerator to @@ -161,7 +160,7 @@ def __deepcopy__(self, memo): clone._sparse = self.sparse clone._derived = self.derived - clone._tallies = OrderedDict() + clone._tallies = {} for tally_type, tally in self.tallies.items(): clone.tallies[tally_type] = copy.deepcopy(tally, memo) @@ -972,7 +971,7 @@ class ChiDelayed(MDGXS): the multi-group cross section estimator : 'analog' The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`ChiDelayed.tally_keys` property and values are instances of :class:`openmc.Tally`. @@ -1491,7 +1490,7 @@ class DelayedNuFissionXS(MDGXS): the multi-group cross section estimator : {'tracklength', 'analog'} The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`DelayedNuFissionXS.tally_keys` property and values are instances of :class:`openmc.Tally`. @@ -1630,7 +1629,7 @@ class Beta(MDGXS): the multi-group cross section estimator : {'tracklength', 'analog'} The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`Beta.tally_keys` property and values are instances of :class:`openmc.Tally`. @@ -1823,7 +1822,7 @@ class DecayRate(MDGXS): the multi-group cross section estimator : {'tracklength', 'analog'} The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`DecayRate.tally_keys` property and values are instances of :class:`openmc.Tally`. @@ -2137,7 +2136,7 @@ class MatrixMDGXS(MDGXS): the multi-group cross section estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section rxn_rate_tally : openmc.Tally Derived tally for the reaction rate tally used in the numerator to @@ -2735,7 +2734,7 @@ class DelayedNuFissionMatrixXS(MatrixMDGXS): the multi-group cross section estimator : 'analog' The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`DelayedNuFissionXS.tally_keys` property and values are instances of :class:`openmc.Tally`. diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 264cbc8924c..b6840d079ce 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -1,4 +1,3 @@ -from collections import OrderedDict import copy from numbers import Integral import os @@ -216,7 +215,7 @@ class MGXS: the multi-group cross section estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section rxn_rate_tally : openmc.Tally Derived tally for the reaction rate tally used in the numerator to @@ -319,7 +318,7 @@ def __deepcopy__(self, memo): clone._derived = self.derived clone._mgxs_type = self._mgxs_type - clone._tallies = OrderedDict() + clone._tallies = {} for tally_type, tally in self.tallies.items(): clone.tallies[tally_type] = copy.deepcopy(tally, memo) @@ -576,7 +575,7 @@ def tallies(self): if self._tallies is None: # Initialize a collection of Tallies - self._tallies = OrderedDict() + self._tallies ={} # Create a domain Filter object filter_type = _DOMAIN_TO_FILTER[self.domain_type] @@ -1690,7 +1689,7 @@ def merge(self, other): merged_mgxs.nuclides = self.nuclides + other.nuclides # Null base tallies but merge reaction rate and cross section tallies - merged_mgxs._tallies = OrderedDict() + merged_mgxs._tallies ={} merged_mgxs._rxn_rate_tally = self.rxn_rate_tally.merge(other.rxn_rate_tally) merged_mgxs._xs_tally = self.xs_tally.merge(other.xs_tally) @@ -2690,7 +2689,7 @@ class TransportXS(MGXS): the multi-group cross section estimator : 'analog' The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`TransportXS.tally_keys` property and values are instances of :class:`openmc.Tally`. @@ -2931,7 +2930,7 @@ class DiffusionCoefficient(TransportXS): the multi-group cross section estimator : 'analog' The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`TransportXS.tally_keys` property and values are instances of :class:`openmc.Tally`. @@ -3078,7 +3077,7 @@ def get_condensed_xs(self, coarse_groups): diff_coef = transport**(-1) / 3.0 diff_coef *= self.tallies['flux (tracklength)'] flux_tally = condensed_xs.tallies['flux (tracklength)'] - condensed_xs._tallies = OrderedDict() + condensed_xs._tallies = {} condensed_xs._tallies[self._rxn_type] = diff_coef condensed_xs._tallies['flux (tracklength)'] = flux_tally condensed_xs._rxn_rate_tally = diff_coef @@ -3394,7 +3393,7 @@ class FissionXS(MGXS): the multi-group cross section estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`FissionXS.tally_keys` property and values are instances of :class:`openmc.Tally`. @@ -3613,7 +3612,7 @@ class ScatterXS(MGXS): the multi-group cross section estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`ScatterXS.tally_keys` property and values are instances of :class:`openmc.Tally`. @@ -3923,7 +3922,7 @@ class ScatterMatrixXS(MatrixMGXS): the multi-group cross section estimator : 'analog' The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`ScatterMatrixXS.tally_keys` property and values are instances of :class:`openmc.Tally`. @@ -5193,7 +5192,7 @@ class NuFissionMatrixXS(MatrixMGXS): the multi-group cross section estimator : 'analog' The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`NuFissionMatrixXS.tally_keys` property and values are instances of :class:`openmc.Tally`. @@ -5354,7 +5353,7 @@ class Chi(MGXS): the multi-group cross section estimator : 'analog' The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`Chi.tally_keys` property and values are instances of :class:`openmc.Tally`. @@ -5922,7 +5921,7 @@ class MeshSurfaceMGXS(MGXS): the multi-group cross section estimator : {'analog'} The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section rxn_rate_tally : openmc.Tally Derived tally for the reaction rate tally used in the numerator to @@ -6309,7 +6308,7 @@ class Current(MeshSurfaceMGXS): the multi-group cross section estimator : {'analog'} The tally estimator used to compute the multi-group cross section - tallies : collections.OrderedDict + tallies : dict OpenMC tallies needed to compute the multi-group cross section. The keys are strings listed in the :attr:`TotalXS.tally_keys` property and values are instances of :class:`openmc.Tally`. diff --git a/openmc/plots.py b/openmc/plots.py index 5bc8a2e6e8a..6afe3ea33b8 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -11,7 +11,7 @@ import openmc.checkvalue as cv from openmc.checkvalue import PathLike -from ._xml import clean_indentation, get_elem_tuple, reorder_attributes +from ._xml import clean_indentation, get_elem_tuple, reorder_attributes, get_text from .mixin import IDManagerMixin _BASES = ['xy', 'xz', 'yz'] @@ -482,6 +482,9 @@ def to_xml_element(self): """ element = ET.Element("plot") + element.set("id", str(self._id)) + if len(self._name) > 0: + element.set("name", str(self.name)) if self._filename is not None: element.set("filename", self._filename) element.set("color_by", self._color_by) @@ -782,7 +785,6 @@ def to_xml_element(self): """ element = super().to_xml_element() - element.set("id", str(self._id)) element.set("type", self._type) if self._type == 'slice': @@ -843,12 +845,14 @@ def from_xml_element(cls, elem): """ plot_id = int(elem.get("id")) - plot = cls(plot_id) + name = get_text(elem, 'name', '') + plot = cls(plot_id, name) if "filename" in elem.keys(): plot.filename = elem.get("filename") plot.color_by = elem.get("color_by") plot.type = elem.get("type") - plot.basis = elem.get("basis") + if plot.type == 'slice': + plot.basis = elem.get("basis") plot.origin = get_elem_tuple(elem, "origin", float) plot.width = get_elem_tuple(elem, "width", float) @@ -1176,7 +1180,6 @@ def to_xml_element(self): """ element = super().to_xml_element() - element.set("id", str(self._id)) element.set("type", "projection") subelement = ET.SubElement(element, "camera_position") @@ -1276,7 +1279,7 @@ def from_xml_element(cls, elem): tmp = elem.find("orthographic_width") if tmp is not None: - self.orthographic_width = float(tmp) + plot.orthographic_width = float(tmp) plot.pixels = get_elem_tuple(elem, "pixels") plot.camera_position = get_elem_tuple(elem, "camera_position", float) diff --git a/openmc/region.py b/openmc/region.py index af9df1ba1e0..a2d3b1d0cbb 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -1,5 +1,4 @@ from abc import ABC, abstractmethod -from collections import OrderedDict from collections.abc import MutableSequence from copy import deepcopy @@ -47,17 +46,17 @@ def get_surfaces(self, surfaces=None): Parameters ---------- - surfaces: collections.OrderedDict, optional + surfaces : dict, optional Dictionary mapping surface IDs to :class:`openmc.Surface` instances Returns ------- - surfaces: collections.OrderedDict + surfaces : dict Dictionary mapping surface IDs to :class:`openmc.Surface` instances """ if surfaces is None: - surfaces = OrderedDict() + surfaces = {} for region in self: surfaces = region.get_surfaces(surfaces) return surfaces @@ -590,17 +589,17 @@ def get_surfaces(self, surfaces=None): Parameters ---------- - surfaces: collections.OrderedDict, optional + surfaces : dict, optional Dictionary mapping surface IDs to :class:`openmc.Surface` instances Returns ------- - surfaces: collections.OrderedDict + surfaces : dict Dictionary mapping surface IDs to :class:`openmc.Surface` instances """ if surfaces is None: - surfaces = OrderedDict() + surfaces = {} for region in self.node: surfaces = region.get_surfaces(surfaces) return surfaces diff --git a/openmc/settings.py b/openmc/settings.py index 7937cd2e7a3..ef20d54dbbc 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -51,14 +51,16 @@ class Settings: create_fission_neutrons : bool Indicate whether fission neutrons should be created or not. cutoff : dict - Dictionary defining weight cutoff and energy cutoff. The dictionary may - have six keys, 'weight', 'weight_avg', 'energy_neutron', 'energy_photon', - 'energy_electron', and 'energy_positron'. Value for 'weight' + Dictionary defining weight cutoff, energy cutoff and time cutoff. The + dictionary may have ten keys, 'weight', 'weight_avg', 'energy_neutron', + 'energy_photon', 'energy_electron', 'energy_positron', 'time_neutron', + 'time_photon', 'time_electron', and 'time_positron'. Value for 'weight' should be a float indicating weight cutoff below which particle undergo Russian roulette. Value for 'weight_avg' should be a float indicating - weight assigned to particles that are not killed after Russian - roulette. Value of energy should be a float indicating energy in eV - below which particle type will be killed. + weight assigned to particles that are not killed after Russian roulette. + Value of energy should be a float indicating energy in eV below which + particle type will be killed. Value of time should be a float in + seconds. Particles will be killed exactly at the specified time. delayed_photon_scaling : bool Indicate whether to scale the fission photon yield by (EGP + EGD)/EGP where EGP is the energy release of prompt photons and EGD is the energy @@ -133,6 +135,8 @@ class Settings: Number of particles per generation photon_transport : bool Whether to use photon transport. + plot_seed : int + Initial seed for randomly generated plot colors. ptables : bool Determine whether probability tables are used. resonance_scattering : dict @@ -270,6 +274,7 @@ def __init__(self, **kwargs): self._confidence_intervals = None self._electron_treatment = None self._photon_transport = None + self._plot_seed = None self._ptables = None self._seed = None self._survival_biasing = None @@ -504,6 +509,16 @@ def photon_transport(self, photon_transport: bool): cv.check_type('photon transport', photon_transport, bool) self._photon_transport = photon_transport + @property + def plot_seed(self): + return self._plot_seed + + @plot_seed.setter + def plot_seed(self, seed): + cv.check_type('random plot color seed', seed, Integral) + cv.check_greater_than('random plot color seed', seed, 0) + self._plot_seed = seed + @property def seed(self) -> int: return self._seed @@ -1118,6 +1133,11 @@ def _create_photon_transport_subelement(self, root): element = ET.SubElement(root, "photon_transport") element.text = str(self._photon_transport).lower() + def _create_plot_seed_subelement(self, root): + if self._plot_seed is not None: + element = ET.SubElement(root, "plot_seed") + element.text = str(self._plot_seed) + def _create_ptables_subelement(self, root): if self._ptables is not None: element = ET.SubElement(root, "ptables") @@ -1491,6 +1511,11 @@ def _photon_transport_from_xml_element(self, root): if text is not None: self.photon_transport = text in ('true', '1') + def _plot_seed_from_xml_element(self, root): + text = get_text(root, 'plot_seed') + if text is not None: + self.plot_seed = int(text) + def _ptables_from_xml_element(self, root): text = get_text(root, 'ptables') if text is not None: @@ -1511,7 +1536,8 @@ def _cutoff_from_xml_element(self, root): if elem is not None: self.cutoff = {} for key in ('energy_neutron', 'energy_photon', 'energy_electron', - 'energy_positron', 'weight', 'weight_avg'): + 'energy_positron', 'weight', 'weight_avg', 'time_neutron', + 'time_photon', 'time_electron', 'time_positron'): value = get_text(elem, key) if value is not None: self.cutoff[key] = float(value) @@ -1706,6 +1732,7 @@ def to_xml_element(self, mesh_memo=None): self._create_energy_mode_subelement(element) self._create_max_order_subelement(element) self._create_photon_transport_subelement(element) + self._create_plot_seed_subelement(element) self._create_ptables_subelement(element) self._create_seed_subelement(element) self._create_survival_biasing_subelement(element) @@ -1802,6 +1829,7 @@ def from_xml_element(cls, elem, meshes=None): settings._energy_mode_from_xml_element(elem) settings._max_order_from_xml_element(elem) settings._photon_transport_from_xml_element(elem) + settings._plot_seed_from_xml_element(elem) settings._ptables_from_xml_element(elem) settings._seed_from_xml_element(elem) settings._survival_biasing_from_xml_element(elem) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 058bee7da69..09a2f8752ff 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -1,15 +1,17 @@ +from __future__ import annotations +import typing from abc import ABC, abstractmethod from collections.abc import Iterable -from math import pi, cos +from math import cos, pi from numbers import Real -import lxml.etree as ET +import lxml.etree as ET import numpy as np import openmc.checkvalue as cv from .._xml import get_text -from .univariate import Univariate, Uniform, PowerLaw from ..mesh import MeshBase +from .univariate import PowerLaw, Uniform, Univariate class UnitSphere(ABC): @@ -177,7 +179,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate isotropic distribution from an XML element Parameters @@ -209,7 +211,7 @@ class Monodirectional(UnitSphere): """ - def __init__(self, reference_uvw=[1., 0., 0.]): + def __init__(self, reference_uvw: typing.Sequence[float] = [1., 0., 0.]): super().__init__(reference_uvw) def to_xml_element(self): @@ -228,7 +230,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate monodirectional distribution from an XML element Parameters @@ -304,7 +306,12 @@ class CartesianIndependent(Spatial): """ - def __init__(self, x, y, z): + def __init__( + self, + x: openmc.stats.Univariate, + y: openmc.stats.Univariate, + z: openmc.stats.Univariate + ): self.x = x self.y = y self.z = z @@ -353,7 +360,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate spatial distribution from an XML element Parameters @@ -477,7 +484,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate spatial distribution from an XML element Parameters @@ -599,7 +606,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate spatial distribution from an XML element Parameters @@ -772,7 +779,12 @@ class Box(Spatial): """ - def __init__(self, lower_left, upper_right, only_fissionable=False): + def __init__( + self, + lower_left: typing.Sequence[float], + upper_right: typing.Sequence[float], + only_fissionable: bool = False + ): self.lower_left = lower_left self.upper_right = upper_right self.only_fissionable = only_fissionable @@ -826,7 +838,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate box distribution from an XML element Parameters @@ -865,7 +877,7 @@ class Point(Spatial): """ - def __init__(self, xyz=(0., 0., 0.)): + def __init__(self, xyz: typing.Sequence[float] = (0., 0., 0.)): self.xyz = xyz @property @@ -894,7 +906,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate point distribution from an XML element Parameters @@ -912,8 +924,13 @@ def from_xml_element(cls, elem): return cls(xyz) -def spherical_uniform(r_outer, r_inner=0.0, thetas=(0., pi), phis=(0., 2*pi), - origin=(0., 0., 0.)): +def spherical_uniform( + r_outer: float, + r_inner: float = 0.0, + thetas: typing.Sequence[float] = (0., pi), + phis: typing.Sequence[float] = (0., 2*pi), + origin: typing.Sequence[float] = (0., 0., 0.) + ): """Return a uniform spatial distribution over a spherical shell. This function provides a uniform spatial distribution over a spherical @@ -926,15 +943,15 @@ def spherical_uniform(r_outer, r_inner=0.0, thetas=(0., pi), phis=(0., 2*pi), ---------- r_outer : float Outer radius of the spherical shell in [cm] - r_inner : float, optional + r_inner : float Inner radius of the spherical shell in [cm] - thetas : iterable of float, optional + thetas : iterable of float Starting and ending theta coordinates (angle relative to the z-axis) in radius in a reference frame centered at `origin` - phis : iterable of float, optional + phis : iterable of float Starting and ending phi coordinates (azimuthal angle) in radians in a reference frame centered at `origin` - origin: iterable of float, optional + origin: iterable of float Coordinates (x0, y0, z0) of the center of the spherical reference frame for the distribution. diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index b5720e23cde..ff1af2d7ba6 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -1,19 +1,19 @@ +import math +import typing from abc import ABC, abstractmethod from collections import defaultdict from collections.abc import Iterable from copy import deepcopy -import math from numbers import Real from warnings import warn -import lxml.etree as ET +import lxml.etree as ET import numpy as np import openmc.checkvalue as cv from .._xml import get_text from ..mixin import EqualityMixin - _INTERPOLATION_SCHEMES = [ 'histogram', 'linear-linear', @@ -66,7 +66,7 @@ def from_xml_element(cls, elem): return Mixture.from_xml_element(elem) @abstractmethod - def sample(n_samples=1, seed=None): + def sample(n_samples: int = 1, seed: typing.Optional[int] = None): """Sample the univariate distribution Parameters @@ -186,7 +186,7 @@ def to_xml_element(self, element_name): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate discrete distribution from an XML element Parameters @@ -206,7 +206,11 @@ def from_xml_element(cls, elem): return cls(x, p) @classmethod - def merge(cls, dists, probs): + def merge( + cls, + dists: typing.Sequence['openmc.stats.Discrete'], + probs: typing.Sequence[int] + ): """Merge multiple discrete distributions into a single distribution .. versionadded:: 0.13.1 @@ -271,7 +275,7 @@ class Uniform(Univariate): """ - def __init__(self, a=0.0, b=1.0): + def __init__(self, a: float = 0.0, b: float = 1.0): self.a = a self.b = b @@ -306,7 +310,7 @@ def sample(self, n_samples=1, seed=None): np.random.seed(seed) return np.random.uniform(self.a, self.b, n_samples) - def to_xml_element(self, element_name): + def to_xml_element(self, element_name: str): """Return XML representation of the uniform distribution Parameters @@ -326,7 +330,7 @@ def to_xml_element(self, element_name): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate uniform distribution from an XML element Parameters @@ -372,7 +376,7 @@ class PowerLaw(Univariate): """ - def __init__(self, a=0.0, b=1.0, n=0): + def __init__(self, a: float = 0.0, b: float = 1.0, n: float = 0.): self.a = a self.b = b self.n = n @@ -415,7 +419,7 @@ def sample(self, n_samples=1, seed=None): span = self.b**pwr - offset return np.power(offset + xi * span, 1/pwr) - def to_xml_element(self, element_name): + def to_xml_element(self, element_name: str): """Return XML representation of the power law distribution Parameters @@ -435,7 +439,7 @@ def to_xml_element(self, element_name): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate power law distribution from an XML element Parameters @@ -493,12 +497,12 @@ def sample(self, n_samples=1, seed=None): return self.sample_maxwell(self.theta, n_samples) @staticmethod - def sample_maxwell(t, n_samples): + def sample_maxwell(t, n_samples: int): r1, r2, r3 = np.random.rand(3, n_samples) c = np.cos(0.5 * np.pi * r3) return -t * (np.log(r1) + np.log(r2) * c * c) - def to_xml_element(self, element_name): + def to_xml_element(self, element_name: str): """Return XML representation of the Maxwellian distribution Parameters @@ -518,7 +522,7 @@ def to_xml_element(self, element_name): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate Maxwellian distribution from an XML element Parameters @@ -593,7 +597,7 @@ def sample(self, n_samples=1, seed=None): aab = self.a * self.a * self.b return w + 0.25*aab + u*np.sqrt(aab*w) - def to_xml_element(self, element_name): + def to_xml_element(self, element_name: str): """Return XML representation of the Watt distribution Parameters @@ -613,7 +617,7 @@ def to_xml_element(self, element_name): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate Watt distribution from an XML element Parameters @@ -683,7 +687,7 @@ def sample(self, n_samples=1, seed=None): np.random.seed(seed) return np.random.normal(self.mean_value, self.std_dev, n_samples) - def to_xml_element(self, element_name): + def to_xml_element(self, element_name: str): """Return XML representation of the Normal distribution Parameters @@ -703,7 +707,7 @@ def to_xml_element(self, element_name): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate Normal distribution from an XML element Parameters @@ -721,7 +725,7 @@ def from_xml_element(cls, elem): return cls(*map(float, params)) -def muir(e0, m_rat, kt): +def muir(e0: float, m_rat: float, kt: float): """Generate a Muir energy spectrum The Muir energy spectrum is a normal distribution, but for convenience @@ -793,8 +797,13 @@ class Tabular(Univariate): """ - def __init__(self, x, p, interpolation='linear-linear', - ignore_negative=False): + def __init__( + self, + x: typing.Sequence[float], + p: typing.Sequence[float], + interpolation: str = 'linear-linear', + ignore_negative: bool = False + ): self._ignore_negative = ignore_negative self.x = x self.p = p @@ -886,7 +895,7 @@ def normalize(self): """Normalize the probabilities stored on the distribution""" self.p /= self.cdf().max() - def sample(self, n_samples=1, seed=None): + def sample(self, n_samples: int = 1, seed: typing.Optional[int] = None): np.random.seed(seed) xi = np.random.rand(n_samples) @@ -947,7 +956,7 @@ def sample(self, n_samples=1, seed=None): assert all(samples_out < self.x[-1]) return samples_out - def to_xml_element(self, element_name): + def to_xml_element(self, element_name: str): """Return XML representation of the tabular distribution Parameters @@ -971,7 +980,7 @@ def to_xml_element(self, element_name): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate tabular distribution from an XML element Parameters @@ -1028,7 +1037,7 @@ class Legendre(Univariate): """ - def __init__(self, coefficients): + def __init__(self, coefficients: typing.Sequence[float]): self.coefficients = coefficients self._legendre_poly = None @@ -1082,7 +1091,11 @@ class Mixture(Univariate): """ - def __init__(self, probability, distribution): + def __init__( + self, + probability: typing.Sequence[float], + distribution: typing.Sequence['openmc.Univariate'] + ): self.probability = probability self.distribution = distribution @@ -1140,7 +1153,7 @@ def normalize(self): norm = sum(self.probability) self.probability = [val / norm for val in self.probability] - def to_xml_element(self, element_name): + def to_xml_element(self, element_name: str): """Return XML representation of the mixture distribution .. versionadded:: 0.13.0 @@ -1167,7 +1180,7 @@ def to_xml_element(self, element_name): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate mixture distribution from an XML element .. versionadded:: 0.13.0 @@ -1207,7 +1220,10 @@ def integral(self): ]) -def combine_distributions(dists, probs): +def combine_distributions( + dists: typing.Sequence['openmc.Univariate'], + probs: typing.Sequence[float] +): """Combine distributions with specified probabilities This function can be used to combine multiple instances of diff --git a/openmc/surface.py b/openmc/surface.py index 1f35e7bff49..6f19de5caae 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1,5 +1,4 @@ from abc import ABC, abstractmethod -from collections import OrderedDict from collections.abc import Iterable from copy import deepcopy import math @@ -2549,17 +2548,17 @@ def get_surfaces(self, surfaces=None): Parameters ---------- - surfaces: collections.OrderedDict, optional + surfaces : dict, optional Dictionary mapping surface IDs to :class:`openmc.Surface` instances Returns ------- - surfaces: collections.OrderedDict + surfaces : dict Dictionary mapping surface IDs to :class:`openmc.Surface` instances """ if surfaces is None: - surfaces = OrderedDict() + surfaces = {} surfaces[self.surface.id] = self.surface return surfaces diff --git a/openmc/trigger.py b/openmc/trigger.py index 86d17dfefd0..8a165526cf6 100644 --- a/openmc/trigger.py +++ b/openmc/trigger.py @@ -29,7 +29,7 @@ class Trigger(EqualityMixin): """ - def __init__(self, trigger_type, threshold): + def __init__(self, trigger_type: str, threshold: float): self.trigger_type = trigger_type self.threshold = threshold self._scores = [] @@ -92,7 +92,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element): """Generate trigger object from an XML element Parameters diff --git a/openmc/universe.py b/openmc/universe.py index c490ac2989c..299709acc19 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -1,7 +1,6 @@ import math import typing from abc import ABC, abstractmethod -from collections import OrderedDict from collections.abc import Iterable from copy import deepcopy from numbers import Integral, Real @@ -46,7 +45,7 @@ def __init__(self, universe_id=None, name=''): # Keys - Cell IDs # Values - Cells - self._cells = OrderedDict() + self._cells = {} def __repr__(self): string = 'Universe\n' @@ -100,13 +99,13 @@ def get_all_universes(self): Returns ------- - universes : collections.OrderedDict + universes : dict Dictionary whose keys are universe IDs and values are :class:`Universe` instances """ # Append all Universes within each Cell to the dictionary - universes = OrderedDict() + universes = {} for cell in self.get_all_cells().values(): universes.update(cell.get_all_universes()) @@ -169,7 +168,7 @@ def clone(self, clone_materials=True, clone_regions=True, memo=None): clone = self._partial_deepcopy() # Clone all cells for the universe clone - clone._cells = OrderedDict() + clone._cells = {} for cell in self._cells.values(): clone.add_cell(cell.clone(clone_materials, clone_regions, memo)) @@ -199,7 +198,7 @@ class Universe(UniverseBase): Unique identifier of the universe name : str Name of the universe - cells : collections.OrderedDict + cells : dict Dictionary whose keys are cell IDs and values are :class:`Cell` instances volume : float @@ -370,8 +369,8 @@ def plot(self, origin=None, width=None, pixels=40000, Returns ------- - matplotlib.image.AxesImage - Resulting image + matplotlib.axes.Axes + Axes containing resulting image """ import matplotlib.image as mpimg @@ -423,7 +422,7 @@ def plot(self, origin=None, width=None, pixels=40000, model = openmc.Model() model.geometry = openmc.Geometry(self) if seed is not None: - model.settings.seed = seed + model.settings.plot_seed = seed # Determine whether any materials contains macroscopic data and if # so, set energy mode accordingly @@ -483,7 +482,7 @@ def plot(self, origin=None, width=None, pixels=40000, # add legend showing which colors represent which material # or cell if that was requested if legend: - if plot.colors is None: + if plot.colors == {}: raise ValueError("Must pass 'colors' dictionary if you " "are adding a legend via legend=True.") @@ -522,7 +521,8 @@ def plot(self, origin=None, width=None, pixels=40000, 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) + axes.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs) + return axes def add_cell(self, cell): """Add a cell to the universe. @@ -610,12 +610,12 @@ def get_nuclide_densities(self): Returns ------- - nuclides : collections.OrderedDict + nuclides : dict Dictionary whose keys are nuclide names and values are 2-tuples of (nuclide, density) """ - nuclides = OrderedDict() + nuclides = {} if self._atoms: volume = self.volume @@ -637,13 +637,13 @@ def get_all_cells(self, memo=None): Returns ------- - cells : collections.OrderedDict + cells : dict Dictionary whose keys are cell IDs and values are :class:`Cell` instances """ - cells = OrderedDict() + cells = {} if memo and self in memo: return cells @@ -665,13 +665,13 @@ def get_all_materials(self, memo=None): Returns ------- - materials : collections.OrderedDict + materials : dict Dictionary whose keys are material IDs and values are :class:`Material` instances """ - materials = OrderedDict() + materials = {} # Append all Cells in each Cell in the Universe to the dictionary cells = self.get_all_cells(memo) @@ -882,10 +882,10 @@ def material_names(self): return sorted(set(material_tags_ascii)) def get_all_cells(self, memo=None): - return OrderedDict() + return {} def get_all_materials(self, memo=None): - return OrderedDict() + return {} def _n_geom_elements(self, geom_type): """ diff --git a/openmc/volume.py b/openmc/volume.py index 72fe89fbc4a..90882267d7b 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -1,4 +1,3 @@ -from collections import OrderedDict from collections.abc import Iterable, Mapping from numbers import Real, Integral import lxml.etree as ET @@ -122,6 +121,10 @@ def __init__(self, domains, samples, lower_left=None, upper_right=None): raise ValueError('Could not automatically determine bounding box ' 'for stochastic volume calculation.') + if np.isinf(self.lower_left).any() or np.isinf(self.upper_right).any(): + raise ValueError('Lower-left and upper-right bounding box ' + 'coordinates must be finite.') + @property def ids(self): return self._ids @@ -281,7 +284,7 @@ def from_hdf5(cls, filename): volumes[domain_id] = volume nucnames = group['nuclides'][()] atoms_ = group['atoms'][()] - atom_dict = OrderedDict() + atom_dict = {} for name_i, atoms_i in zip(nucnames, atoms_): atom_dict[name_i.decode()] = ufloat(*atoms_i) atoms[domain_id] = atom_dict diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 3a0a8a354fd..6ea71fa1cb9 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -166,6 +166,7 @@ def __repr__(self) -> str: string += '{: <16}=\t{}\n'.format('\tMesh', self.mesh) string += '{: <16}=\t{}\n'.format('\tParticle Type', self._particle_type) string += '{: <16}=\t{}\n'.format('\tEnergy Bounds', self._energy_bounds) + string += '{: <16}=\t{}\n'.format('\tMax lower bound ratio', self.max_lower_bound_ratio) string += '{: <16}=\t{}\n'.format('\tLower WW Bounds', self._lower_ww_bounds) string += '{: <16}=\t{}\n'.format('\tUpper WW Bounds', self._upper_ww_bounds) string += '{: <16}=\t{}\n'.format('\tSurvival Ratio', self._survival_ratio) @@ -588,13 +589,19 @@ def wwinp_to_wws(path: PathLike) -> List[WeightWindows]: mesh = RectilinearMesh() mesh.x_grid, mesh.y_grid, mesh.z_grid = grids elif nwg == 2: - mesh = CylindricalMesh() - mesh.r_grid, mesh.z_grid, mesh.phi_grid = grids - mesh.origin = xyz0 + mesh = CylindricalMesh( + r_grid=grids[0], + z_grid=grids[1], + phi_grid=grids[2], + origin = xyz0, + ) elif nwg == 3: - mesh = SphericalMesh() - mesh.r_grid, mesh.theta_grid, mesh.phi_grid = grids - mesh.origin = xyz0 + mesh = SphericalMesh( + r_grid=grids[0], + theta_grid=grids[1], + phi_grid=grids[2], + origin = xyz0 + ) # extract weight window values from array wws = [] diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 4d63f089266..35522d3cafc 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -182,10 +182,11 @@ void DAGUniverse::init_geometry() model::cell_map[c->id_] = model::cells.size(); } else { warning(fmt::format("DAGMC Cell IDs: {}", dagmc_ids_for_dim(3))); - fatal_error(fmt::format("DAGMC Universe {} contains a cell with ID {}, which " - "already exists elsewhere in the geometry. Setting auto_geom_ids " - "to True when initiating the DAGMC Universe may " - "resolve this issue", + fatal_error(fmt::format( + "DAGMC Universe {} contains a cell with ID {}, which " + "already exists elsewhere in the geometry. Setting auto_geom_ids " + "to True when initiating the DAGMC Universe may " + "resolve this issue", this->id_, c->id_)); } @@ -395,7 +396,7 @@ bool DAGUniverse::find_cell(Particle& p) const // cells, place it in the implicit complement bool found = Universe::find_cell(p); if (!found && model::universe_map[this->id_] != model::root_universe) { - p.coord(p.n_coord() - 1).cell = implicit_complement_idx(); + p.lowest_coord().cell = implicit_complement_idx(); found = true; } return found; @@ -581,7 +582,7 @@ std::pair DAGCell::distance( p->history().reset(); } - const auto& univ = model::universes[p->coord(p->n_coord() - 1).universe]; + const auto& univ = model::universes[p->lowest_coord().universe]; DAGUniverse* dag_univ = static_cast(univ.get()); if (!dag_univ) @@ -609,7 +610,8 @@ std::pair DAGCell::distance( // 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}; + 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 @@ -739,7 +741,8 @@ void check_dagmc_root_univ() } } -int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ) { +int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ) +{ auto surfp = dynamic_cast(model::surfaces[surf - 1].get()); auto cellp = dynamic_cast(model::cells[curr_cell].get()); auto univp = static_cast(model::universes[univ].get()); @@ -750,11 +753,12 @@ int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ) { cellp->dagmc_ptr()->entity_by_index(3, cellp->dag_index()); moab::EntityHandle new_vol; - moab::ErrorCode rval = cellp->dagmc_ptr()->next_vol(surf_handle, curr_vol, new_vol); - if (rval != moab::MB_SUCCESS) return -1; + moab::ErrorCode rval = + cellp->dagmc_ptr()->next_vol(surf_handle, curr_vol, new_vol); + if (rval != moab::MB_SUCCESS) + return -1; - return cellp->dagmc_ptr()->index_by_handle(new_vol) + - univp->cell_idx_offset_; + return cellp->dagmc_ptr()->index_by_handle(new_vol) + univp->cell_idx_offset_; } } // namespace openmc diff --git a/src/event.cpp b/src/event.cpp index 1db5c99d78e..e9490a77f13 100644 --- a/src/event.cpp +++ b/src/event.cpp @@ -112,6 +112,8 @@ void process_advance_particle_events() int64_t buffer_idx = simulation::advance_particle_queue[i].idx; Particle& p = simulation::particles[buffer_idx]; p.event_advance(); + if (!p.alive()) + continue; if (p.collision_distance() > p.boundary().distance) { simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx}); } else { diff --git a/src/finalize.cpp b/src/finalize.cpp index 59294b28c9d..f3fe20c2743 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -78,6 +78,7 @@ int openmc_finalize() settings::electron_treatment = ElectronTreatment::LED; settings::delayed_photon_scaling = true; settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0}; + settings::time_cutoff = {INFTY, INFTY, INFTY, INFTY}; settings::entropy_on = false; settings::event_based = false; settings::gen_per_batch = 1; diff --git a/src/geometry.cpp b/src/geometry.cpp index 29ca8b1bf8f..5e015c1bc0f 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -110,7 +110,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list) i_cell = *it; // Make sure the search cell is in the same universe. - int i_universe = p.coord(p.n_coord() - 1).universe; + int i_universe = p.lowest_coord().universe; if (model::cells[i_cell]->universe_ != i_universe) continue; @@ -119,7 +119,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list) Direction u {p.u_local()}; auto surf = p.surface(); if (model::cells[i_cell]->contains(r, u, surf)) { - p.coord(p.n_coord() - 1).cell = i_cell; + p.lowest_coord().cell = i_cell; found = true; break; } @@ -145,7 +145,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list) // code below this conditional, we set i_cell back to C_NONE to indicate // that. if (i_cell == C_NONE) { - int i_universe = p.coord(p.n_coord() - 1).universe; + int i_universe = p.lowest_coord().universe; const auto& univ {model::universes[i_universe]}; found = univ->find_cell(p); } @@ -153,7 +153,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list) if (!found) { return found; } - i_cell = p.coord(p.n_coord() - 1).cell; + i_cell = p.lowest_coord().cell; // Announce the cell that the particle is entering. if (found && (settings::verbosity >= 10 || p.trace())) { @@ -289,7 +289,7 @@ bool neighbor_list_find_cell(Particle& p) bool exhaustive_find_cell(Particle& p) { - int i_universe = p.coord(p.n_coord() - 1).universe; + int i_universe = p.lowest_coord().universe; if (i_universe == C_NONE) { p.coord(0).universe = model::root_universe; p.n_coord() = 1; @@ -306,7 +306,7 @@ bool exhaustive_find_cell(Particle& p) void cross_lattice(Particle& p, const BoundaryInfo& boundary) { - auto& coord {p.coord(p.n_coord() - 1)}; + auto& coord {p.lowest_coord()}; auto& lat {*model::lattices[coord.lattice]}; if (settings::verbosity >= 10 || p.trace()) { @@ -344,7 +344,7 @@ void cross_lattice(Particle& p, const BoundaryInfo& boundary) } else { // Find cell in next lattice element. - p.coord(p.n_coord() - 1).universe = lat[coord.lattice_i]; + p.lowest_coord().universe = lat[coord.lattice_i]; bool found = exhaustive_find_cell(p); if (!found) { @@ -472,7 +472,7 @@ extern "C" int openmc_find_cell( return OPENMC_E_GEOMETRY; } - *index = p.coord(p.n_coord() - 1).cell; + *index = p.lowest_coord().cell; *instance = p.cell_instance(); return 0; } diff --git a/src/mgxs.cpp b/src/mgxs.cpp index 472beda3506..9eae5a7da7f 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -5,10 +5,6 @@ #include #include -#ifdef _OPENMP -#include -#endif - #include "xtensor/xadapt.hpp" #include "xtensor/xmath.hpp" #include "xtensor/xsort.hpp" @@ -46,15 +42,6 @@ void Mgxs::init(const std::string& in_name, double in_awr, n_azi = in_azimuthal.size(); polar = in_polar; azimuthal = in_azimuthal; - - // Set the cross section index cache -#ifdef _OPENMP - int n_threads = omp_get_max_threads(); -#else - int n_threads = 1; -#endif - cache.resize(n_threads); - // vector.resize() will value-initialize the members of cache[:] } //============================================================================== @@ -425,18 +412,10 @@ void Mgxs::combine(const vector& micros, const vector& scalars, //============================================================================== -double Mgxs::get_xs( - MgxsType xstype, int gin, const int* gout, const double* mu, const int* dg) +double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu, + const int* dg, int t, int a) { - // This method assumes that the temperature and angle indices are set -#ifdef _OPENMP - int tid = omp_get_thread_num(); - XsData* xs_t = &xs[cache[tid].t]; - int a = cache[tid].a; -#else - XsData* xs_t = &xs[cache[0].t]; - int a = cache[0].a; -#endif + XsData* xs_t = &xs[t]; double val; switch (xstype) { case MgxsType::TOTAL: @@ -538,19 +517,14 @@ double Mgxs::get_xs( //============================================================================== -void Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed) +void Mgxs::sample_fission_energy( + int gin, int& dg, int& gout, uint64_t* seed, int t, int a) { - // This method assumes that the temperature and angle indices are set -#ifdef _OPENMP - int tid = omp_get_thread_num(); -#else - int tid = 0; -#endif - XsData* xs_t = &xs[cache[tid].t]; - double nu_fission = xs_t->nu_fission(cache[tid].a, gin); + XsData* xs_t = &xs[t]; + double nu_fission = xs_t->nu_fission(a, gin); // Find the probability of having a prompt neutron - double prob_prompt = xs_t->prompt_nu_fission(cache[tid].a, gin); + double prob_prompt = xs_t->prompt_nu_fission(a, gin); // sample random numbers double xi_pd = prn(seed) * nu_fission; @@ -566,7 +540,7 @@ void Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed) // sample the outgoing energy group double prob_gout = 0.; for (gout = 0; gout < num_groups; ++gout) { - prob_gout += xs_t->chi_prompt(cache[tid].a, gin, gout); + prob_gout += xs_t->chi_prompt(a, gin, gout); if (xi_gout < prob_gout) break; } @@ -576,7 +550,7 @@ void Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed) // get the delayed group for (dg = 0; dg < num_delayed_groups; ++dg) { - prob_prompt += xs_t->delayed_nu_fission(cache[tid].a, dg, gin); + prob_prompt += xs_t->delayed_nu_fission(a, dg, gin); if (xi_pd < prob_prompt) break; } @@ -587,7 +561,7 @@ void Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed) // sample the outgoing energy group double prob_gout = 0.; for (gout = 0; gout < num_groups; ++gout) { - prob_gout += xs_t->chi_delayed(cache[tid].a, dg, gin, gout); + prob_gout += xs_t->chi_delayed(a, dg, gin, gout); if (xi_gout < prob_gout) break; } @@ -597,35 +571,39 @@ void Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed) //============================================================================== void Mgxs::sample_scatter( - int gin, int& gout, double& mu, double& wgt, uint64_t* seed) + int gin, int& gout, double& mu, double& wgt, uint64_t* seed, int t, int a) { - // This method assumes that the temperature and angle indices are set // Sample the data -#ifdef _OPENMP - int tid = omp_get_thread_num(); -#else - int tid = 0; -#endif - xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt, seed); + xs[t].scatter[a]->sample(gin, gout, mu, wgt, seed); } //============================================================================== void Mgxs::calculate_xs(Particle& p) { - // Set our indices -#ifdef _OPENMP - int tid = omp_get_thread_num(); -#else - int tid = 0; -#endif - set_temperature_index(p.sqrtkT()); - set_angle_index(p.u_local()); - XsData* xs_t = &xs[cache[tid].t]; - p.macro_xs().total = xs_t->total(cache[tid].a, p.g()); - p.macro_xs().absorption = xs_t->absorption(cache[tid].a, p.g()); + // If the material is different, then we need to do a full lookup + if (p.material() != p.mg_xs_cache().material) { + set_temperature_index(p); + set_angle_index(p); + p.mg_xs_cache().material = p.material(); + } else { + // If material is the same, but temperature is different, need to + // find the new temperature index + if (p.sqrtkT() != p.mg_xs_cache().sqrtkT) { + set_temperature_index(p); + } + // If the material is the same, but angle is different, need to + // find the new angle index + if (p.u_local() != p.mg_xs_cache().u) { + set_angle_index(p); + } + } + int temperature = p.mg_xs_cache().t; + int angle = p.mg_xs_cache().a; + p.macro_xs().total = xs[temperature].total(angle, p.g()); + p.macro_xs().absorption = xs[temperature].absorption(angle, p.g()); p.macro_xs().nu_fission = - fissionable ? xs_t->nu_fission(cache[tid].a, p.g()) : 0.; + fissionable ? xs[temperature].nu_fission(angle, p.g()) : 0.; } //============================================================================== @@ -643,32 +621,26 @@ bool Mgxs::equiv(const Mgxs& that) //============================================================================== -void Mgxs::set_temperature_index(double sqrtkT) +int Mgxs::get_temperature_index(double sqrtkT) const { - // See if we need to find the new index -#ifdef _OPENMP - int tid = omp_get_thread_num(); -#else - int tid = 0; -#endif - if (sqrtkT != cache[tid].sqrtkT) { - cache[tid].t = xt::argmin(xt::abs(kTs - sqrtkT * sqrtkT))[0]; - cache[tid].sqrtkT = sqrtkT; - } + return xt::argmin(xt::abs(kTs - sqrtkT * sqrtkT))[0]; } //============================================================================== -void Mgxs::set_angle_index(Direction u) +void Mgxs::set_temperature_index(Particle& p) { - // See if we need to find the new index -#ifdef _OPENMP - int tid = omp_get_thread_num(); -#else - int tid = 0; -#endif - if (!is_isotropic && ((u.x != cache[tid].u) || (u.y != cache[tid].v) || - (u.z != cache[tid].w))) { + p.mg_xs_cache().t = get_temperature_index(p.sqrtkT()); + p.mg_xs_cache().sqrtkT = p.sqrtkT(); +} + +//============================================================================== + +int Mgxs::get_angle_index(const Direction& u) const +{ + if (is_isotropic) { + return 0; + } else { // convert direction to polar and azimuthal angles double my_pol = std::acos(u.z); double my_azi = std::atan2(u.y, u.x); @@ -679,12 +651,18 @@ void Mgxs::set_angle_index(Direction u) delta_angle = 2. * PI / n_azi; int a = std::floor((my_azi + PI) / delta_angle); - cache[tid].a = n_azi * p + a; + return n_azi * p + a; + } +} + +//============================================================================== - // store this direction as the last one used - cache[tid].u = u.x; - cache[tid].v = u.y; - cache[tid].w = u.z; +void Mgxs::set_angle_index(Particle& p) +{ + // See if we need to find the new index + if (!is_isotropic) { + p.mg_xs_cache().a = get_angle_index(p.u_local()); + p.mg_xs_cache().u = p.u_local(); } } diff --git a/src/particle.cpp b/src/particle.cpp index 200bcc6d6dd..ceef073c480 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -65,6 +65,13 @@ double Particle::speed() const return C_LIGHT * std::sqrt(1 - inv_gamma * inv_gamma); } +void Particle::move_distance(double length) +{ + for (int j = 0; j < n_coord(); ++j) { + coord(j).r += length * coord(j).u; + } +} + void Particle::create_secondary( double wgt, Direction u, double E, ParticleType type) { @@ -140,7 +147,7 @@ void Particle::event_calculate_xs() // If the cell hasn't been determined based on the particle's location, // initiate a search for the current cell. This generally happens at the // beginning of the history and again for any secondary particles - if (coord(n_coord() - 1).cell == C_NONE) { + if (lowest_coord().cell == C_NONE) { if (!exhaustive_find_cell(*this)) { mark_as_lost( "Could not find the cell containing particle " + std::to_string(id())); @@ -149,7 +156,7 @@ void Particle::event_calculate_xs() // Set birth cell attribute if (cell_born() == C_NONE) - cell_born() = coord(n_coord() - 1).cell; + cell_born() = lowest_coord().cell; } // Write particle track. @@ -203,11 +210,25 @@ void Particle::event_advance() double distance = std::min(boundary().distance, collision_distance()); // Advance particle in space and time + // Short-term solution until the surface source is revised and we can use + // this->move_distance(distance) for (int j = 0; j < n_coord(); ++j) { coord(j).r += distance * coord(j).u; } this->time() += distance / this->speed(); + // Kill particle if its time exceeds the cutoff + bool hit_time_boundary = false; + double time_cutoff = settings::time_cutoff[static_cast(type())]; + if (time() > time_cutoff) { + double dt = time() - time_cutoff; + time() = time_cutoff; + + double push_back_distance = speed() * dt; + this->move_distance(-push_back_distance); + hit_time_boundary = true; + } + // Score track-length tallies if (!model::active_tracklength_tallies.empty()) { score_tracklength_tally(*this, distance); @@ -223,6 +244,11 @@ void Particle::event_advance() if (!model::active_tallies.empty()) { score_track_derivative(*this, distance); } + + // Set particle weight to zero if it hit the time boundary + if (hit_time_boundary) { + wgt() = 0.0; + } } void Particle::event_cross_surface() @@ -366,7 +392,7 @@ void Particle::event_revive_from_secondary() // Since the birth cell of the particle has not been set we // have to determine it before the energy of the secondary particle can be // removed from the pulse-height of this cell. - if (coord(n_coord() - 1).cell == C_NONE) { + if (lowest_coord().cell == C_NONE) { if (!exhaustive_find_cell(*this)) { mark_as_lost("Could not find the cell containing particle " + std::to_string(id())); @@ -374,7 +400,7 @@ void Particle::event_revive_from_secondary() } // Set birth cell attribute if (cell_born() == C_NONE) - cell_born() = coord(n_coord() - 1).cell; + cell_born() = lowest_coord().cell; } pht_secondary_particles(); } @@ -430,7 +456,7 @@ void Particle::pht_collision_energy() // determine index of cell in pulse_height_cells auto it = std::find(model::pulse_height_cells.begin(), - model::pulse_height_cells.end(), coord(n_coord() - 1).cell); + model::pulse_height_cells.end(), lowest_coord().cell); if (it != model::pulse_height_cells.end()) { int index = std::distance(model::pulse_height_cells.begin(), it); @@ -509,7 +535,7 @@ void Particle::cross_surface() material_last() = material(); sqrtkT_last() = sqrtkT(); // set new cell value - coord(n_coord() - 1).cell = i_cell; + lowest_coord().cell = i_cell; cell_instance() = 0; material() = model::cells[i_cell]->material_[0]; sqrtkT() = model::cells[i_cell]->sqrtkT_[0]; diff --git a/src/physics.cpp b/src/physics.cpp index a2fdaf2e92d..28c316c6878 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -106,7 +106,7 @@ void sample_neutron_reaction(Particle& p) const auto& nuc {data::nuclides[i_nuclide]}; - if (nuc->fissionable_) { + if (nuc->fissionable_ && p.neutron_xs(i_nuclide).fission > 0.0) { auto& rx = sample_fission(i_nuclide, p); if (settings::run_mode == RunMode::EIGENVALUE) { create_fission_sites(p, i_nuclide, rx); diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index ab1fdef028c..43e3cfa7455 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -72,8 +72,8 @@ void sample_reaction(Particle& p) void scatter(Particle& p) { - data::mg.macro_xs_[p.material()].sample_scatter( - p.g_last(), p.g(), p.mu(), p.wgt(), p.current_seed()); + data::mg.macro_xs_[p.material()].sample_scatter(p.g_last(), p.g(), p.mu(), + p.wgt(), p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a); // Rotate the angle p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); @@ -145,7 +145,7 @@ void create_fission_sites(Particle& p) int dg; int gout; data::mg.macro_xs_[p.material()].sample_fission_energy( - p.g(), dg, gout, p.current_seed()); + p.g(), dg, gout, p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a); // Store the energy and delayed groups on the fission bank site.E = gout; diff --git a/src/plot.cpp b/src/plot.cpp index 3bcae05ac53..e2f07c47cf0 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -58,7 +58,7 @@ void IdData::set_value(size_t y, size_t x, const Particle& p, int level) } // set material data - Cell* c = model::cells.at(p.coord(p.n_coord() - 1).cell).get(); + Cell* c = model::cells.at(p.lowest_coord().cell).get(); if (p.material() == MATERIAL_VOID) { data_(y, x, 2) = MATERIAL_VOID; return; @@ -79,7 +79,7 @@ PropertyData::PropertyData(size_t h_res, size_t v_res) void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level) { - Cell* c = model::cells.at(p.coord(p.n_coord() - 1).cell).get(); + Cell* c = model::cells.at(p.lowest_coord().cell).get(); data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN; if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) { Material* m = model::materials.at(p.material()).get(); @@ -1327,9 +1327,8 @@ void ProjectionPlot::create_output() const // edges on the model boundary for the same cell. if (first_inside_model) { this_line_segments[tid][horiz].emplace_back( - color_by_ == PlotColorBy::mats - ? p.material() - : p.coord(p.n_coord() - 1).cell, + color_by_ == PlotColorBy::mats ? p.material() + : p.lowest_coord().cell, 0.0, first_surface); first_inside_model = false; } @@ -1339,7 +1338,7 @@ void ProjectionPlot::create_output() const auto dist = distance_to_boundary(p); this_line_segments[tid][horiz].emplace_back( color_by_ == PlotColorBy::mats ? p.material() - : p.coord(p.n_coord() - 1).cell, + : p.lowest_coord().cell, dist.distance, std::abs(dist.surface_index)); // Advance particle diff --git a/src/settings.cpp b/src/settings.cpp index 1c039e9ba90..ac56320e8e6 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -22,6 +22,7 @@ #include "openmc/mesh.h" #include "openmc/message_passing.h" #include "openmc/output.h" +#include "openmc/plot.h" #include "openmc/random_lcg.h" #include "openmc/simulation.h" #include "openmc/source.h" @@ -95,6 +96,7 @@ int64_t max_particles_in_flight {100000}; ElectronTreatment electron_treatment {ElectronTreatment::TTB}; array energy_cutoff {0.0, 1000.0, 0.0, 0.0}; +array time_cutoff {INFTY, INFTY, INFTY, INFTY}; int legendre_to_tabular_points {C_NONE}; int max_order {0}; int n_log_bins {8000}; @@ -390,6 +392,12 @@ void read_settings_xml(pugi::xml_node root) } } + // Copy plotting random number seed if specified + if (check_for_node(root, "plot_seed")) { + auto seed = std::stoll(get_node_value(root, "plot_seed")); + model::plotter_seed = seed; + } + // Copy random number seed if specified if (check_for_node(root, "seed")) { auto seed = std::stoll(get_node_value(root, "seed")); @@ -533,6 +541,18 @@ void read_settings_xml(pugi::xml_node root) energy_cutoff[3] = std::stod(get_node_value(node_cutoff, "energy_positron")); } + if (check_for_node(node_cutoff, "time_neutron")) { + time_cutoff[0] = std::stod(get_node_value(node_cutoff, "time_neutron")); + } + if (check_for_node(node_cutoff, "time_photon")) { + time_cutoff[1] = std::stod(get_node_value(node_cutoff, "time_photon")); + } + if (check_for_node(node_cutoff, "time_electron")) { + time_cutoff[2] = std::stod(get_node_value(node_cutoff, "time_electron")); + } + if (check_for_node(node_cutoff, "time_positron")) { + time_cutoff[3] = std::stod(get_node_value(node_cutoff, "time_positron")); + } } // Particle trace diff --git a/src/tallies/filter_cell_instance.cpp b/src/tallies/filter_cell_instance.cpp index 59f28c7d66f..3e78a5bbed9 100644 --- a/src/tallies/filter_cell_instance.cpp +++ b/src/tallies/filter_cell_instance.cpp @@ -72,7 +72,7 @@ void CellInstanceFilter::set_cell_instances(gsl::span instances) void CellInstanceFilter::get_all_bins( const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - gsl::index index_cell = p.coord(p.n_coord() - 1).cell; + gsl::index index_cell = p.lowest_coord().cell; gsl::index instance = p.cell_instance(); if (cells_.count(index_cell) > 0) { diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 47fe4d1be7b..e798161ec2f 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -907,7 +907,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, int m; switch (score_bin) { - // clang-format off + // clang-format off case N_GAMMA: m = 0; break; case N_P: m = 1; break; case N_A: m = 2; break; @@ -1595,10 +1595,13 @@ void score_general_mg(Particle& p, int i_tally, int start_index, auto& macro_xs = data::mg.macro_xs_[p.material()]; // Find the temperature and angle indices of interest - macro_xs.set_angle_index(p_u); + int macro_t = p.mg_xs_cache().t; + int macro_a = macro_xs.get_angle_index(p_u); + int nuc_t = 0; + int nuc_a = 0; if (i_nuclide >= 0) { - nuc_xs.set_temperature_index(p.sqrtkT()); - nuc_xs.set_angle_index(p_u); + nuc_t = nuc_xs.get_temperature_index(p.sqrtkT()); + nuc_a = nuc_xs.get_angle_index(p_u); } for (auto i = 0; i < tally.scores_.size(); ++i) { @@ -1626,12 +1629,14 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // use the weight of the particle entering the collision as the score score = flux * p.wgt_last(); if (i_nuclide >= 0) { - score *= atom_density * nuc_xs.get_xs(MgxsType::TOTAL, p_g) / - macro_xs.get_xs(MgxsType::TOTAL, p_g); + score *= atom_density * + nuc_xs.get_xs(MgxsType::TOTAL, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::TOTAL, p_g, macro_t, macro_a); } } else { if (i_nuclide >= 0) { - score = atom_density * flux * nuc_xs.get_xs(MgxsType::TOTAL, p_g); + score = atom_density * flux * + nuc_xs.get_xs(MgxsType::TOTAL, p_g, nuc_t, nuc_a); } else { score = p.macro_xs().total * flux; } @@ -1646,17 +1651,21 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // to count 'events' exactly for the inverse velocity score = flux * p.wgt_last(); if (i_nuclide >= 0) { - score *= nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) / - macro_xs.get_xs(MgxsType::TOTAL, p_g); + score *= + nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::TOTAL, p_g, macro_t, macro_a); } else { - score *= macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) / - macro_xs.get_xs(MgxsType::TOTAL, p_g); + score *= + macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g, macro_t, macro_a) / + macro_xs.get_xs(MgxsType::TOTAL, p_g, macro_t, macro_a); } } else { if (i_nuclide >= 0) { - score = flux * nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g); + score = + flux * nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g, nuc_t, nuc_a); } else { - score = flux * macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g); + score = flux * macro_xs.get_xs( + MgxsType::INVERSE_VELOCITY, p_g, macro_t, macro_a); } } break; @@ -1672,18 +1681,18 @@ void score_general_mg(Particle& p, int i_tally, int start_index, if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last(), &p.g(), - &p.mu(), nullptr) / + &p.mu(), nullptr, nuc_t, nuc_a) / macro_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last(), &p.g(), - &p.mu(), nullptr); + &p.mu(), nullptr, macro_t, macro_a); } } else { if (i_nuclide >= 0) { - score = - atom_density * flux * - nuc_xs.get_xs(MgxsType::SCATTER, p_g, nullptr, &p.mu(), nullptr); + score = atom_density * flux * + nuc_xs.get_xs(MgxsType::SCATTER, p_g, nullptr, &p.mu(), + nullptr, nuc_t, nuc_a); } else { - score = flux * macro_xs.get_xs( - MgxsType::SCATTER, p_g, nullptr, &p.mu(), nullptr); + score = flux * macro_xs.get_xs(MgxsType::SCATTER, p_g, nullptr, + &p.mu(), nullptr, macro_t, macro_a); } } break; @@ -1703,16 +1712,17 @@ void score_general_mg(Particle& p, int i_tally, int start_index, if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::NU_SCATTER_FMU, p.g_last(), &p.g(), - &p.mu(), nullptr) / + &p.mu(), nullptr, nuc_t, nuc_a) / macro_xs.get_xs(MgxsType::NU_SCATTER_FMU, p.g_last(), &p.g(), - &p.mu(), nullptr); + &p.mu(), nullptr, macro_t, macro_a); } } else { if (i_nuclide >= 0) { - score = - atom_density * flux * nuc_xs.get_xs(MgxsType::NU_SCATTER, p_g); + score = atom_density * flux * + nuc_xs.get_xs(MgxsType::NU_SCATTER, p_g, nuc_t, nuc_a); } else { - score = flux * macro_xs.get_xs(MgxsType::NU_SCATTER, p_g); + score = + flux * macro_xs.get_xs(MgxsType::NU_SCATTER, p_g, macro_t, macro_a); } } break; @@ -1732,13 +1742,14 @@ void score_general_mg(Particle& p, int i_tally, int start_index, score = p.wgt_last() * flux; } if (i_nuclide >= 0) { - score *= atom_density * nuc_xs.get_xs(MgxsType::ABSORPTION, p_g) / - macro_xs.get_xs(MgxsType::ABSORPTION, p_g); + score *= atom_density * + nuc_xs.get_xs(MgxsType::ABSORPTION, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a); } } else { if (i_nuclide >= 0) { - score = - atom_density * flux * nuc_xs.get_xs(MgxsType::ABSORPTION, p_g); + score = atom_density * flux * + nuc_xs.get_xs(MgxsType::ABSORPTION, p_g, nuc_t, nuc_a); } else { score = p.macro_xs().absorption * flux; } @@ -1762,17 +1773,20 @@ void score_general_mg(Particle& p, int i_tally, int start_index, score = p.wgt_last() * flux; } if (i_nuclide >= 0) { - score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / - macro_xs.get_xs(MgxsType::ABSORPTION, p_g); + score *= atom_density * + nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a); } else { - score *= macro_xs.get_xs(MgxsType::FISSION, p_g) / - macro_xs.get_xs(MgxsType::ABSORPTION, p_g); + score *= macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a) / + macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a); } } else { if (i_nuclide >= 0) { - score = atom_density * flux * nuc_xs.get_xs(MgxsType::FISSION, p_g); + score = atom_density * flux * + nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a); } else { - score = flux * macro_xs.get_xs(MgxsType::FISSION, p_g); + score = + flux * macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a); } } break; @@ -1793,11 +1807,14 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // nu-fission score = wgt_absorb * flux; if (i_nuclide >= 0) { - score *= atom_density * nuc_xs.get_xs(MgxsType::NU_FISSION, p_g) / - macro_xs.get_xs(MgxsType::ABSORPTION, p_g); + score *= + atom_density * + nuc_xs.get_xs(MgxsType::NU_FISSION, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a); } else { - score *= macro_xs.get_xs(MgxsType::NU_FISSION, p_g) / - macro_xs.get_xs(MgxsType::ABSORPTION, p_g); + score *= + macro_xs.get_xs(MgxsType::NU_FISSION, p_g, macro_t, macro_a) / + macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a); } } else { // Skip any non-fission events @@ -1810,16 +1827,18 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // score. score = simulation::keff * p.wgt_bank() * flux; if (i_nuclide >= 0) { - score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / - macro_xs.get_xs(MgxsType::FISSION, p_g); + score *= atom_density * + nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a); } } } else { if (i_nuclide >= 0) { - score = - atom_density * flux * nuc_xs.get_xs(MgxsType::NU_FISSION, p_g); + score = atom_density * flux * + nuc_xs.get_xs(MgxsType::NU_FISSION, p_g, nuc_t, nuc_a); } else { - score = flux * macro_xs.get_xs(MgxsType::NU_FISSION, p_g); + score = + flux * macro_xs.get_xs(MgxsType::NU_FISSION, p_g, macro_t, macro_a); } } break; @@ -1840,12 +1859,15 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // prompt-nu-fission score = wgt_absorb * flux; if (i_nuclide >= 0) { - score *= atom_density * - nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g) / - macro_xs.get_xs(MgxsType::ABSORPTION, p_g); + score *= + atom_density * + nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a); } else { - score *= macro_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g) / - macro_xs.get_xs(MgxsType::ABSORPTION, p_g); + score *= + macro_xs.get_xs( + MgxsType::PROMPT_NU_FISSION, p_g, macro_t, macro_a) / + macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a); } } else { // Skip any non-fission events @@ -1861,16 +1883,18 @@ void score_general_mg(Particle& p, int i_tally, int start_index, auto prompt_frac = 1. - n_delayed / static_cast(p.n_bank()); score = simulation::keff * p.wgt_bank() * prompt_frac * flux; if (i_nuclide >= 0) { - score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / - macro_xs.get_xs(MgxsType::FISSION, p_g); + score *= atom_density * + nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a); } } } else { if (i_nuclide >= 0) { score = atom_density * flux * - nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g); + nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g, nuc_t, nuc_a); } else { - score = flux * macro_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g); + score = flux * macro_xs.get_xs( + MgxsType::PROMPT_NU_FISSION, p_g, macro_t, macro_a); } } break; @@ -1889,7 +1913,8 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - double abs_xs = macro_xs.get_xs(MgxsType::ABSORPTION, p_g); + double abs_xs = + macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a); if (abs_xs > 0.) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; @@ -1902,11 +1927,11 @@ void score_general_mg(Particle& p, int i_tally, int start_index, score = wgt_absorb * flux; if (i_nuclide >= 0) { score *= nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, - nullptr, nullptr, &d) / + nullptr, nullptr, &d, nuc_t, nuc_a) / abs_xs; } else { score *= macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, - nullptr, nullptr, &d) / + nullptr, nullptr, &d, macro_t, macro_a) / abs_xs; } score_fission_delayed_dg( @@ -1919,11 +1944,13 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // delayed-nu-fission xs to the absorption xs score = wgt_absorb * flux; if (i_nuclide >= 0) { - score *= - nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g) / abs_xs; + score *= nuc_xs.get_xs( + MgxsType::DELAYED_NU_FISSION, p_g, nuc_t, nuc_a) / + abs_xs; } else { - score *= - macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g) / abs_xs; + score *= macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, + macro_t, macro_a) / + abs_xs; } } } @@ -1948,8 +1975,10 @@ void score_general_mg(Particle& p, int i_tally, int start_index, score = simulation::keff * p.wgt_bank() / p.n_bank() * p.n_delayed_bank(d - 1) * flux; if (i_nuclide >= 0) { - score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / - macro_xs.get_xs(MgxsType::FISSION, p_g); + score *= + atom_density * + nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a); } score_fission_delayed_dg( i_tally, d_bin, score, score_index, p.filter_matches()); @@ -1962,8 +1991,10 @@ void score_general_mg(Particle& p, int i_tally, int start_index, score = simulation::keff * p.wgt_bank() / p.n_bank() * n_delayed * flux; if (i_nuclide >= 0) { - score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / - macro_xs.get_xs(MgxsType::FISSION, p_g); + score *= + atom_density * + nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a); } } } @@ -1978,10 +2009,10 @@ void score_general_mg(Particle& p, int i_tally, int start_index, if (i_nuclide >= 0) { score = flux * atom_density * nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, - nullptr, &d); + nullptr, &d, nuc_t, nuc_a); } else { score = flux * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, - nullptr, nullptr, &d); + nullptr, nullptr, &d, macro_t, macro_a); } score_fission_delayed_dg( i_tally, d_bin, score, score_index, p.filter_matches()); @@ -1989,10 +2020,12 @@ void score_general_mg(Particle& p, int i_tally, int start_index, continue; } else { if (i_nuclide >= 0) { - score = flux * atom_density * - nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g); + score = + flux * atom_density * + nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nuc_t, nuc_a); } else { - score = flux * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g); + score = flux * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, + macro_t, macro_a); } } } @@ -2004,7 +2037,8 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - double abs_xs = macro_xs.get_xs(MgxsType::ABSORPTION, p_g); + double abs_xs = + macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a); if (abs_xs > 0) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; @@ -2016,16 +2050,16 @@ void score_general_mg(Particle& p, int i_tally, int start_index, auto d = filt.groups()[d_bin] - 1; score = wgt_absorb * flux; if (i_nuclide >= 0) { - score *= nuc_xs.get_xs( - MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * + score *= nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, + nullptr, &d, nuc_t, nuc_a) * nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, - nullptr, nullptr, &d) / + nullptr, nullptr, &d, nuc_t, nuc_a) / abs_xs; } else { - score *= macro_xs.get_xs( - MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * + score *= macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, + nullptr, &d, macro_t, macro_a) * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, - nullptr, nullptr, &d) / + nullptr, nullptr, &d, macro_t, macro_a) / abs_xs; } score_fission_delayed_dg( @@ -2041,17 +2075,17 @@ void score_general_mg(Particle& p, int i_tally, int start_index, for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) { if (i_nuclide >= 0) { score += wgt_absorb * flux * - nuc_xs.get_xs( - MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * + nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, + nullptr, &d, nuc_t, nuc_a) * nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, - nullptr, nullptr, &d) / + nullptr, nullptr, &d, nuc_t, nuc_a) / abs_xs; } else { score += wgt_absorb * flux * - macro_xs.get_xs( - MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * + macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, + nullptr, &d, macro_t, macro_a) * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, - nullptr, nullptr, &d) / + nullptr, nullptr, &d, macro_t, macro_a) / abs_xs; } } @@ -2074,15 +2108,16 @@ void score_general_mg(Particle& p, int i_tally, int start_index, auto d = bank.delayed_group - 1; if (d != -1) { if (i_nuclide >= 0) { - score += simulation::keff * atom_density * bank.wgt * flux * - nuc_xs.get_xs( - MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * - nuc_xs.get_xs(MgxsType::FISSION, p_g) / - macro_xs.get_xs(MgxsType::FISSION, p_g); + score += + simulation::keff * atom_density * bank.wgt * flux * + nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d, + nuc_t, nuc_a) * + nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a); } else { score += simulation::keff * bank.wgt * flux * - macro_xs.get_xs( - MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d); + macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, + nullptr, &d, macro_t, macro_a); } if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; @@ -2112,17 +2147,17 @@ void score_general_mg(Particle& p, int i_tally, int start_index, for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin] - 1; if (i_nuclide >= 0) { - score = - atom_density * flux * - nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * - nuc_xs.get_xs( - MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); + score = atom_density * flux * + nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, + &d, nuc_t, nuc_a) * + nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, + nullptr, &d, nuc_t, nuc_a); } else { score = flux * - macro_xs.get_xs( - MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * + macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, + nullptr, &d, macro_t, macro_a) * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, - nullptr, nullptr, &d); + nullptr, nullptr, &d, macro_t, macro_a); } score_fission_delayed_dg( i_tally, d_bin, score, score_index, p.filter_matches()); @@ -2132,17 +2167,17 @@ void score_general_mg(Particle& p, int i_tally, int start_index, score = 0.; for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) { if (i_nuclide >= 0) { - score += - atom_density * flux * - nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * - nuc_xs.get_xs( - MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); + score += atom_density * flux * + nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, + nullptr, &d, nuc_t, nuc_a) * + nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, + nullptr, &d, nuc_t, nuc_a); } else { score += flux * - macro_xs.get_xs( - MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * + macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, + nullptr, &d, macro_t, macro_a) * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, - nullptr, nullptr, &d); + nullptr, nullptr, &d, macro_t, macro_a); } } } @@ -2166,18 +2201,21 @@ void score_general_mg(Particle& p, int i_tally, int start_index, score = p.wgt_last() * flux; } if (i_nuclide >= 0) { - score *= atom_density * nuc_xs.get_xs(MgxsType::KAPPA_FISSION, p_g) / - macro_xs.get_xs(MgxsType::ABSORPTION, p_g); + score *= atom_density * + nuc_xs.get_xs(MgxsType::KAPPA_FISSION, p_g, nuc_t, nuc_a) / + macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a); } else { - score *= macro_xs.get_xs(MgxsType::KAPPA_FISSION, p_g) / - macro_xs.get_xs(MgxsType::ABSORPTION, p_g); + score *= + macro_xs.get_xs(MgxsType::KAPPA_FISSION, p_g, macro_t, macro_a) / + macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a); } } else { if (i_nuclide >= 0) { - score = - atom_density * flux * nuc_xs.get_xs(MgxsType::KAPPA_FISSION, p_g); + score = atom_density * flux * + nuc_xs.get_xs(MgxsType::KAPPA_FISSION, p_g, nuc_t, nuc_a); } else { - score = flux * macro_xs.get_xs(MgxsType::KAPPA_FISSION, p_g); + score = flux * macro_xs.get_xs( + MgxsType::KAPPA_FISSION, p_g, macro_t, macro_a); } } break; diff --git a/src/universe.cpp b/src/universe.cpp index f8f9a82d896..67e1d1354c3 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -41,18 +41,17 @@ bool Universe::find_cell(Particle& p) const const auto& cells { !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())}; - for (auto it = cells.begin(); it != cells.end(); it++) { - int32_t i_cell = *it; - int32_t i_univ = p.coord(p.n_coord() - 1).universe; + Position r {p.r_local()}; + Position u {p.u_local()}; + auto surf = p.surface(); + int32_t i_univ = p.lowest_coord().universe; + + for (auto i_cell : cells) { if (model::cells[i_cell]->universe_ != i_univ) continue; - - // Check if this cell contains the particle; - Position r {p.r_local()}; - Direction u {p.u_local()}; - auto surf = p.surface(); + // Check if this cell contains the particle if (model::cells[i_cell]->contains(r, u, surf)) { - p.coord(p.n_coord() - 1).cell = i_cell; + p.lowest_coord().cell = i_cell; return true; } } diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 0b1ea357505..8949d5c20d3 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -534,8 +534,21 @@ int openmc_calculate_volumes() // Display domain volumes for (int j = 0; j < vol_calc.domain_ids_.size(); j++) { - write_message(4, "{}{}: {} +/- {} cm^3", domain_type, - vol_calc.domain_ids_[j], results[j].volume[0], results[j].volume[1]); + std::string region_name {""}; + if (vol_calc.domain_type_ == VolumeCalculation::TallyDomain::CELL) { + int cell_idx = model::cell_map[vol_calc.domain_ids_[j]]; + region_name = model::cells[cell_idx]->name(); + } else if (vol_calc.domain_type_ == + VolumeCalculation::TallyDomain::MATERIAL) { + int mat_idx = model::material_map[vol_calc.domain_ids_[j]]; + region_name = model::materials[mat_idx]->name(); + } + if (region_name.size()) + region_name.insert(0, " "); // prepend space for formatting + + write_message(4, "{}{}{}: {} +/- {} cm^3", domain_type, + vol_calc.domain_ids_[j], region_name, results[j].volume[0], + results[j].volume[1]); } // Write volumes to HDF5 file diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index a1f8731f8e5..fe2fcb5cb7a 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -1024,6 +1024,86 @@ extern "C" int openmc_weight_windows_set_bounds(int32_t index, return 0; } +extern "C" int openmc_weight_windows_get_survival_ratio( + int32_t index, double* ratio) +{ + if (int err = verify_ww_index(index)) + return err; + const auto& wws = variance_reduction::weight_windows[index]; + *ratio = wws->survival_ratio(); + return 0; +} + +extern "C" int openmc_weight_windows_set_survival_ratio( + int32_t index, double ratio) +{ + if (int err = verify_ww_index(index)) + return err; + const auto& wws = variance_reduction::weight_windows[index]; + wws->survival_ratio() = ratio; + std::cout << "Survival ratio: " << wws->survival_ratio() << std::endl; + return 0; +} + +extern "C" int openmc_weight_windows_get_max_lower_bound_ratio( + int32_t index, double* lb_ratio) +{ + if (int err = verify_ww_index(index)) + return err; + const auto& wws = variance_reduction::weight_windows[index]; + *lb_ratio = wws->max_lower_bound_ratio(); + return 0; +} + +extern "C" int openmc_weight_windows_set_max_lower_bound_ratio( + int32_t index, double lb_ratio) +{ + if (int err = verify_ww_index(index)) + return err; + const auto& wws = variance_reduction::weight_windows[index]; + wws->max_lower_bound_ratio() = lb_ratio; + return 0; +} + +extern "C" int openmc_weight_windows_get_weight_cutoff( + int32_t index, double* cutoff) +{ + if (int err = verify_ww_index(index)) + return err; + const auto& wws = variance_reduction::weight_windows[index]; + *cutoff = wws->weight_cutoff(); + return 0; +} + +extern "C" int openmc_weight_windows_set_weight_cutoff( + int32_t index, double cutoff) +{ + if (int err = verify_ww_index(index)) + return err; + const auto& wws = variance_reduction::weight_windows[index]; + wws->weight_cutoff() = cutoff; + return 0; +} + +extern "C" int openmc_weight_windows_get_max_split( + int32_t index, int* max_split) +{ + if (int err = verify_ww_index(index)) + return err; + const auto& wws = variance_reduction::weight_windows[index]; + *max_split = wws->max_split(); + return 0; +} + +extern "C" int openmc_weight_windows_set_max_split(int32_t index, int max_split) +{ + if (int err = verify_ww_index(index)) + return err; + const auto& wws = variance_reduction::weight_windows[index]; + wws->max_split() = max_split; + return 0; +} + extern "C" int openmc_extend_weight_windows( int32_t n, int32_t* index_start, int32_t* index_end) { diff --git a/tests/regression_tests/filter_mesh/test.py b/tests/regression_tests/filter_mesh/test.py index b047dd252e1..62f0ed015d3 100644 --- a/tests/regression_tests/filter_mesh/test.py +++ b/tests/regression_tests/filter_mesh/test.py @@ -60,11 +60,11 @@ def model(): recti_mesh_exp_vols = np.multiply.outer(dxdy, dz) np.testing.assert_allclose(recti_mesh.volumes, recti_mesh_exp_vols) - cyl_mesh = openmc.CylindricalMesh() - cyl_mesh.r_grid = np.linspace(0, 7.5, 18) - cyl_mesh.phi_grid = np.linspace(0, 2*pi, 19) - cyl_mesh.z_grid = np.linspace(-7.5, 7.5, 17) - + cyl_mesh = openmc.CylindricalMesh( + r_grid=np.linspace(0, 7.5, 18), + phi_grid=np.linspace(0, 2*pi, 19), + z_grid=np.linspace(-7.5, 7.5, 17), + ) dr = 0.5 * np.diff(np.linspace(0, 7.5, 18)**2) dp = np.full(cyl_mesh.dimension[1], 2*pi / 18) dz = np.full(cyl_mesh.dimension[2], 15 / 16) @@ -72,10 +72,11 @@ def model(): cyl_mesh_exp_vols = np.multiply.outer(drdp, dz) np.testing.assert_allclose(cyl_mesh.volumes, cyl_mesh_exp_vols) - sph_mesh = openmc.SphericalMesh() - sph_mesh.r_grid = np.linspace(0, 7.5, 18) - sph_mesh.theta_grid = np.linspace(0, pi, 9) - sph_mesh.phi_grid = np.linspace(0, 2*pi, 19) + sph_mesh = openmc.SphericalMesh( + r_grid=np.linspace(0, 7.5, 18), + theta_grid=np.linspace(0, pi, 9), + phi_grid=np.linspace(0, 2*pi, 19) + ) dr = np.diff(np.linspace(0, 7.5, 18)**3) / 3 dt = np.diff(-np.cos(np.linspace(0, pi, 9))) dp = np.full(sph_mesh.dimension[2], 2*pi / 18) diff --git a/tests/regression_tests/mg_temperature_multi/__init__.py b/tests/regression_tests/mg_temperature_multi/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/mg_temperature_multi/inputs_true.dat b/tests/regression_tests/mg_temperature_multi/inputs_true.dat new file mode 100644 index 00000000000..b2ad5063c2b --- /dev/null +++ b/tests/regression_tests/mg_temperature_multi/inputs_true.dat @@ -0,0 +1,53 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 10 + 5 + + + -0.63 -0.63 -1 0.63 0.63 1 + + + multi-group + + + + 1 + + + 2 + + + 1 + flux + + + 2 + flux + + + diff --git a/tests/regression_tests/mg_temperature_multi/results_true.dat b/tests/regression_tests/mg_temperature_multi/results_true.dat new file mode 100644 index 00000000000..80c2f562250 --- /dev/null +++ b/tests/regression_tests/mg_temperature_multi/results_true.dat @@ -0,0 +1,8 @@ +k-combined: +1.337115E+00 7.221249E-03 +tally 1: +2.549066E+01 +1.299987E+02 +tally 2: +9.276778E+01 +1.721791E+03 diff --git a/tests/regression_tests/mg_temperature_multi/test.py b/tests/regression_tests/mg_temperature_multi/test.py new file mode 100755 index 00000000000..2090f1b0c2c --- /dev/null +++ b/tests/regression_tests/mg_temperature_multi/test.py @@ -0,0 +1,166 @@ +import os + +import numpy as np +import openmc +import openmc.mgxs + +from tests.testing_harness import PyAPITestHarness + + +def create_library(): + # Instantiate the energy group data + egroups = [1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6, 1.0e6, 20.0e6] + groups = openmc.mgxs.EnergyGroups(egroups) + + # Instantiate the 7-group (C5G7) cross section data + uo2_xsdata = openmc.XSdata('UO2', groups, temperatures=[294.0, 600.0]) + uo2_xsdata.order = 0 + scatter_matrix = np.array([[ + [0.1275370, 0.0423780, 0.0000094, 0.0000000, 0.0000000, 0.0000000, 0.0000000], + [0.0000000, 0.3244560, 0.0016314, 0.0000000, 0.0000000, 0.0000000, 0.0000000], + [0.0000000, 0.0000000, 0.4509400, 0.0026792, 0.0000000, 0.0000000, 0.0000000], + [0.0000000, 0.0000000, 0.0000000, 0.4525650, 0.0055664, 0.0000000, 0.0000000], + [0.0000000, 0.0000000, 0.0000000, 0.0001253, 0.2714010, 0.0102550, 0.0000000], + [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0012968, 0.2658020, 0.0168090], + [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0085458, 0.2730800] + ]]) + scatter_matrix = np.rollaxis(scatter_matrix, 0, 3) + + # Original C5G7 data + uo2_xsdata.set_total([0.1779492, 0.3298048, 0.4803882, 0.5543674, 0.3118013, 0.3951678, 0.5644058], temperature=294.0) + uo2_xsdata.set_absorption([8.0248E-03, 3.7174E-03, 2.6769E-02, 9.6236E-02, 3.0020E-02, 1.1126E-01, 2.8278E-01], temperature=294.0) + uo2_xsdata.set_scatter_matrix(scatter_matrix, temperature=294.0) + uo2_xsdata.set_fission([7.21206E-03, 8.19301E-04, 6.45320E-03, 1.85648E-02, 1.78084E-02, 8.30348E-02, 2.16004E-01], temperature=294.0) + uo2_xsdata.set_nu_fission([2.005998E-02, 2.027303E-03, 1.570599E-02, 4.518301E-02, 4.334208E-02, 2.020901E-01, 5.257105E-01], temperature=294.0) + uo2_xsdata.set_chi([5.8791E-01, 4.1176E-01, 3.3906E-04, 1.1761E-07, 0.0000E+00, 0.0000E+00, 0.0000E+00], temperature=294.0) + + # Altered C5G7 data (permuted Chi) + uo2_xsdata.set_total([0.1779492, 0.3298048, 0.4803882, 0.5543674, 0.3118013, 0.3951678, 0.5644058], temperature=600.0) + uo2_xsdata.set_absorption([8.0248E-03, 3.7174E-03, 2.6769E-02, 9.6236E-02, 3.0020E-02, 1.1126E-01, 2.8278E-01], temperature=600.0) + uo2_xsdata.set_scatter_matrix(scatter_matrix, temperature=600.0) + uo2_xsdata.set_fission([7.21206E-03, 8.19301E-04, 6.45320E-03, 1.85648E-02, 1.78084E-02, 8.30348E-02, 2.16004E-01], temperature=600.0) + uo2_xsdata.set_nu_fission([2.005998E-02, 2.027303E-03, 1.570599E-02, 4.518301E-02, 4.334208E-02, 2.020901E-01, 5.257105E-01], temperature=600.0) + uo2_xsdata.set_chi([4.1176E-01, 5.8791E-01, 3.3906E-04, 1.1761E-07, 0.0000E+00, 0.0000E+00, 0.0000E+00], temperature=600.0) + + h2o_xsdata = openmc.XSdata('LWTR', groups) + h2o_xsdata.order = 0 + h2o_xsdata.set_total([0.15920605, 0.412969593, 0.59030986, 0.58435, + 0.718, 1.2544497, 2.650379]) + h2o_xsdata.set_absorption([6.0105E-04, 1.5793E-05, 3.3716E-04, + 1.9406E-03, 5.7416E-03, 1.5001E-02, + 3.7239E-02]) + scatter_matrix = np.array([[ + [0.0444777, 0.1134000, 0.0007235, 0.0000037, 0.0000001, 0.0000000, 0.0000000], + [0.0000000, 0.2823340, 0.1299400, 0.0006234, 0.0000480, 0.0000074, 0.0000010], + [0.0000000, 0.0000000, 0.3452560, 0.2245700, 0.0169990, 0.0026443, 0.0005034], + [0.0000000, 0.0000000, 0.0000000, 0.0910284, 0.4155100, 0.0637320, 0.0121390], + [0.0000000, 0.0000000, 0.0000000, 0.0000714, 0.1391380, 0.5118200, 0.0612290], + [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200], + [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000] + ]]) + scatter_matrix = np.rollaxis(scatter_matrix, 0, 3) + h2o_xsdata.set_scatter_matrix(scatter_matrix) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata]) + mg_cross_sections_file.export_to_hdf5() + + +class MGXSTestHarness(PyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +def test_mg_temperature_multi(): + ############################################################################### + # Create multigroup data + create_library() + + ############################################################################### + # Create materials for the problem + + # Instantiate some Macroscopic Data + uo2_data = openmc.Macroscopic('UO2') + h2o_data = openmc.Macroscopic('LWTR') + + # Instantiate some Materials and register the appropriate Macroscopic objects + uo2 = openmc.Material(name='UO2 fuel') + uo2.set_density('macro', 1.0) + uo2.add_macroscopic(uo2_data) + + water = openmc.Material(name='Water') + water.set_density('macro', 1.0) + water.add_macroscopic(h2o_data) + + # Instantiate a Materials collection and export to XML + materials = openmc.Materials([uo2, water]) + materials.cross_sections = "mgxs.h5" + + ############################################################################### + # Define problem geometry + + # Create a surface for the fuel outer radius + fuel_ir = openmc.ZCylinder(r=0.25, name='Fuel IR') + fuel_or = openmc.ZCylinder(r=0.54, name='Fuel OR') + + # Create a region represented as the inside of a rectangular prism + pitch = 1.26 + box = openmc.rectangular_prism(pitch, pitch, boundary_type='reflective') + + # Instantiate Cells + fuel_inner = openmc.Cell(fill=uo2, region=-fuel_ir, name='fuel inner') + fuel_inner.temperature = 600.0 + fuel_outer = openmc.Cell(fill=uo2, region=+fuel_ir & -fuel_or, name='fuel outer') + fuel_outer.temperature = 294.0 + moderator = openmc.Cell(fill=water, region=+fuel_or & box, name='moderator') + + # Create a geometry with the two cells and export to XML + geometry = openmc.Geometry([fuel_inner, fuel_outer, moderator]) + + ############################################################################### + # Define problem settings + + # Instantiate a Settings object, set all runtime parameters, and export to XML + settings = openmc.Settings() + settings.energy_mode = "multi-group" + settings.batches = 10 + settings.inactive = 5 + settings.particles = 1000 + + # Create an initial uniform spatial source distribution over fissionable zones + lower_left = (-pitch/2, -pitch/2, -1) + upper_right = (pitch/2, pitch/2, 1) + uniform_dist = openmc.stats.Box(lower_left, upper_right, only_fissionable=True) + settings.source = openmc.IndependentSource(space=uniform_dist) + + ############################################################################### + # Define tallies + + # Instantiate the energy group data + egroups = [1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6, 1.0e6, 20.0e6] + + inner_filter = openmc.CellFilter(fuel_inner) + outer_filter = openmc.CellFilter(fuel_outer) + energy_filter = openmc.EnergyFilter(egroups) + + inner_tally = openmc.Tally(name="inner tally") + inner_tally.filters = [energy_filter] + inner_tally.filters = [inner_filter] + inner_tally.scores = ['flux'] + + outer_tally = openmc.Tally(name="outer tally") + outer_tally.filters = [energy_filter] + outer_tally.filters = [outer_filter] + outer_tally.scores = ['flux'] + + # Instantiate a Tallies collection and export to XML + tallies = openmc.Tallies([inner_tally, outer_tally]) + + # Generate model and run test + model = openmc.Model(geometry, materials, settings, tallies) + + harness = MGXSTestHarness('statepoint.10.h5', model=model) + harness.main() diff --git a/tests/regression_tests/pulse_height/test.py b/tests/regression_tests/pulse_height/test.py index c6459a8662b..90d960f6640 100644 --- a/tests/regression_tests/pulse_height/test.py +++ b/tests/regression_tests/pulse_height/test.py @@ -30,9 +30,11 @@ def sphere_model(): model.settings.batches = 5 model.settings.particles = 100 model.settings.photon_transport = True - model.settings.source = openmc.Source(space=openmc.stats.Point(), - energy=openmc.stats.Discrete([1e6],[1]), - particle='photon') + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point(), + energy=openmc.stats.Discrete([1e6], [1]), + particle='photon' + ) # Define tallies tally = openmc.Tally(name="pht tally") diff --git a/tests/regression_tests/time_cutoff/__init__.py b/tests/regression_tests/time_cutoff/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/time_cutoff/inputs_true.dat b/tests/regression_tests/time_cutoff/inputs_true.dat new file mode 100644 index 00000000000..7d3e94f50f9 --- /dev/null +++ b/tests/regression_tests/time_cutoff/inputs_true.dat @@ -0,0 +1,34 @@ + + + + + + + + + + fixed source + 100 + 10 + + + 0.0 0.0 0.0 + + + 10000.0 1.0 + + + + 1e-07 + + + + + 0.0 1e-07 2e-07 + + + 1 + flux + + + diff --git a/tests/regression_tests/time_cutoff/results_true.dat b/tests/regression_tests/time_cutoff/results_true.dat new file mode 100644 index 00000000000..1cafceb772e --- /dev/null +++ b/tests/regression_tests/time_cutoff/results_true.dat @@ -0,0 +1,5 @@ +tally 1: +2.000000E+03 +4.000000E+05 +0.000000E+00 +0.000000E+00 diff --git a/tests/regression_tests/time_cutoff/test.py b/tests/regression_tests/time_cutoff/test.py new file mode 100755 index 00000000000..8e13d57b73e --- /dev/null +++ b/tests/regression_tests/time_cutoff/test.py @@ -0,0 +1,42 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def time_model(): + model = openmc.Model() + time_cutoff = 1e-7 + + # A single sphere + s1 = openmc.Sphere(r=200, boundary_type='vacuum') + sphere = openmc.Cell() + sphere.region = -s1 + model.geometry = openmc.Geometry([sphere]) + + # Set the running parameters + settings_file = openmc.Settings() + settings_file.run_mode = 'fixed source' + settings_file.batches = 10 + settings_file.particles = 100 + settings_file.cutoff = {'time_neutron': time_cutoff} + settings_file.source = openmc.source.Source(space=openmc.stats.Point(), + energy=openmc.stats.Discrete([1e4], [1])) + model.settings = settings_file + + # Tally flux under time cutoff + tallies = openmc.Tallies() + tally = openmc.Tally() + tally.scores = ['flux'] + time_filter = openmc.TimeFilter([0, time_cutoff, 2*time_cutoff]) + tally.filters = [time_filter] + tallies.append(tally) + model.tallies = tallies + + return model + + +def test_time_cutoff(time_model): + harness = PyAPITestHarness('statepoint.10.h5', time_model) + harness.main() diff --git a/tests/unit_tests/mesh_to_vtk/test.py b/tests/unit_tests/mesh_to_vtk/test.py index 36928676cff..3ee6c0298ce 100644 --- a/tests/unit_tests/mesh_to_vtk/test.py +++ b/tests/unit_tests/mesh_to_vtk/test.py @@ -32,17 +32,19 @@ def diff_file(file1, file2): rect_mesh.y_grid = np.geomspace(5., 20., 10) rect_mesh.z_grid = np.linspace(1, 100, 20) -cyl_mesh = openmc.CylindricalMesh() -cyl_mesh.origin = (10, 10, -10) -cyl_mesh.r_grid = np.linspace(0, 5, 5) -cyl_mesh.phi_grid = np.linspace(0, 2 * np.pi, 4) -cyl_mesh.z_grid = np.linspace(0, 2, 4) - -sphere_mesh = openmc.SphericalMesh() -sphere_mesh.origin = (10, 10, -10) -sphere_mesh.r_grid = np.linspace(0, 5, 3) -sphere_mesh.theta_grid = np.linspace(0, 0.5 * np.pi, 4) -sphere_mesh.phi_grid = np.linspace(0, 2*np.pi, 8) +cyl_mesh = openmc.CylindricalMesh( + origin=(10, 10, -10), + r_grid=np.linspace(0, 5, 5), + phi_grid=np.linspace(0, 2 * np.pi, 4), + z_grid=np.linspace(0, 2, 4), +) + +sphere_mesh = openmc.SphericalMesh( + origin=(10, 10, -10), + r_grid=np.linspace(0, 5, 3), + theta_grid=np.linspace(0, 0.5 * np.pi, 4), + phi_grid=np.linspace(0, 2*np.pi, 8), +) def mesh_data(mesh_dims): diff --git a/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py b/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py index 4f23ef8591f..baa54498438 100644 --- a/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py +++ b/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py @@ -49,15 +49,17 @@ def model(): np.concatenate((-rectilinear_mesh.y_grid[::-1], rectilinear_mesh.y_grid)) rectilinear_mesh.z_grid = np.linspace(-10, 10, 11) -cylinder_mesh = openmc.CylindricalMesh() -cylinder_mesh.r_grid = np.linspace(0, 10, 23) +cylinder_mesh = openmc.CylindricalMesh( + r_grid=np.linspace(0, 10, 23), + z_grid=np.linspace(0, 1, 15) +) cylinder_mesh.phi_grid = np.linspace(0, np.pi, 21) -cylinder_mesh.z_grid = np.linspace(0, 1, 15) -spherical_mesh = openmc.SphericalMesh() -spherical_mesh.r_grid = np.linspace(1, 10, 30) -spherical_mesh.phi_grid = np.linspace(0, 0.8*np.pi, 25) -spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, 15) +spherical_mesh = openmc.SphericalMesh( + r_grid=np.linspace(1, 10, 30), + phi_grid=np.linspace(0, 0.8*np.pi, 25), + theta_grid=np.linspace(0, np.pi / 2, 15), +) MESHES = [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh] @@ -143,16 +145,17 @@ def test_write_data_to_vtk_size_mismatch(mesh): mesh.write_data_to_vtk(filename="out.vtk", datasets={"label": data}) def test_write_data_to_vtk_round_trip(run_in_tmpdir): - cmesh = openmc.CylindricalMesh() - cmesh.r_grid = (0.0, 1.0, 2.0) - cmesh.z_grid = (0.0, 2.0, 4.0, 5.0) - cmesh.phi_grid = (0.0, 3.0, 6.0) - - smesh = openmc.SphericalMesh() - smesh.r_grid = (0.0, 1.0, 2.0) - smesh.theta_grid = (0.0, 2.0, 4.0, 5.0) - smesh.phi_grid = (0.0, 3.0, 6.0) + cmesh = openmc.CylindricalMesh( + r_grid=(0.0, 1.0, 2.0), + z_grid=(0.0, 2.0, 4.0, 5.0), + phi_grid=(0.0, 3.0, 6.0), + ) + smesh = openmc.SphericalMesh( + r_grid=(0.0, 1.0, 2.0), + theta_grid=(0.0, 2.0, 4.0, 5.0), + phi_grid=(0.0, 3.0, 6.0), + ) rmesh = openmc.RegularMesh() rmesh.lower_left = (0.0, 0.0, 0.0) rmesh.upper_right = (1.0, 3.0, 5.0) @@ -276,11 +279,11 @@ def in_geom(cell): def test_sphere_mesh_coordinates(run_in_tmpdir): - mesh = openmc.SphericalMesh() - mesh.r_grid = np.linspace(0.1, 10, 30) - mesh.phi_grid = np.linspace(0, 1.5*np.pi, 25) - mesh.theta_grid = np.linspace(0, np.pi / 2, 15) - + mesh = openmc.SphericalMesh( + r_grid=np.linspace(0.1, 10, 30), + phi_grid=np.linspace(0, 1.5*np.pi, 25), + theta_grid=np.linspace(0, np.pi / 2, 15), + ) # write the data to a VTK file (no data) vtk_filename = 'test.vtk' mesh.write_data_to_vtk(vtk_filename, {}) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index 22b6ed8169e..df0c5d2bb53 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -33,11 +33,11 @@ def model(): settings.run_mode = 'fixed source' # build - mesh = openmc.CylindricalMesh() - mesh.phi_grid = np.linspace(0, 2*np.pi, 21) - mesh.z_grid = np.linspace(-geom_size, geom_size, 11) - mesh.r_grid = np.linspace(0, geom_size, geom_size) - + mesh = openmc.CylindricalMesh( + phi_grid=np.linspace(0, 2*np.pi, 21), + z_grid=np.linspace(-geom_size, geom_size, 11), + r_grid=np.linspace(0, geom_size, geom_size) + ) tally = openmc.Tally() mesh_filter = openmc.MeshFilter(mesh) @@ -127,10 +127,11 @@ def void_coincident_geom_model(): settings.particles = 1000 model.settings = settings - mesh = openmc.CylindricalMesh() - mesh.r_grid = np.linspace(0, 250, 501) - mesh.z_grid = [-250, 250] - mesh.phi_grid = np.linspace(0, 2*np.pi, 2) + mesh = openmc.CylindricalMesh( + r_grid=np.linspace(0, 250, 501), + z_grid=[-250, 250], + phi_grid=np.linspace(0, 2*np.pi, 2), + ) mesh_filter = openmc.MeshFilter(mesh) tally = openmc.Tally() diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index 457a6fea294..190affea99d 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -15,6 +15,27 @@ def res(): / 'test_reference.h5') return openmc.deplete.Results(filename) +def test_get_activity(res): + """Tests evaluating activity""" + t, a = res.get_activity("1") + + t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0]) + a_ref = np.array( + [1.25167956e+06, 3.71938527e+11, 4.43264300e+11, 3.55547176e+11]) + + np.testing.assert_allclose(t, t_ref) + np.testing.assert_allclose(a, a_ref) + + # Check by_nuclide + a_xe135_ref = np.array( + [2.106574218e+05, 1.227519888e+11, 1.177491828e+11, 1.031986176e+11]) + t_nuc, a_nuc = res.get_activity("1", by_nuclide=True) + + a_xe135 = np.array([a_nuc_i["Xe135"] for a_nuc_i in a_nuc]) + + np.testing.assert_allclose(t_nuc, t_ref) + np.testing.assert_allclose(a_xe135, a_xe135_ref) + def test_get_atoms(res): """Tests evaluating single nuclide concentration.""" @@ -43,6 +64,31 @@ def test_get_atoms(res): assert t_hour == pytest.approx(t_ref / (60 * 60)) +def test_get_decay_heat(res): + """Tests evaluating decay heat.""" + # Set chain file for testing + openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml' + + t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0]) + dh_ref = np.array( + [1.27933813e-09, 5.85347232e-03, 7.38773010e-03, 5.79954067e-03]) + + t, dh = res.get_decay_heat("1") + + np.testing.assert_allclose(t, t_ref) + np.testing.assert_allclose(dh, dh_ref) + + # Check by nuclide + dh_xe135_ref = np.array( + [1.27933813e-09, 7.45481920e-04, 7.15099509e-04, 6.26732849e-04]) + t_nuc, dh_nuc = res.get_decay_heat("1", by_nuclide=True) + + dh_nuc_xe135 = np.array([dh_nuc_i["Xe135"] for dh_nuc_i in dh_nuc]) + + np.testing.assert_allclose(t_nuc, t_ref) + np.testing.assert_allclose(dh_nuc_xe135, dh_xe135_ref) + + def test_get_mass(res): """Tests evaluating single nuclide concentration.""" t, n = res.get_mass("1", "Xe135") diff --git a/tests/unit_tests/test_filter_mesh.py b/tests/unit_tests/test_filter_mesh.py index c93bfa03770..36beed25d1c 100644 --- a/tests/unit_tests/test_filter_mesh.py +++ b/tests/unit_tests/test_filter_mesh.py @@ -22,15 +22,17 @@ def test_spherical_mesh_estimators(run_in_tmpdir): model.settings.inactive = 10 model.settings.batches = 20 - sph_mesh = openmc.SphericalMesh() - sph_mesh.r_grid = np.linspace(0.0, 5.0**3, 20)**(1/3) + sph_mesh = openmc.SphericalMesh( + r_grid=np.linspace(0.0, 5.0**3, 20)**(1/3) + ) tally1 = openmc.Tally() tally1.filters = [openmc.MeshFilter(sph_mesh)] tally1.scores = ['flux'] tally1.estimator = 'collision' - sph_mesh = openmc.SphericalMesh() - sph_mesh.r_grid = np.linspace(0.0, 5.0**3, 20)**(1/3) + sph_mesh = openmc.SphericalMesh( + r_grid=np.linspace(0.0, 5.0**3, 20)**(1/3) + ) tally2 = openmc.Tally() tally2.filters = [openmc.MeshFilter(sph_mesh)] tally2.scores = ['flux'] @@ -75,17 +77,19 @@ def test_cylindrical_mesh_estimators(run_in_tmpdir): model.settings.inactive = 10 model.settings.batches = 20 - cyl_mesh = openmc.CylindricalMesh() - cyl_mesh.r_grid = np.linspace(0.0, 5.0**3, 20)**(1/3) - cyl_mesh.z_grid = [-5., 5.] + cyl_mesh = openmc.CylindricalMesh( + r_grid=np.linspace(0.0, 5.0**3, 20)**(1/3), + z_grid=[-5., 5.] + ) tally1 = openmc.Tally() tally1.filters = [openmc.MeshFilter(cyl_mesh)] tally1.scores = ['flux'] tally1.estimator = 'collision' - cyl_mesh = openmc.CylindricalMesh() - cyl_mesh.r_grid = np.linspace(0.0, 5.0**3, 20)**(1/3) - cyl_mesh.z_grid = [-5., 5.] + cyl_mesh = openmc.CylindricalMesh( + r_grid=np.linspace(0.0, 5.0**3, 20)**(1/3), + z_grid=[-5., 5.] + ) tally2 = openmc.Tally() tally2.filters = [openmc.MeshFilter(cyl_mesh)] tally2.scores = ['flux'] @@ -133,10 +137,11 @@ def test_cylindrical_mesh_coincident(scale, run_in_tmpdir): model.settings.batches = 10 model.settings.inactive = 0 - cyl_mesh = openmc.CylindricalMesh() - cyl_mesh.r_grid = [0., 1.25*scale] - cyl_mesh.phi_grid = [0., 2*math.pi] - cyl_mesh.z_grid = [-1e10, 1e10] + cyl_mesh = openmc.CylindricalMesh( + r_grid=[0., 1.25*scale], + phi_grid=[0., 2*math.pi], + z_grid=[-1e10, 1e10] + ) cyl_mesh_filter = openmc.MeshFilter(cyl_mesh) cell_filter = openmc.CellFilter([cell1]) @@ -183,10 +188,12 @@ def test_spherical_mesh_coincident(scale, run_in_tmpdir): model.settings.batches = 10 model.settings.inactive = 0 - sph_mesh = openmc.SphericalMesh() - sph_mesh.r_grid = [0., 1.25*scale] - sph_mesh.phi_grid = [0., 2*math.pi] - sph_mesh.theta_grid = [0., math.pi] + sph_mesh = openmc.SphericalMesh( + r_grid=[0., 1.25*scale], + phi_grid=[0., 2*math.pi], + theta_grid=[0., math.pi], + ) + sph_mesh_filter = openmc.MeshFilter(sph_mesh) cell_filter = openmc.CellFilter([cell1]) @@ -227,10 +234,11 @@ def test_get_reshaped_data(run_in_tmpdir): model.settings.inactive = 10 model.settings.batches = 20 - sph_mesh = openmc.SphericalMesh() - sph_mesh.r_grid = np.linspace(0.0, 5.0**3, 20)**(1/3) - sph_mesh.theta_grid = np.linspace(0, math.pi, 4) - sph_mesh.phi_grid = np.linspace(0, 2*math.pi, 3) + sph_mesh = openmc.SphericalMesh( + r_grid=np.linspace(0.0, 5.0**3, 20)**(1/3), + theta_grid=np.linspace(0, math.pi, 4), + phi_grid=np.linspace(0, 2*math.pi, 3) + ) tally1 = openmc.Tally() efilter = openmc.EnergyFilter([0, 1e5, 1e8]) meshfilter = openmc.MeshFilter(sph_mesh) diff --git a/tests/unit_tests/test_geometry.py b/tests/unit_tests/test_geometry.py index 6fa34d26d01..ebb185642fd 100644 --- a/tests/unit_tests/test_geometry.py +++ b/tests/unit_tests/test_geometry.py @@ -362,8 +362,13 @@ def get_cyl_cell(r1, r2, z1, z2, fill): clad = get_cyl_cell(r1, r2, z1, z2, m2) water = get_cyl_cell(r2, r3, z1, z2, m3) root = openmc.Universe(cells=[fuel, clad, water]) - geom = openmc.Geometry(root) - geom.merge_surfaces=True + geom = openmc.Geometry(root=root, merge_surfaces=True, surface_precision=11) + assert geom.merge_surfaces is True + geom.merge_surfaces = False + assert geom.merge_surfaces is False + assert geom.surface_precision == 11 + geom.surface_precision = 10 + assert geom.surface_precision == 10 model = openmc.model.Model(geometry=geom, materials=openmc.Materials([m1, m2, m3])) diff --git a/tests/unit_tests/test_lattice_discretization.py b/tests/unit_tests/test_lattice_discretization.py index b21d7e13e0e..c7779ea8bf5 100644 --- a/tests/unit_tests/test_lattice_discretization.py +++ b/tests/unit_tests/test_lattice_discretization.py @@ -34,8 +34,8 @@ def test_discretization_clone_only_some_materials(rlat2): fuel1 = next(iter(rlat_clone.get_universe((0, 0)).cells.values())).fill rlat_clone.discretize(materials_to_clone=[fuel1]) - assert next(reversed(rlat_clone.get_universe((0, 0)).cells.values())).fill\ - == next(reversed(rlat_clone.get_universe((1, 0)).cells.values())).fill + assert next(reversed(list(rlat_clone.get_universe((0, 0)).cells.values()))).fill\ + == next(reversed(list(rlat_clone.get_universe((1, 0)).cells.values()))).fill assert next(iter(rlat_clone.get_universe((0, 0)).cells.values())).fill\ != next(iter(rlat_clone.get_universe((1, 0)).cells.values())).fill diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 4a09df0455d..480104154d5 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -1,6 +1,9 @@ -import openmc -import pytest +from math import pi + import numpy as np +import pytest +import openmc + @pytest.mark.parametrize("val_left,val_right", [(0, 0), (-1., -1.), (2.0, 2)]) def test_raises_error_when_flat(val_left, val_right): @@ -35,12 +38,107 @@ def test_raises_error_when_flat(val_left, val_right): mesh.lower_left = [-25, -25, val_left] -def test_mesh_bounding_box(): +def test_regular_mesh_bounding_box(): mesh = openmc.RegularMesh() - mesh.lower_left = [-2, -3 ,-5] + mesh.lower_left = [-2, -3, -5] mesh.upper_right = [2, 3, 5] bb = mesh.bounding_box assert isinstance(bb, openmc.BoundingBox) - np.testing.assert_array_equal(bb.lower_left, np.array([-2, -3 ,-5])) - np.testing.assert_array_equal(bb.upper_right, np.array([2, 3, 5])) + np.testing.assert_array_equal(bb.lower_left, (-2, -3 ,-5)) + np.testing.assert_array_equal(bb.upper_right, (2, 3, 5)) + + +def test_cylindrical_mesh_bounding_box(): + # test with mesh at origin (0, 0, 0) + mesh = openmc.CylindricalMesh( + r_grid=[0.1, 0.2, 0.5, 1.], + z_grid=[0.1, 0.2, 0.4, 0.6, 1.], + origin=(0, 0, 0) + ) + np.testing.assert_array_equal(mesh.upper_right, (1, 1, 1)) + np.testing.assert_array_equal(mesh.lower_left, (-1, -1, -1)) + bb = mesh.bounding_box + assert isinstance(bb, openmc.BoundingBox) + np.testing.assert_array_equal(bb.lower_left, (-1, -1, -1)) + np.testing.assert_array_equal(bb.upper_right, (1, 1, 1)) + + # test with mesh at origin (3, 5, 7) + mesh.origin = (3, 5, 7) + np.testing.assert_array_equal(mesh.upper_right, (4, 6, 8)) + np.testing.assert_array_equal(mesh.lower_left, (2, 4, 6)) + bb = mesh.bounding_box + assert isinstance(bb, openmc.BoundingBox) + np.testing.assert_array_equal(bb.lower_left, (2, 4, 6)) + np.testing.assert_array_equal(bb.upper_right, (4, 6, 8)) + + +def test_spherical_mesh_bounding_box(): + # test with mesh at origin (0, 0, 0) + mesh = openmc.SphericalMesh([0.1, 0.2, 0.5, 1.], origin=(0., 0., 0.)) + np.testing.assert_array_equal(mesh.upper_right, (1, 1, 1)) + np.testing.assert_array_equal(mesh.lower_left, (-1, -1, -1)) + bb = mesh.bounding_box + assert isinstance(bb, openmc.BoundingBox) + np.testing.assert_array_equal(bb.lower_left, (-1, -1, -1)) + np.testing.assert_array_equal(bb.upper_right, (1, 1, 1)) + + # test with mesh at origin (3, 5, 7) + mesh.origin = (3, 5, 7) + np.testing.assert_array_equal(mesh.upper_right, (4, 6, 8)) + np.testing.assert_array_equal(mesh.lower_left, (2, 4, 6)) + bb = mesh.bounding_box + assert isinstance(bb, openmc.BoundingBox) + np.testing.assert_array_equal(bb.lower_left, (2, 4, 6)) + np.testing.assert_array_equal(bb.upper_right, (4, 6, 8)) + + +def test_SphericalMesh_initiation(): + # test defaults + mesh = openmc.SphericalMesh(r_grid=(0, 10)) + assert (mesh.origin == np.array([0, 0, 0])).all() + assert mesh.r_grid == (0, 10) + assert (mesh.theta_grid == np.array([0, pi])).all() + assert (mesh.phi_grid == np.array([0, 2*pi])).all() + + # test setting on creation + mesh = openmc.SphericalMesh( + origin=(1, 2, 3), + r_grid=(0, 2), + theta_grid=(1, 3), + phi_grid=(2, 4) + ) + assert (mesh.origin == np.array([1, 2, 3])).all() + assert mesh.r_grid == (0., 2.) + assert mesh.theta_grid == (1, 3) + assert mesh.phi_grid == (2, 4) + + # test attribute changing + mesh.r_grid = (0, 11) + assert (mesh.r_grid == np.array([0., 11.])).all() + + +def test_CylindricalMesh_initiation(): + # test defaults + mesh = openmc.CylindricalMesh(r_grid=(0, 10), z_grid=(0, 10)) + assert (mesh.origin == np.array([0, 0, 0])).all() + assert mesh.r_grid == (0, 10) + assert mesh.phi_grid == (0, 2*pi) + assert mesh.z_grid == (0, 10) + + # test setting on creation + mesh = openmc.CylindricalMesh( + origin=(1, 2, 3), + r_grid=(0, 2), + z_grid=(1, 3), + phi_grid=(2, 4) + ) + assert (mesh.origin == np.array([1, 2, 3])).all() + assert (mesh.r_grid == np.array([0., 2.])).all() + assert mesh.z_grid == (1, 3) + assert (mesh.phi_grid == np.array([2, 4])).all() + # test attribute changing + mesh.r_grid = (0., 10.) + assert (mesh.r_grid == np.array([0, 10.])).all() + mesh.z_grid = (0., 4.) + assert (mesh.z_grid == np.array([0, 4.])).all() diff --git a/tests/unit_tests/test_mesh_from_domain.py b/tests/unit_tests/test_mesh_from_domain.py index b4edae196f8..015c8f88f0d 100644 --- a/tests/unit_tests/test_mesh_from_domain.py +++ b/tests/unit_tests/test_mesh_from_domain.py @@ -9,7 +9,7 @@ def test_reg_mesh_from_cell(): surface = openmc.Sphere(r=10, x0=2, y0=3, z0=5) cell = openmc.Cell(region=-surface) - mesh = openmc.RegularMesh.from_domain(cell, dimension=[7, 11, 13]) + mesh = openmc.RegularMesh.from_domain(domain=cell, dimension=[7, 11, 13]) assert isinstance(mesh, openmc.RegularMesh) assert np.array_equal(mesh.dimension, (7, 11, 13)) assert np.array_equal(mesh.lower_left, cell.bounding_box[0]) @@ -23,7 +23,7 @@ def test_cylindrical_mesh_from_cell(): z_surface_1 = openmc.ZPlane(z0=30) z_surface_2 = openmc.ZPlane(z0=0) cell = openmc.Cell(region=-cy_surface & -z_surface_1 & +z_surface_2) - mesh = openmc.CylindricalMesh.from_domain(cell, dimension=[2, 4, 3]) + mesh = openmc.CylindricalMesh.from_domain(domain=cell, dimension=[2, 4, 3]) assert isinstance(mesh, openmc.CylindricalMesh) assert np.array_equal(mesh.dimension, (2, 4, 3)) @@ -38,7 +38,7 @@ def test_reg_mesh_from_region(): surface = openmc.Sphere(r=1, x0=-5, y0=-3, z0=-2) region = -surface - mesh = openmc.RegularMesh.from_domain(region) + mesh = openmc.RegularMesh.from_domain(domain=region) assert isinstance(mesh, openmc.RegularMesh) assert np.array_equal(mesh.dimension, (10, 10, 10)) # default values assert np.array_equal(mesh.lower_left, region.bounding_box[0]) @@ -53,7 +53,7 @@ def test_cylindrical_mesh_from_region(): z_surface_2 = openmc.ZPlane(z0=-30) cell = openmc.Cell(region=-cy_surface & -z_surface_1 & +z_surface_2) mesh = openmc.CylindricalMesh.from_domain( - cell, + domain=cell, dimension=(6, 2, 3), phi_grid_bounds=(0., np.pi) ) diff --git a/tests/unit_tests/test_plots.py b/tests/unit_tests/test_plots.py index 5e5e1f22c12..922e4d9dc95 100644 --- a/tests/unit_tests/test_plots.py +++ b/tests/unit_tests/test_plots.py @@ -199,3 +199,25 @@ def test_plots(run_in_tmpdir): assert new_plots[0].colors == p1.colors assert new_plots[0].mask_components == p1.mask_components assert new_plots[1].origin == p2.origin + + +def test_voxel_plot_roundtrip(): + # Define a voxel plot and create XML element + plot = openmc.Plot(name='my voxel plot') + plot.type = 'voxel' + plot.filename = 'voxel1' + plot.pixels = (50, 50, 50) + plot.origin = (0., 0., 0.) + plot.width = (75., 75., 75.) + plot.color_by = 'material' + elem = plot.to_xml_element() + + # Read back from XML and make sure it hasn't changed + new_plot = plot.from_xml_element(elem) + assert new_plot.name == plot.name + assert new_plot.filename == plot.filename + assert new_plot.type == plot.type + assert new_plot.pixels == plot.pixels + assert new_plot.origin == plot.origin + assert new_plot.width == plot.width + assert new_plot.color_by == plot.color_by diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index d75a50477a3..4b2cc18240f 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -23,10 +23,13 @@ def test_export_to_xml(run_in_tmpdir): s.surf_source_write = {'surface_ids': [2], 'max_particles': 200} s.confidence_intervals = True s.ptables = True + s.plot_seed = 100 s.survival_biasing = True s.cutoff = {'weight': 0.25, 'weight_avg': 0.5, 'energy_neutron': 1.0e-5, 'energy_photon': 1000.0, 'energy_electron': 1.0e-5, - 'energy_positron': 1.0e-5} + 'energy_positron': 1.0e-5, 'time_neutron': 1.0e-5, + 'time_photon': 1.0e-5, 'time_electron': 1.0e-5, + 'time_positron': 1.0e-5} mesh = openmc.RegularMesh() mesh.lower_left = (-10., -10., -10.) mesh.upper_right = (10., 10., 10.) @@ -82,11 +85,14 @@ def test_export_to_xml(run_in_tmpdir): assert s.surf_source_write == {'surface_ids': [2], 'max_particles': 200} assert s.confidence_intervals assert s.ptables + assert s.plot_seed == 100 assert s.seed == 17 assert s.survival_biasing assert s.cutoff == {'weight': 0.25, 'weight_avg': 0.5, 'energy_neutron': 1.0e-5, 'energy_photon': 1000.0, - 'energy_electron': 1.0e-5, 'energy_positron': 1.0e-5} + 'energy_electron': 1.0e-5, 'energy_positron': 1.0e-5, + 'time_neutron': 1.0e-5, 'time_photon': 1.0e-5, + 'time_electron': 1.0e-5, 'time_positron': 1.0e-5} assert isinstance(s.entropy_mesh, openmc.RegularMesh) assert s.entropy_mesh.lower_left == [-10., -10., -10.] assert s.entropy_mesh.upper_right == [10., 10., 10.] @@ -108,7 +114,7 @@ def test_export_to_xml(run_in_tmpdir): 'energy_min': 1.0, 'energy_max': 1000.0, 'nuclides': ['U235', 'U238', 'Pu239']} assert s.create_fission_neutrons - assert not s.create_delayed_neutrons + assert not s.create_delayed_neutrons assert s.log_grid_bins == 2000 assert not s.photon_transport assert s.electron_treatment == 'led' diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py index 80f65bb420e..b491014cc70 100644 --- a/tests/unit_tests/test_spherical_mesh.py +++ b/tests/unit_tests/test_spherical_mesh.py @@ -33,11 +33,11 @@ def model(): settings.run_mode = 'fixed source' # build - mesh = openmc.SphericalMesh() - mesh.phi_grid = np.linspace(0, 2*np.pi, 13) - mesh.theta_grid = np.linspace(0, np.pi, 7) - mesh.r_grid = np.linspace(0, geom_size, geom_size) - + mesh = openmc.SphericalMesh( + phi_grid=np.linspace(0, 2*np.pi, 13), + theta_grid=np.linspace(0, np.pi, 7), + r_grid=np.linspace(0, geom_size, geom_size), + ) tally = openmc.Tally() mesh_filter = openmc.MeshFilter(mesh) @@ -135,8 +135,7 @@ def void_coincident_geom_model(): settings.particles = 5000 model.settings = settings - mesh = openmc.SphericalMesh() - mesh.r_grid = np.linspace(0, 250, 501) + mesh = openmc.SphericalMesh(r_grid=np.linspace(0, 250, 501)) mesh_filter = openmc.MeshFilter(mesh) tally = openmc.Tally() diff --git a/tests/unit_tests/test_universe.py b/tests/unit_tests/test_universe.py index f4cf43d919b..630e66df3ab 100644 --- a/tests/unit_tests/test_universe.py +++ b/tests/unit_tests/test_universe.py @@ -61,14 +61,17 @@ def test_plot(run_in_tmpdir, sphere_model): } for basis in ('xy', 'yz', 'xz'): - pincell.geometry.root_universe.plot( + plot = pincell.geometry.root_universe.plot( colors=mat_colors, color_by="material", legend=True, pixels=(10, 10), basis=basis, - outline=True + outline=True, + axis_units='m' ) + assert plot.xaxis.get_label().get_text() == f'{basis[0]} [m]' + assert plot.yaxis.get_label().get_text() == f'{basis[1]} [m]' # model with no inf values in bounding box m = sphere_model.materials[0] @@ -77,7 +80,7 @@ def test_plot(run_in_tmpdir, sphere_model): colors = {m: 'limegreen'} for basis in ('xy', 'yz', 'xz'): - univ.plot( + plot = univ.plot( colors=colors, color_by="cell", legend=False, @@ -85,6 +88,17 @@ def test_plot(run_in_tmpdir, sphere_model): basis=basis, outline=False ) + assert plot.xaxis.get_label().get_text() == f'{basis[0]} [cm]' + assert plot.yaxis.get_label().get_text() == f'{basis[1]} [cm]' + + msg = "Must pass 'colors' dictionary if you are adding a legend via legend=True." + # This plot call should fail as legend is True but colors is None + with pytest.raises(ValueError, match=msg): + univ.plot( + color_by="cell", + legend=True, + pixels=100, + ) def test_get_nuclides(uo2): diff --git a/tests/unit_tests/test_volume.py b/tests/unit_tests/test_volume.py new file mode 100644 index 00000000000..f49bb35018c --- /dev/null +++ b/tests/unit_tests/test_volume.py @@ -0,0 +1,15 @@ +import numpy as np +import pytest + +import openmc + + +def test_infinity_handling(): + surf1 = openmc.Sphere(boundary_type="vacuum") + cell1 = openmc.Cell(region=-surf1) + + lower_left = (-2, -np.inf, -2) + upper_right = (np.inf, 2, 2) + + with pytest.raises(ValueError, match="must be finite"): + openmc.VolumeCalculation([cell1], 100, lower_left, upper_right) diff --git a/tests/unit_tests/weightwindows/test_ww_gen.py b/tests/unit_tests/weightwindows/test_ww_gen.py index a1aaca5b0eb..072e332cb75 100644 --- a/tests/unit_tests/weightwindows/test_ww_gen.py +++ b/tests/unit_tests/weightwindows/test_ww_gen.py @@ -206,19 +206,37 @@ def test_ww_import_export(run_in_tmpdir, model): lb_before, up_before = ww.bounds + # set some additional weight windows properties after transport + ww.survival_ratio = 0.7 + assert ww.survival_ratio == 0.7 + + ww.weight_cutoff = 1e-10 + assert ww.weight_cutoff == 1e-10 + + ww.max_lower_bound_ratio = 200.0 + assert ww.max_lower_bound_ratio == 200.0 + + ww.max_split = 26000 + assert ww.max_split == 26000 + openmc.lib.export_weight_windows() assert Path('weight_windows.h5').exists() openmc.lib.import_weight_windows('weight_windows.h5') - ww = openmc.lib.weight_windows[2] + imported_ww = openmc.lib.weight_windows[2] - lb_after, up_after = ww.bounds + lb_after, up_after = imported_ww.bounds assert np.allclose(lb_before, lb_after) assert np.allclose(up_before, up_after) + assert ww.survival_ratio == imported_ww.survival_ratio + assert ww.max_lower_bound_ratio == imported_ww.max_lower_bound_ratio + assert ww.weight_cutoff == imported_ww.weight_cutoff + assert ww.max_split == imported_ww.max_split + openmc.lib.finalize() diff --git a/tools/ci/gha-install.sh b/tools/ci/gha-install.sh index fb81fcea2d1..2cef974d346 100755 --- a/tools/ci/gha-install.sh +++ b/tools/ci/gha-install.sh @@ -40,7 +40,7 @@ if [[ $MPI == 'y' ]]; then export CC=mpicc export HDF5_MPI=ON export HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/mpich - pip install wheel cython + pip install wheel "cython<3.0" pip install --no-binary=h5py --no-build-isolation h5py fi