diff --git a/Dockerfile b/Dockerfile index ae59ec08943..fa223bd2284 100644 --- a/Dockerfile +++ b/Dockerfile @@ -91,7 +91,7 @@ RUN cd $HOME \ RUN if [ "$build_dagmc" = "on" ]; then \ # Install addition packages required for DAGMC apt-get -y install libeigen3-dev libnetcdf-dev libtbb-dev libglfw3-dev \ - && pip install --upgrade numpy cython \ + && pip install --upgrade numpy "cython<3.0" \ # Clone and install EMBREE && mkdir -p $HOME/EMBREE && cd $HOME/EMBREE \ && git clone --single-branch -b ${EMBREE_TAG} --depth 1 ${EMBREE_REPO} \ 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 a52de1b9048..db90d93dc0d 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -508,14 +508,15 @@ attributes/sub-elements: as a unique pointer. This function can be used to pass parameters through to the source from the XML, if needed. - More documentation on how to build sources can be found in :ref:`custom_source`. + More documentation on how to build sources can be found in + :ref:`compiled_source`. :parameters: If this attribute is given, it indicated that the source type is ``compiled``. Its value provides the parameters to pass through to the class generated using the ``library`` parameter. More documentation on how to build parametrized sources can be found in - :ref:`parameterized_custom_source`. + :ref:`parameterized_compiled_source`. :space: An element specifying the spatial distribution of source sites. This element diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index ffd521a7583..afecaff5333 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -21,7 +21,10 @@ Simulation Settings :nosignatures: :template: myclass.rst - openmc.Source + openmc.SourceBase + openmc.IndependentSource + openmc.FileSource + openmc.CompiledSource openmc.SourceParticle openmc.VolumeCalculation openmc.Settings @@ -42,9 +45,6 @@ Material Specification :nosignatures: :template: myclass.rst - openmc.Nuclide - openmc.Element - openmc.Macroscopic openmc.Material openmc.Materials diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index dc9370ab595..8731a9a1383 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -57,6 +57,16 @@ The :class:`CoupledOperator` and :class:`IndependentOperator` classes must also have some knowledge of how nuclides transmute and decay. This is handled by the :class:`Chain` class. +The :class:`IndependentOperator` class requires a set of fluxes and microscopic +cross sections. The following function can be used to generate this information: + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myfunction.rst + + get_microxs_and_flux + Minimal Example --------------- diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 439e3030528..261900ce61a 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -197,43 +197,54 @@ across all material instances. Transport-independent depletion =============================== -.. warning:: - - This feature is still under heavy development and has yet to be rigorously - verified. API changes and feature additions are possible and likely in - the near future. - -This category of operator uses one-group microscopic cross sections to obtain -transmutation reaction rates. The cross sections are pre-calculated, so there is -no need for direct coupling between a transport-independent operator and a -transport solver. The :mod:`openmc.deplete` module offers a single -transport-independent operator, :class:`~openmc.deplete.IndependentOperator`, -and only one operator is needed since, in theory, any transport code could -calcuate the one-group microscopic cross sections. +This category of operator uses multigroup microscopic cross sections along with +multigroup flux spectra to obtain transmutation reaction rates. The cross +sections are pre-calculated, so there is no need for direct coupling between a +transport-independent operator and a transport solver. The :mod:`openmc.deplete` +module offers a single transport-independent operator, +:class:`~openmc.deplete.IndependentOperator`, and only one operator is needed +since, in theory, any transport code could calculate the multigroup microscopic +cross sections. The :class:`~openmc.deplete.IndependentOperator` class has two +constructors. The default constructor requires a :class:`openmc.Materials` +instance, a list of multigroup flux arrays, and a list of +:class:`~openmc.deplete.MicroXS` instances containing multigroup microscopic +cross sections in units of barns. This might look like the following:: + + materials = openmc.Materials([m1, m2, m3]) + ... -The :class:`~openmc.deplete.IndependentOperator` class has two constructors. -The default constructor requires a :class:`openmc.Materials` instance, a -:class:`~openmc.deplete.MicroXS` instance containing one-group microscoic cross -sections in units of barns, and a path to a depletion chain file:: + # Assign fluxes (generated from any code) + flux_m1 = numpy.array([...]) + flux_m2 = numpy.array([...]) + flux_m3 = numpy.array([...]) + fluxes = [flux_m1, flux_m2, flux_m3] - materials = openmc.Materials() - ... + # Assign microscopic cross sections + micro_m1 = openmc.deplete.MicroXS.from_csv('xs_m1.csv') + micro_m2 = openmc.deplete.MicroXS.from_csv('xs_m2.csv') + micro_m3 = openmc.deplete.MicroXS.from_csv('xs_m3.csv') + micros = [micro_m1, micro_m2, micro_m3] - # load in the microscopic cross sections - micro_xs = openmc.deplete.MicroXS.from_csv(micro_xs_path) + # Create operator + op = openmc.deplete.IndependentOperator(materials, fluxes, micros) - op = openmc.deplete.IndependentOperator(materials, micro_xs, chain_file) +For more details on the :class:`~openmc.deplete.MicroXS` class, including how to +use OpenMC's transport solver to generate microscopic cross sections and fluxes +for use with :class:`~openmc.deplete.IndependentOperator`, see :ref:`micros`. .. note:: - The same statements from :ref:`coupled-depletion` about which - materials are depleted and the requirement for depletable materials to have - a specified volume also apply here. + The same statements from :ref:`coupled-depletion` about which materials are + depleted and the requirement for depletable materials to have a specified + volume also apply here. An alternate constructor, :meth:`~openmc.deplete.IndependentOperator.from_nuclides`, accepts a volume and dictionary of nuclide concentrations in place of the :class:`openmc.Materials` -instance:: +instance. Note that while the normal constructor allows multiple materials to be +depleted with a single operator, the +:meth:`~openmc.deplete.IndependentOperator.from_nuclides` classmethod only works +for a single material:: nuclides = {'U234': 8.92e18, 'U235': 9.98e20, @@ -244,6 +255,7 @@ instance:: volume = 0.5 op = openmc.deplete.IndependentOperator.from_nuclides(volume, nuclides, + flux, micro_xs, chain_file, nuc_units='atom/cm3') @@ -253,18 +265,21 @@ transport-depletion calculation and follow the same steps from there. .. note:: - Ideally, one-group cross section data should be available for every - reaction in the depletion chain. If cross section data is not present for - a nuclide in the depletion chain with at least one reaction, that reaction - will not be simulated. + Ideally, multigroup cross section data should be available for every reaction + in the depletion chain. If cross section data is not present for a nuclide in + the depletion chain with at least one reaction, that reaction will not be + simulated. + +.. _micros: Loading and Generating Microscopic Cross Sections ------------------------------------------------- -As mentioned earlier, any transport code could be used to calculate one-group -microscopic cross sections. The :mod:`openmc.deplete` module provides the -:class:`~openmc.deplete.MicroXS` class, which contains methods to read in -pre-calculated cross sections from a ``.csv`` file or from data arrays:: +As mentioned above, any transport code could be used to calculate multigroup +microscopic cross sections and fluxes. The :mod:`openmc.deplete` module provides +the :class:`~openmc.deplete.MicroXS` class, which can either be instantiated +from pre-calculated cross sections in a ``.csv`` file or from data arrays +directly:: micro_xs = MicroXS.from_csv(micro_xs_path) @@ -273,37 +288,31 @@ pre-calculated cross sections from a ``.csv`` file or from data arrays:: data = np.array([[0.1, 0.2], [0.3, 0.4], [0.01, 0.5]]) - micro_xs = MicroXS.from_array(nuclides, reactions, data) + micro_xs = MicroXS(data, nuclides, reactions) .. important:: - Both :meth:`~openmc.deplete.MicroXS.from_csv()` and - :meth:`~openmc.deplete.MicroXS.from_array()` assume the cross section values - provided are in barns by defualt, but have no way of verifying this. Make - sure your cross sections are in the correct units before passing to a + The cross section values are assumed to be in units of barns. Make sure your + cross sections are in the correct units before passing to a :class:`~openmc.deplete.IndependentOperator` object. -The :class:`~openmc.deplete.MicroXS` class also contains a method to generate one-group microscopic cross sections using OpenMC's transport solver. The -:meth:`~openmc.deplete.MicroXS.from_model()` method will produce a -:class:`~openmc.deplete.MicroXS` instance with microscopic cross section data in -units of barns:: - - import openmc +Additionally, a convenience function, +:func:`~openmc.deplete.get_microxs_and_flux`, can provide the needed fluxes and +cross sections using OpenMC's transport solver:: - model = openmc.Model.from_xml() + model = openmc.Model() + ... - micro_xs = openmc.deplete.MicroXS.from_model(model, - model.materials[0], - chain_file) + fluxes, micros = openmc.deplete.get_microxs_and_flux(model, materials) -If you are running :meth:`~openmc.deplete.MicroXS.from_model()` on a cluster +If you are running :func:`~openmc.deplete.get_microxs_and_flux` on a cluster where temporary files are created on a local filesystem that is not shared across nodes, you'll need to set an environment variable pointing to a local directoy so that each MPI process knows where to store output files used to calculate the microscopic cross sections. In order of priority, they are -:envvar:`TMPDIR`. :envvar:`TEMP`, and :envvar:`TMP`. Users interested in -further details can read the documentation for the `tempfile `_ module. - +:envvar:`TMPDIR`. :envvar:`TEMP`, and :envvar:`TMP`. Users interested in further +details can read the documentation for the `tempfile +`_ module. Caveats ------- @@ -325,24 +334,30 @@ normalizing reaction rates: the time integrator is a flux, and obtains the reaction rates by multiplying the cross sections by the ``source-rate``. 2. ``fission-q`` normalization, which uses the ``power`` or ``power_density`` - provided by the time integrator to obtain reaction rates by computing a value - for the flux based on this power. The equation we use for this calculation is + provided by the time integrator to obtain normalized reaction rates by + computing a normalization factor as the ratio of the user-specified power to + the "observed" power based on fission reaction rates. The equation for the + normalization factor is .. math:: :label: fission-q - \phi = \frac{P}{\sum\limits_i (Q_i \sigma^f_i N_i)} + f = \frac{P}{\sum\limits_m \sum\limits_i \left(Q_i N_{i,m} \sum\limits_g + \sigma^f_{i,g,m} \phi_{g,m} \right)} where :math:`P` is the power, :math:`Q_i` is the fission Q value for nuclide - :math:`i`, :math:`\sigma_i^f` is the microscopic fission cross section for - nuclide :math:`i`, and :math:`N_i` is the number of atoms of nuclide - :math:`i`. This equation makes the same assumptions and issues as discussed - in :ref:`energy-deposition`. Unfortunately, the proposed solution in that - section does not apply here since we are decoupled from transport code. - However, there is a method to converge to a more accurate value for flux by - using substeps during time integration. `This paper - `_ provides a good discussion - of this method. + :math:`i`, :math:`\sigma_{i,g,m}^f` is the microscopic fission cross section + for nuclide :math:`i` in energy group :math:`g` for material :math:`m`, + :math:`\phi_{g,m}` is the neutron flux in group :math:`g` for material + :math:`m`, and :math:`N_{i,m}` is the number of atoms of nuclide :math:`i` + for material :math:`m`. Reaction rates are then multiplied by :math:`f` so + that the total fission power matches :math:`P`. This equation makes the same + assumptions and issues as discussed in :ref:`energy-deposition`. + Unfortunately, the proposed solution in that section does not apply here + since we are decoupled from transport code. However, there is a method to + converge to a more accurate value for flux by using substeps during time + integration. `This paper `_ + provides a good discussion of this method. .. warning:: @@ -359,14 +374,14 @@ useful for running many different cases of a particular scenario. A transport-independent depletion simulation using ``fission-q`` normalization will sum the fission energy values across all materials into :math:`Q_i` in Equation :math:numref:`fission-q`, and Equation :math:numref:`fission-q` -provides the flux we use to calculate the reaction rates in each material. +provides the normalization factor applied to reaction rates in each material. This can be useful for running a scenario with multiple depletable materials that are part of the same reactor. This behavior may change in the future. Time integration ~~~~~~~~~~~~~~~~ -The values of the one-group microscopic cross sections passed to +The values of the microscopic cross sections passed to :class:`openmc.deplete.IndependentOperator` are fixed for the entire depletion simulation. This implicit assumption may produce inaccurate results for certain scenarios. diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 69089d14d33..81fa78991a7 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -169,16 +169,19 @@ External Source Distributions External source distributions can be specified through the :attr:`Settings.source` attribute. If you have a single external source, you can -create an instance of :class:`openmc.Source` and use it to set the -:attr:`Settings.source` attribute. If you have multiple external sources with -varying source strengths, :attr:`Settings.source` should be set to a list of -:class:`openmc.Source` objects. - -The :class:`openmc.Source` class has four main attributes that one can set: -:attr:`Source.space`, which defines the spatial distribution, -:attr:`Source.angle`, which defines the angular distribution, -:attr:`Source.energy`, which defines the energy distribution, and -:attr:`Source.time`, which defines the time distribution. +create an instance of any of the subclasses of :class:`openmc.SourceBase` +(:class:`openmc.IndependentSource`, :class:`openmc.FileSource`, +:class:`openmc.CompiledSource`) and use it to set the :attr:`Settings.source` +attribute. If you have multiple external sources with varying source strengths, +:attr:`Settings.source` should be set to a list of :class:`openmc.SourceBase` +objects. + +The :class:`openmc.IndependentSource` class is the primary class for defining +source distributions and has four main attributes that one can set: +:attr:`IndependentSource.space`, which defines the spatial distribution, +:attr:`IndependentSource.angle`, which defines the angular distribution, +:attr:`IndependentSource.energy`, which defines the energy distribution, and +:attr:`IndependentSource.time`, which defines the time distribution. The spatial distribution can be set equal to a sub-class of :class:`openmc.stats.Spatial`; common choices are :class:`openmc.stats.Point` or @@ -232,10 +235,10 @@ would run:: source.time = openmc.stats.Uniform(0, 1e-6) settings.source = source -The :class:`openmc.Source` class also has a :attr:`Source.strength` attribute -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:: +All subclasses of :class:`openmc.SourceBase` have a :attr:`SourceBase.strength` +attribute 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.IndependentSource() src1.strength = 0.7 @@ -247,9 +250,9 @@ time and another that should be sampled 30% of the time:: settings.source = [src1, src2] -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:: +Finally, the :attr:`IndependentSource.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.IndependentSource() source.particle = 'photon' @@ -263,10 +266,10 @@ For a full list of all classes related to statistical distributions, see File-based Sources ------------------ -OpenMC can use a pregenerated HDF5 source file by specifying the ``filename`` -argument to :class:`openmc.Source`:: +OpenMC can use a pregenerated HDF5 source file through the +:class:`openmc.FileSource` class:: - settings.source = openmc.IndependentSource(filename='source.h5') + settings.source = openmc.FileSource('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 @@ -286,10 +289,10 @@ attribute:: In this example, at most 10,000 source particles are stored when particles cross surfaces with IDs of 1, 2, or 3. -.. _custom_source: +.. _compiled_source: -Custom Sources --------------- +Compiled Sources +---------------- It is often the case that one may wish to simulate a complex source distribution that is not possible to represent with the classes described above. For these @@ -305,7 +308,7 @@ below. #include "openmc/source.h" #include "openmc/particle.h" - class CustomSource : public openmc::Source + class CompiledSource : public openmc::Source { openmc::SourceSite sample(uint64_t* seed) const { @@ -327,9 +330,9 @@ below. } }; - extern "C" std::unique_ptr openmc_create_source(std::string parameters) + extern "C" std::unique_ptr openmc_create_source(std::string parameters) { - return std::make_unique(); + return std::make_unique(); } The above source creates monodirectional 14.08 MeV neutrons that are distributed @@ -356,19 +359,21 @@ OpenMC shared library. This can be done by writing a CMakeLists.txt file: target_link_libraries(source OpenMC::libopenmc) After running ``cmake`` and ``make``, you will have a libsource.so (or .dylib) -file in your build directory. Setting the :attr:`openmc.Source.library` -attribute to the path of this shared library will indicate that it should be -used for sampling source particles at runtime. +file in your build directory. You can then use this as an external source during +an OpenMC run by passing the path of the shared library to the +:class:`openmc.CompiledSource` class, which is then set as the +:attr:`Settings.source` attribute:: -.. _parameterized_custom_source: + settings.source = openmc.CompiledSource('libsource.so') -Custom Parameterized Sources ----------------------------- +.. _parameterized_compiled_source: -Some custom sources may have values (parameters) that can be changed between -runs. This is supported by using the ``openmc_create_source()`` function to -pass parameters defined in the :attr:`openmc.Source.parameters` attribute to -the source class when it is created: +Parameterized Compiled Sources +------------------------------ + +Some compiled sources may have values (parameters) that can be changed between +runs. This is supported by using the ``openmc_create_source()`` function to pass +parameters to the source class when it is created: .. code-block:: c++ @@ -377,9 +382,9 @@ the source class when it is created: #include "openmc/source.h" #include "openmc/particle.h" - class CustomSource : public openmc::Source { + class CompiledSource : public openmc::Source { public: - CustomSource(double energy) : energy_{energy} { } + CompiledSource(double energy) : energy_{energy} { } // Samples from an instance of this class. openmc::SourceSite sample(uint64_t* seed) const @@ -404,13 +409,16 @@ the source class when it is created: double energy_; }; - extern "C" std::unique_ptr openmc_create_source(std::string parameter) { + extern "C" std::unique_ptr openmc_create_source(std::string parameter) { double energy = std::stod(parameter); - return std::make_unique(energy); + return std::make_unique(energy); } -As with the basic custom source functionality, the custom source library -location must be provided in the :attr:`openmc.Source.library` attribute. +When creating an instance of the :class:`openmc.CompiledSource` class, you will +need to pass both the path of the shared library as well as the parameters as a +string, which gets passed down to the ``openmc_create_source()`` function:: + + settings.source = openmc.CompiledSource('libsource.so', '3.5e6') .. _usersguide_entropy: diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 28ce6795fc0..02328af58ca 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -318,9 +318,9 @@ The following tables show all valid scores: | |NJOY's HEATR module. | +----------------------+---------------------------------------------------+ |pulse-height |The energy deposited by an entire photon's history | - | |(including its progeny). Units are eV per source | + | |(including its progeny). Units are eV per source | | |particle. Note that this score can only be combined| - | |with a cell filter and an energy filter. | + | |with a cell filter and an energy filter. | +----------------------+---------------------------------------------------+ .. _usersguide_tally_normalization: @@ -337,7 +337,7 @@ it is usually straightforward to convert units if the source rate is known. For example, if the system being modeled includes a source that is emitting 10\ :sup:`4` neutrons per second, the tally results just need to be multipled by 10\ :sup:`4`. This can either be done manually or using the -:attr:`openmc.Source.strength` attribute. +:attr:`openmc.SourceBase.strength` attribute. For a :math:`k`\ -eigenvalue calculation, normalizing tally results is not as simple because the source rate is not actually known. Instead, we typically know diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 03b046bbbe1..65c68ad6457 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -22,6 +22,10 @@ class Distribution { public: virtual ~Distribution() = default; virtual double sample(uint64_t* seed) const = 0; + + //! Return integral of distribution + //! \return Integral of distribution + virtual double integral() const { return 1.0; }; }; using UPtrDist = unique_ptr; @@ -51,11 +55,13 @@ class DiscreteIndex { // Properties const vector& prob() const { return prob_; } const vector& alias() const { return alias_; } + double integral() const { return integral_; } private: vector prob_; //!< Probability of accepting the uniformly sampled bin, //!< mapped to alias method table vector alias_; //!< Alias table + double integral_; //!< Integral of distribution //! Normalize distribution so that probabilities sum to unity void normalize(); @@ -78,6 +84,8 @@ class Discrete : public Distribution { //! \return Sampled value double sample(uint64_t* seed) const override; + double integral() const override { return di_.integral(); }; + // Properties const vector& x() const { return x_; } const vector& prob() const { return di_.prob(); } @@ -219,17 +227,19 @@ class Tabular : public Distribution { //! \return Sampled value double sample(uint64_t* seed) const override; - // x property + // properties vector& x() { return x_; } const vector& x() const { return x_; } const vector& p() const { return p_; } Interpolation interp() const { return interp_; } + double integral() const override { return integral_; }; private: vector x_; //!< tabulated independent variable vector p_; //!< tabulated probability density vector c_; //!< cumulative distribution at tabulated values Interpolation interp_; //!< interpolation rule + double integral_; //!< Integral of distribution //! Initialize tabulated probability density function //! \param x Array of values for independent variable @@ -272,12 +282,15 @@ class Mixture : public Distribution { //! \return Sampled value double sample(uint64_t* seed) const override; + double integral() const override { return integral_; } + private: // Storrage for probability + distribution using DistPair = std::pair; vector - distribution_; //!< sub-distributions + cummulative probabilities + distribution_; //!< sub-distributions + cummulative probabilities + double integral_; //!< integral of distribution }; } // namespace openmc diff --git a/openmc/cell.py b/openmc/cell.py index 8572a640d09..9a744910c80 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -618,7 +618,10 @@ def plot(self, *args, **kwargs): Axes containing resulting image """ - return openmc.Universe(cells=[self]).plot(*args, **kwargs) + # Create dummy universe but preserve used_ids + u = openmc.Universe(cells=[self], universe_id=openmc.Universe.next_id + 1) + openmc.Universe.used_ids.remove(u.id) + return u.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/thermal.py b/openmc/data/thermal.py index b6e700a0cbd..994a8ed2ae8 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -676,6 +676,15 @@ def from_ace(cls, ace_or_filename, name=None): # here, because they are equiprobable, so the order # doesn't matter. mu.sort() + + # Older versions of NJOY had a bug, and the discrete + # scattering angles could sometimes be less than -1 or + # greater than 1. We check for this here, and warn users. + if mu[0] < -1. or mu[-1] > 1.: + warn('S(a,b) scattering angle for incident energy index ' + f'{i} and exit energy index {j} outside of the ' + 'interval [-1, 1].') + p_mu = 1. / n_mu * np.ones(n_mu) mu_ij = Discrete(mu, p_mu) mu_ij.c = np.cumsum(p_mu) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 9a1fc6a3fa4..1c9d96f679c 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -196,13 +196,13 @@ def reset_tally_means(self): """ self._rate_tally_means_cache = None - def get_material_rates(self, mat_id, nuc_index, rx_index): + def get_material_rates(self, mat_index, nuc_index, rx_index): """Return an array of reaction rates for a material Parameters ---------- - mat_id : int - Unique ID for the requested material + mat_index : int + Index for the material nuc_index : iterable of int Index for each nuclide in :attr:`nuclides` in the desired reaction rate matrix @@ -216,7 +216,7 @@ def get_material_rates(self, mat_id, nuc_index, rx_index): reaction rates in this material """ self._results_cache.fill(0.0) - full_tally_res = self.rate_tally_means[mat_id] + full_tally_res = self.rate_tally_means[mat_index] for i_tally, (i_nuc, i_rx) in enumerate(product(nuc_index, rx_index)): self._results_cache[i_nuc, i_rx] = full_tally_res[i_tally] diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index 75727840669..aa71942dfd8 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -5,8 +5,10 @@ """ +from __future__ import annotations +from collections.abc import Iterable import copy -from itertools import product +from typing import List, Set import numpy as np from uncertainties import ufloat @@ -37,14 +39,20 @@ class IndependentOperator(OpenMCOperator): .. versionadded:: 0.13.1 + .. versionchanged:: 0.13.4 + Arguments updated to include list of fluxes and microscopic cross + sections. + Parameters ---------- materials : openmc.Materials Materials to deplete. - micro_xs : MicroXS - One-group microscopic cross sections in [b]. If the - :class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation will - be run. + fluxes : list of numpy.ndarray + Flux in each group in [n-cm/src] for each domain + micros : list of MicroXS + Cross sections in [b] for each domain. If the + :class:`~openmc.deplete.MicroXS` object is empty, a decay-only + calculation will be run. chain_file : str Path to the depletion chain XML file. Defaults to ``openmc.config['chain_file']``. @@ -109,7 +117,8 @@ class IndependentOperator(OpenMCOperator): def __init__(self, materials, - micro_xs, + fluxes, + micros, chain_file=None, keff=None, normalization_mode='fission-q', @@ -120,7 +129,7 @@ def __init__(self, fission_yield_opts=None): # Validate micro-xs parameters check_type('materials', materials, openmc.Materials) - check_type('micro_xs', micro_xs, MicroXS) + check_type('micros', micros, Iterable, MicroXS) if keff is not None: check_type('keff', keff, tuple, float) keff = ufloat(*keff) @@ -132,10 +141,15 @@ def __init__(self, helper_kwargs = {'normalization_mode': normalization_mode, 'fission_yield_opts': fission_yield_opts} - cross_sections = micro_xs * 1e-24 + # Sort fluxes and micros in same order that materials get sorted + index_sort = np.argsort([mat.id for mat in materials]) + fluxes = [fluxes[i] for i in index_sort] + micros = [micros[i] for i in index_sort] + + self.fluxes = fluxes super().__init__( materials, - cross_sections, + micros, chain_file, prev_results, fission_q=fission_q, @@ -145,6 +159,7 @@ def __init__(self, @classmethod def from_nuclides(cls, volume, nuclides, + flux, micro_xs, chain_file=None, nuc_units='atom/b-cm', @@ -163,10 +178,11 @@ def from_nuclides(cls, volume, nuclides, nuclides : dict of str to float Dictionary with nuclide names as keys and nuclide concentrations as values. + flux : numpy.ndarray + Flux in each group in [n-cm/src] micro_xs : MicroXS - One-group microscopic cross sections in [b]. If the - :class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation - will be run. + Cross sections in [b]. If the :class:`~openmc.deplete.MicroXS` + object is empty, a decay-only calculation will be run. chain_file : str, optional Path to the depletion chain XML file. Defaults to ``openmc.config['chain_file']``. @@ -203,8 +219,11 @@ def from_nuclides(cls, volume, nuclides, """ check_type('nuclides', nuclides, dict, str) materials = cls._consolidate_nuclides_to_material(nuclides, nuc_units, volume) + fluxes = [flux] + micros = [micro_xs] return cls(materials, - micro_xs, + fluxes, + micros, chain_file, keff=keff, normalization_mode=normalization_mode, @@ -256,9 +275,9 @@ def _load_previous_results(self): self.prev_res.append(new_res) - def _get_nuclides_with_data(self, cross_sections): + def _get_nuclides_with_data(self, cross_sections: List[MicroXS]) -> Set[str]: """Finds nuclides with cross section data""" - return set(cross_sections.index) + return set(cross_sections[0].nuclides) class _IndependentRateHelper(ReactionRateHelper): """Class for generating one-group reaction rates with flux and @@ -285,7 +304,7 @@ class _IndependentRateHelper(ReactionRateHelper): """ - def __init__(self, op): + def __init__(self, op: IndependentOperator): rates = op.reaction_rates super().__init__(rates.n_nuc, rates.n_react) @@ -301,13 +320,13 @@ def reset_tally_means(self): """Unused in this case""" pass - def get_material_rates(self, mat_id, nuc_index, react_index): + def get_material_rates(self, mat_index, nuc_index, react_index): """Return 2D array of [nuclide, reaction] reaction rates Parameters ---------- - mat_id : int - Unique ID for the requested material + mat_index : int + Index for the material nuc_index : list of str Ordering of desired nuclides react_index : list of str @@ -315,20 +334,18 @@ def get_material_rates(self, mat_id, nuc_index, react_index): """ self._results_cache.fill(0.0) - # Get volume in units of [b-cm] - volume_b_cm = 1e24 * self._op.number.get_mat_volume(mat_id) + # Get flux and microscopic cross sections from operator + flux = self._op.fluxes[mat_index] + xs = self._op.cross_sections[mat_index] - for i_nuc, i_react in product(nuc_index, react_index): + for i_nuc in nuc_index: nuc = self.nuc_ind_map[i_nuc] - rx = self.rx_ind_map[i_react] - - # OK, this is kind of weird, but we multiply by volume in [b-cm] - # only because OpenMCOperator._calculate_reaction_rates has to - # divide it out later. It might make more sense to account for - # the source rate (flux) here rather than in the normalization - # helper. - self._results_cache[i_nuc, i_react] = \ - self._op.cross_sections[rx][nuc] * volume_b_cm + for i_rx in react_index: + rx = self.rx_ind_map[i_rx] + + # Determine reaction rate by multiplying xs in [b] by flux + # in [n-cm/src] to give [(reactions/src)*b-cm/atom] + self._results_cache[i_nuc, i_rx] = (xs[nuc, rx] * flux).sum() return self._results_cache diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index 050389cc26d..59285e213b2 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -1,15 +1,17 @@ """MicroXS module -A pandas.DataFrame storing microscopic cross section data with -nuclide names as row indices and reaction names as column indices. +A class for storing microscopic cross section data that can be used with the +IndependentOperator class for depletion. """ +from __future__ import annotations import tempfile +from typing import List, Tuple, Iterable, Optional, Union -from pandas import DataFrame, read_csv, Series +import pandas as pd import numpy as np -from openmc.checkvalue import check_type, check_value, check_iterable_type +from openmc.checkvalue import check_type, check_value, check_iterable_type, PathLike from openmc.exceptions import DataError from openmc import StatePoint import openmc @@ -20,168 +22,178 @@ _valid_rxns.append('fission') -class MicroXS(DataFrame): - """Microscopic cross section data for use in transport-independent depletion. - - .. versionadded:: 0.13.1 +def get_microxs_and_flux( + model: openmc.Model, + domains, + nuclides: Optional[Iterable[str]] = None, + reactions: Optional[Iterable[str]] = None, + energies: Optional[Union[Iterable[float], str]] = None, + chain_file: Optional[PathLike] = None, + run_kwargs=None + ) -> Tuple[List[np.ndarray], List[MicroXS]]: + """Generate a microscopic cross sections and flux from a Model + + .. versionadded:: 0.13.4 + + Parameters + ---------- + model : openmc.Model + OpenMC model object. Must contain geometry, materials, and settings. + domains : list of openmc.Material or openmc.Cell or openmc.Universe + Domains in which to tally reaction rates. + nuclides : list of str + Nuclides to get cross sections for. If not specified, all burnable + nuclides from the depletion chain file are used. + reactions : list of str + Reactions to get cross sections for. If not specified, all neutron + reactions listed in the depletion chain file are used. + energies : iterable of float or str + Energy group boundaries in [eV] or the name of the group structure + chain_file : str, optional + Path to the depletion chain XML file that will be used in depletion + simulation. Used to determine cross sections for materials not + present in the inital composition. Defaults to + ``openmc.config['chain_file']``. + run_kwargs : dict, optional + Keyword arguments passed to :meth:`openmc.Model.run` + + Returns + ------- + list of numpy.ndarray + Flux in each group in [n-cm/src] for each domain + list of MicroXS + Cross section data in [b] for each domain """ + # Save any original tallies on the model + original_tallies = model.tallies - @classmethod - def from_model(cls, - model, - domain, - nuclides=None, - reactions=None, - chain_file=None, - energy_bounds=(0, 20e6), - run_kwargs=None): - """Generate a one-group cross-section dataframe using OpenMC. - - Note that the ``openmc`` executable must be compiled. - - Parameters - ---------- - model : openmc.Model - OpenMC model object. Must contain geometry, materials, and settings. - domain : openmc.Material or openmc.Cell or openmc.Universe - Domain in which to tally reaction rates. - nuclides : list of str - Nuclides to get cross sections for. If not specified, all burnable - nuclides from the depletion chain file are used. - reactions : list of str - Reactions to get cross sections for. If not specified, all neutron - reactions listed in the depletion chain file are used. - chain_file : str, optional - Path to the depletion chain XML file that will be used in depletion - simulation. Used to determine cross sections for materials not - present in the inital composition. Defaults to - ``openmc.config['chain_file']``. - energy_bound : 2-tuple of float, optional - Bounds for the energy group. - run_kwargs : dict, optional - Keyword arguments passed to :meth:`openmc.model.Model.run` - - Returns - ------- - MicroXS - Cross section data in [b] - - """ - # Save any original tallies on the model - original_tallies = model.tallies - - # Determine what reactions and nuclides are available in chain + # Determine what reactions and nuclides are available in chain + if chain_file is None: + chain_file = openmc.config.get('chain_file') if chain_file is None: - chain_file = openmc.config.get('chain_file') - if chain_file is None: - raise DataError( - "No depletion chain specified and could not find depletion " - "chain in openmc.config['chain_file']" - ) - chain = Chain.from_xml(chain_file) - if reactions is None: - reactions = chain.reactions - if not nuclides: - cross_sections = _find_cross_sections(model) - nuclides_with_data = _get_nuclides_with_data(cross_sections) - nuclides = [nuc.name for nuc in chain.nuclides - if nuc.name in nuclides_with_data] - - # Set up the reaction rate and flux tallies - energy_filter = openmc.EnergyFilter(energy_bounds) - if isinstance(domain, openmc.Material): - domain_filter = openmc.MaterialFilter([domain]) - elif isinstance(domain, openmc.Cell): - domain_filter = openmc.CellFilter([domain]) - elif isinstance(domain, openmc.Universe): - domain_filter = openmc.UniverseFilter([domain]) + raise DataError( + "No depletion chain specified and could not find depletion " + "chain in openmc.config['chain_file']" + ) + chain = Chain.from_xml(chain_file) + if reactions is None: + reactions = chain.reactions + if not nuclides: + cross_sections = _find_cross_sections(model) + nuclides_with_data = _get_nuclides_with_data(cross_sections) + nuclides = [nuc.name for nuc in chain.nuclides + if nuc.name in nuclides_with_data] + + # Set up the reaction rate and flux tallies + if energies is None: + energy_filter = openmc.EnergyFilter([0.0, 100.0e6]) + elif isinstance(energies, str): + energy_filter = openmc.EnergyFilter.from_group_structure(energies) + else: + energy_filter = openmc.EnergyFilter(energies) + if isinstance(domains[0], openmc.Material): + domain_filter = openmc.MaterialFilter(domains) + elif isinstance(domains[0], openmc.Cell): + domain_filter = openmc.CellFilter(domains) + elif isinstance(domains[0], openmc.Universe): + domain_filter = openmc.UniverseFilter(domains) + else: + raise ValueError(f"Unsupported domain type: {type(domains[0])}") + + rr_tally = openmc.Tally(name='MicroXS RR') + rr_tally.filters = [domain_filter, energy_filter] + rr_tally.nuclides = nuclides + rr_tally.multiply_density = False + rr_tally.scores = reactions + + flux_tally = openmc.Tally(name='MicroXS flux') + flux_tally.filters = [domain_filter, energy_filter] + flux_tally.scores = ['flux'] + model.tallies = openmc.Tallies([rr_tally, flux_tally]) + + # create temporary run + with tempfile.TemporaryDirectory() as temp_dir: + if run_kwargs is None: + run_kwargs = {} else: - raise ValueError(f"Unsupported domain type: {type(domain)}") - - rr_tally = openmc.Tally(name='MicroXS RR') - rr_tally.filters = [domain_filter, energy_filter] - rr_tally.nuclides = nuclides - rr_tally.multiply_density = False - rr_tally.scores = reactions + run_kwargs = dict(run_kwargs) + run_kwargs.setdefault('cwd', temp_dir) + statepoint_path = model.run(**run_kwargs) - flux_tally = openmc.Tally(name='MicroXS flux') - flux_tally.filters = [domain_filter, energy_filter] - flux_tally.scores = ['flux'] - tallies = openmc.Tallies([rr_tally, flux_tally]) + with StatePoint(statepoint_path) as sp: + rr_tally = sp.tallies[rr_tally.id] + rr_tally._read_results() + flux_tally = sp.tallies[flux_tally.id] + flux_tally._read_results() - model.tallies = tallies + # Get reaction rates and flux values + reaction_rates = rr_tally.get_reshaped_data() # (domains, groups, nuclides, reactions) + flux = flux_tally.get_reshaped_data() # (domains, groups, 1, 1) - # create temporary run - with tempfile.TemporaryDirectory() as temp_dir: - if run_kwargs is None: - run_kwargs = {} - run_kwargs.setdefault('cwd', temp_dir) - statepoint_path = model.run(**run_kwargs) + # Make energy groups last dimension + reaction_rates = np.moveaxis(reaction_rates, 1, -1) # (domains, nuclides, reactions, groups) + flux = np.moveaxis(flux, 1, -1) # (domains, 1, 1, groups) - with StatePoint(statepoint_path) as sp: - rr_tally = sp.tallies[rr_tally.id] - rr_tally._read_results() - flux_tally = sp.tallies[flux_tally.id] - flux_tally._read_results() + # Divide RR by flux to get microscopic cross sections + xs = np.empty_like(reaction_rates) # (domains, nuclides, reactions, groups) + d, _, _, g = np.nonzero(flux) + xs[d, ..., g] = reaction_rates[d, ..., g] / flux[d, :, :, g] - # Get reaction rates and flux values - reaction_rates = rr_tally.mean.sum(axis=0) # (nuclides, reactions) - flux = flux_tally.mean[0, 0, 0] + # Reset tallies + model.tallies = original_tallies - # Divide RR by flux to get microscopic cross sections - xs = reaction_rates / flux + # Create lists where each item corresponds to one domain + fluxes = list(flux.squeeze((1, 2))) + micros = [MicroXS(xs_i, nuclides, reactions) for xs_i in xs] + return fluxes, micros - # Build Series objects - series = {} - for i, rx in enumerate(reactions): - series[rx] = Series(xs[..., i], index=rr_tally.nuclides) - # Revert to the original tallies and materials - model.tallies = original_tallies - - return cls(series).rename_axis('nuclide') +class MicroXS: + """Microscopic cross section data for use in transport-independent depletion. - @classmethod - def from_array(cls, nuclides, reactions, data): - """ - Creates a ``MicroXS`` object from arrays. + .. versionadded:: 0.13.1 - Parameters - ---------- - nuclides : list of str - List of nuclide symbols for that have data for at least one - reaction. - reactions : list of str - List of reactions. All reactions must match those in - :data:`openmc.deplete.chain.REACTIONS` - data : ndarray of floats - Array containing one-group microscopic cross section values for - each nuclide and reaction. Cross section values are assumed to be - in [b]. - - Returns - ------- - MicroXS - """ + .. versionchanged:: 0.13.4 + Class was heavily refactored and no longer subclasses pandas.DataFrame. + + Parameters + ---------- + data : numpy.ndarray of floats + 3D array containing microscopic cross section values for each + nuclide, reaction, and energy group. Cross section values are assumed to + be in [b], and indexed by [nuclide, reaction, energy group] + nuclides : list of str + List of nuclide symbols for that have data for at least one + reaction. + reactions : list of str + List of reactions. All reactions must match those in + :data:`openmc.deplete.chain.REACTIONS` + """ + def __init__(self, data: np.ndarray, nuclides: List[str], reactions: List[str]): # Validate inputs - if data.shape != (len(nuclides), len(reactions)): + if data.shape[:2] != (len(nuclides), len(reactions)): raise ValueError( f'Nuclides list of length {len(nuclides)} and ' f'reactions array of length {len(reactions)} do not ' f'match dimensions of data array of shape {data.shape}') + check_iterable_type('nuclides', nuclides, str) + check_iterable_type('reactions', reactions, str) + check_type('data', data, np.ndarray, expected_iter_type=float) + for reaction in reactions: + check_value('reactions', reaction, _valid_rxns) - cls._validate_micro_xs_inputs( - nuclides, reactions, data) - micro_xs = cls(index=nuclides, columns=reactions, data=data) + self.data = data + self.nuclides = nuclides + self.reactions = reactions - return micro_xs + # TODO: Add a classmethod for generating MicroXS directly from cross section + # data using a known flux spectrum @classmethod def from_csv(cls, csv_file, **kwargs): - """ - Load a ``MicroXS`` object from a ``.csv`` file. + """Load data from a comma-separated values (csv) file. Parameters ---------- @@ -199,17 +211,37 @@ def from_csv(cls, csv_file, **kwargs): if 'float_precision' not in kwargs: kwargs['float_precision'] = 'round_trip' - micro_xs = cls(read_csv(csv_file, index_col=0, **kwargs)) + df = pd.read_csv(csv_file, **kwargs) + df.set_index(['nuclides', 'reactions', 'groups'], inplace=True) + nuclides = list(df.index.unique(level='nuclides')) + reactions = list(df.index.unique(level='reactions')) + groups = list(df.index.unique(level='groups')) + shape = (len(nuclides), len(reactions), len(groups)) + data = df.values.reshape(shape) + return cls(data, nuclides, reactions) - cls._validate_micro_xs_inputs(list(micro_xs.index), - list(micro_xs.columns), - micro_xs.to_numpy()) - return micro_xs + def __getitem__(self, index): + nuc, rx = index + i_nuc = self.nuclides.index(nuc) + i_rx = self.reactions.index(rx) + return self.data[i_nuc, i_rx] + + def to_csv(self, *args, **kwargs): + """Write data to a comma-separated values (csv) file + + Parameters + ---------- + *args + Positional arguments passed to :meth:`pandas.DataFrame.to_csv` + **kwargs + Keyword arguments passed to :meth:`pandas.DataFrame.to_csv` + + """ + groups = self.data.shape[2] + multi_index = pd.MultiIndex.from_product( + [self.nuclides, self.reactions, range(1, groups + 1)], + names=['nuclides', 'reactions', 'groups'] + ) + df = pd.DataFrame({'xs': self.data.flatten()}, index=multi_index) + df.to_csv(*args, **kwargs) - @staticmethod - def _validate_micro_xs_inputs(nuclides, reactions, data): - check_iterable_type('nuclides', nuclides, str) - check_iterable_type('reactions', reactions, str) - check_type('data', data, np.ndarray, expected_iter_type=float) - for reaction in reactions: - check_value('reactions', reaction, _valid_rxns) diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index d3b9ef911cd..319d9870122 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -6,6 +6,7 @@ """ from abc import abstractmethod +from typing import List, Tuple, Dict import numpy as np @@ -31,9 +32,9 @@ class OpenMCOperator(TransportOperator): ---------- materials : openmc.Materials List of all materials in the model - cross_sections : str or pandas.DataFrame - Path to continuous energy cross section library, or object containing - one-group cross-sections. + cross_sections : str or list of MicroXS + Path to continuous energy cross section library, or list of objects + containing cross sections. chain_file : str, optional Path to the depletion chain XML file. Defaults to openmc.config['chain_file']. @@ -60,9 +61,9 @@ class OpenMCOperator(TransportOperator): ---------- materials : openmc.Materials All materials present in the model - cross_sections : str or MicroXS - Path to continuous energy cross section library, or object - containing one-group cross-sections. + cross_sections : str or list of MicroXS + Path to continuous energy cross section library, or list of objects + containing cross sections. output_dir : pathlib.Path Path to output directory to save results. round_number : bool @@ -164,7 +165,7 @@ def _differentiate_burnable_mats(self): """Assign distribmats for each burnable material""" pass - def _get_burnable_mats(self): + def _get_burnable_mats(self) -> Tuple[List[str], Dict[str, float], List[str]]: """Determine depletable materials, volumes, and nuclides Returns diff --git a/openmc/deplete/reaction_rates.py b/openmc/deplete/reaction_rates.py index b806dfad84b..4aaf829a9a5 100644 --- a/openmc/deplete/reaction_rates.py +++ b/openmc/deplete/reaction_rates.py @@ -2,6 +2,8 @@ An ndarray to store reaction rates with string, integer, or slice indexing. """ +from typing import Dict + import numpy as np @@ -51,6 +53,10 @@ class ReactionRates(np.ndarray): # the __array_finalize__ method (discussed here: # https://docs.scipy.org/doc/numpy/user/basics.subclassing.html) + index_mat: Dict[str, int] + index_nuc: Dict[str, int] + index_rx: Dict[str, int] + def __new__(cls, local_mats, nuclides, reactions, from_results=False): # Create appropriately-sized zeroed-out ndarray shape = (len(local_mats), len(nuclides), len(reactions)) diff --git a/openmc/deplete/stepresult.py b/openmc/deplete/stepresult.py index b8648fb0ba4..5b4b8d2d3f9 100644 --- a/openmc/deplete/stepresult.py +++ b/openmc/deplete/stepresult.py @@ -214,11 +214,23 @@ def get_material(self, mat_id): ------- openmc.Material Equivalent material + + Raises + ------ + KeyError + If specified material ID is not found in the StepResult + """ with warnings.catch_warnings(): warnings.simplefilter('ignore', openmc.IDWarning) material = openmc.Material(material_id=int(mat_id)) - vol = self.volume[mat_id] + try: + vol = self.volume[mat_id] + except KeyError as e: + raise KeyError( + f'mat_id {mat_id} not found in StepResult. Available mat_id ' + f'values are {list(self.volume.keys())}' + ) from e for nuc, _ in sorted(self.index_nuc.items(), key=lambda x: x[1]): atoms = self[0, mat_id, nuc] if atoms < 0.0: diff --git a/openmc/geometry.py b/openmc/geometry.py index 165b7f1a282..82a3fe61c5c 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -734,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/mesh.py b/openmc/mesh.py index 5313ebb398b..60a1c0c9398 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1188,7 +1188,8 @@ class CylindricalMesh(StructuredMesh): 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. + 1-D array of mesh boundary points along the z-axis relative to the + origin. 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. @@ -1207,8 +1208,7 @@ class CylindricalMesh(StructuredMesh): name : str Name of the mesh dimension : Iterable of int - The number of mesh cells in each direction (r_grid, - phi_grid, z_grid). + 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 @@ -1218,7 +1218,8 @@ class CylindricalMesh(StructuredMesh): 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. z_grid : numpy.ndarray - 1-D array of mesh boundary points along the z-axis. + 1-D array of mesh boundary points along the z-axis relative to the + origin. origin : numpy.ndarray 1-D array of length 3 the (x,y,z) origin of the mesh in cartesian coordinates @@ -1314,13 +1315,12 @@ 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] + self.origin[2] + self.z_grid[0] )) @property @@ -1413,8 +1413,7 @@ def from_domain( (openmc.Cell, openmc.Region, openmc.Universe, openmc.Geometry), ) - - # loaded once to avoid reading h5m file repeatedly + # loaded once to avoid recalculating bounding box cached_bb = domain.bounding_box max_bounding_box_radius = max( [ @@ -1439,8 +1438,14 @@ def from_domain( cached_bb[1][2], num=dimension[2]+1 ) + origin = (cached_bb.center[0], cached_bb.center[1], z_grid[0]) mesh = cls( - r_grid=r_grid, z_grid=z_grid, phi_grid=phi_grid, mesh_id=mesh_id, name=name + r_grid=r_grid, + z_grid=z_grid, + phi_grid=phi_grid, + mesh_id=mesh_id, + name=name, + origin=origin ) return mesh 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/distribution.cpp b/src/distribution.cpp index 10d404c10f4..9c76147ee27 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -104,10 +104,12 @@ size_t DiscreteIndex::sample(uint64_t* seed) const void DiscreteIndex::normalize() { - // Renormalize density function so that it sums to unity - double norm = std::accumulate(prob_.begin(), prob_.end(), 0.0); + // Renormalize density function so that it sums to unity. Note that we save + // the integral of the distribution so that if it is used as part of another + // distribution (e.g., Mixture), we know its relative strength. + integral_ = std::accumulate(prob_.begin(), prob_.end(), 0.0); for (auto& p_i : prob_) { - p_i /= norm; + p_i /= integral_; } } @@ -300,10 +302,13 @@ void Tabular::init( } } - // Normalize density and distribution functions + // Normalize density and distribution functions. Note that we save the + // integral of the distribution so that if it is used as part of another + // distribution (e.g., Mixture), we know its relative strength. + integral_ = c_[n - 1]; for (int i = 0; i < n; ++i) { - p_[i] = p_[i] / c_[n - 1]; - c_[i] = c_[i] / c_[n - 1]; + p_[i] = p_[i] / integral_; + c_[i] = c_[i] / integral_; } } @@ -379,12 +384,14 @@ Mixture::Mixture(pugi::xml_node node) if (!pair.child("dist")) fatal_error("Mixture pair element does not have a distribution."); - // cummulative sum of probybilities - cumsum += std::stod(pair.attribute("probability").value()); + // cummulative sum of probabilities + double p = std::stod(pair.attribute("probability").value()); - // Save cummulative probybility and distrubution - distribution_.push_back( - std::make_pair(cumsum, distribution_from_xml(pair.child("dist")))); + // Save cummulative probability and distribution + auto dist = distribution_from_xml(pair.child("dist")); + cumsum += p * dist->integral(); + + distribution_.push_back(std::make_pair(cumsum, std::move(dist))); } // Normalize cummulative probabilities to 1 diff --git a/src/geometry.cpp b/src/geometry.cpp index 29ca8b1bf8f..47152e17295 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) { @@ -419,13 +419,14 @@ BoundaryInfo distance_to_boundary(Particle& p) double& d = info.distance; if (d_surf < d_lat - FP_COINCIDENT) { if (d == INFINITY || (d - d_surf) / d >= FP_REL_PRECISION) { + // Update closest distance d = d_surf; // If the cell is not simple, it is possible that both the negative and // positive half-space were given in the region specification. Thus, we // have to explicitly check which half-space the particle would be // traveling into if the surface is crossed - if (c.is_simple()) { + if (c.is_simple() || d == INFTY) { info.surface_index = level_surf_cross; } else { Position r_hit = r + d_surf * u; @@ -472,7 +473,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/nuclide.cpp b/src/nuclide.cpp index 8f73a530930..6b70ad9e323 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -857,9 +857,8 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const // This guarantees the randomness and, at the same time, makes sure we // reuse random numbers for the same nuclide at different temperatures, // therefore preserving correlation of temperature in probability tables. - p.stream() = STREAM_URR_PTABLE; - double r = future_prn(static_cast(index_), *p.current_seed()); - p.stream() = STREAM_TRACKING; + double r = + future_prn(static_cast(index_), p.seeds(STREAM_URR_PTABLE)); // Warning: this assumes row-major order of cdf_values_ int i_low = upper_bound_index(&urr.cdf_values_(i_energy, 0), diff --git a/src/particle.cpp b/src/particle.cpp index a40ea7024c4..ceef073c480 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -147,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())); @@ -156,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. @@ -392,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())); @@ -400,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(); } @@ -456,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); @@ -535,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 28c316c6878..7101a2969ba 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -150,9 +150,7 @@ void sample_neutron_reaction(Particle& p) // Advance URR seed stream 'N' times after energy changes if (p.E() != p.E_last()) { - p.stream() = STREAM_URR_PTABLE; - advance_prn_seed(data::nuclides.size(), p.current_seed()); - p.stream() = STREAM_TRACKING; + advance_prn_seed(data::nuclides.size(), &p.seeds(STREAM_URR_PTABLE)); } // Play russian roulette if survival biasing is turned on 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/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/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/tests/micro_xs_simple.csv b/tests/micro_xs_simple.csv index 71981c718c8..146896aa290 100644 --- a/tests/micro_xs_simple.csv +++ b/tests/micro_xs_simple.csv @@ -1,13 +1,25 @@ -nuclide,"(n,gamma)",fission -U234,22.231989822002454,0.4962074466374984 -U235,10.479008971197121,48.41787337164606 -U238,0.8673334105437321,0.1046788058876236 -U236,8.651710446071224,0.31948392400019293 -O16,7.497851000107522e-05,0.0 -O17,0.0004079227797153372,0.0 -I135,6.842395323713929,0.0 -Xe135,227463.8642699061,0.0 -Xe136,0.023178960347535887,0.0 -Cs135,2.1721665580713623,0.0 -Gd157,12786.099392370175,0.0 -Gd156,3.4006085445846983,0.0 +nuclides,reactions,groups,xs +U234,"(n,gamma)",1,22.23198982200245 +U234,fission,1,0.4962074466374984 +U235,"(n,gamma)",1,10.47900897119712 +U235,fission,1,48.41787337164606 +U238,"(n,gamma)",1,0.8673334105437321 +U238,fission,1,0.1046788058876236 +U236,"(n,gamma)",1,8.651710446071224 +U236,fission,1,0.3194839240001929 +O16,"(n,gamma)",1,7.497851000107522e-05 +O16,fission,1,0.0 +O17,"(n,gamma)",1,0.0004079227797153 +O17,fission,1,0.0 +I135,"(n,gamma)",1,6.842395323713929 +I135,fission,1,0.0 +Xe135,"(n,gamma)",1,227463.8642699061 +Xe135,fission,1,0.0 +Xe136,"(n,gamma)",1,0.0231789603475358 +Xe136,fission,1,0.0 +Cs135,"(n,gamma)",1,2.1721665580713623 +Cs135,fission,1,0.0 +Gd157,"(n,gamma)",1,12786.099392370175 +Gd157,fission,1,0.0 +Gd156,"(n,gamma)",1,3.4006085445846983 +Gd156,fission,1,0.0 diff --git a/tests/regression_tests/deplete_no_transport/test.py b/tests/regression_tests/deplete_no_transport/test.py index 64769c0b0e8..54c29cd5b70 100644 --- a/tests/regression_tests/deplete_no_transport/test.py +++ b/tests/regression_tests/deplete_no_transport/test.py @@ -36,13 +36,16 @@ def chain_file(): return Path(__file__).parents[2] / 'chain_simple.xml' -@pytest.mark.parametrize("multiproc, from_nuclides, normalization_mode, power, flux", [ - (True, True, 'source-rate', None, 1164719970082145.0), - (False, True, 'source-rate', None, 1164719970082145.0), +neutron_per_cm2_sec = 1164719970082145.0 + + +@pytest.mark.parametrize("multiproc, from_nuclides, normalization_mode, power, source_rate", [ + (True, True, 'source-rate', None, 1.0), + (False, True, 'source-rate', None, 1.0), (True, True, 'fission-q', 174, None), (False, True, 'fission-q', 174, None), - (True, False, 'source-rate', None, 1164719970082145.0), - (False, False, 'source-rate', None, 1164719970082145.0), + (True, False, 'source-rate', None, 1.0), + (False, False, 'source-rate', None, 1.0), (True, False, 'fission-q', 174, None), (False, False, 'fission-q', 174, None)]) def test_against_self(run_in_tmpdir, @@ -53,7 +56,7 @@ def test_against_self(run_in_tmpdir, from_nuclides, normalization_mode, power, - flux): + source_rate): """Transport free system test suite. Runs an OpenMC transport-free depletion calculation and verifies @@ -61,8 +64,10 @@ def test_against_self(run_in_tmpdir, """ # Create operator + flux = neutron_per_cm2_sec * fuel.volume op = _create_operator(from_nuclides, fuel, + flux, micro_xs, chain_file, normalization_mode) @@ -75,12 +80,12 @@ def test_against_self(run_in_tmpdir, openmc.deplete.PredictorIntegrator(op, dt, power=power, - source_rates=flux, + source_rates=source_rate, timestep_units='s').integrate() # Get path to test and reference results path_test = op.output_dir / 'depletion_results.h5' - if flux is not None: + if power is None: ref_path = 'test_reference_source_rate.h5' else: ref_path = 'test_reference_fission_q.h5' @@ -99,8 +104,8 @@ def test_against_self(run_in_tmpdir, _assert_same_mats(res_test, res_ref) tol = 1.0e-14 - assert_atoms_equal(res_test, res_ref, tol) - assert_reaction_rates_equal(res_test, res_ref, tol) + assert_atoms_equal(res_ref, res_test, tol) + assert_reaction_rates_equal(res_ref, res_test, tol) @pytest.mark.parametrize("multiproc, dt, time_units, time_type, atom_tol, rx_tol ", [ @@ -123,7 +128,8 @@ def test_against_coupled(run_in_tmpdir, atom_tol, rx_tol): # Create operator - op = _create_operator(False, fuel, micro_xs, chain_file, 'fission-q') + flux = neutron_per_cm2_sec * fuel.volume + op = _create_operator(False, fuel, flux, micro_xs, chain_file, 'fission-q') # Power and timesteps dt = [dt] # single step @@ -151,12 +157,13 @@ def test_against_coupled(run_in_tmpdir, # Assert same mats _assert_same_mats(res_test, res_ref) - assert_atoms_equal(res_test, res_ref, atom_tol) - assert_reaction_rates_equal(res_test, res_ref, rx_tol) + assert_atoms_equal(res_ref, res_test, atom_tol) + assert_reaction_rates_equal(res_ref, res_test, rx_tol) def _create_operator(from_nuclides, fuel, + flux, micro_xs, chain_file, normalization_mode): @@ -167,13 +174,15 @@ def _create_operator(from_nuclides, op = IndependentOperator.from_nuclides(fuel.volume, nuclides, + flux, micro_xs, chain_file, normalization_mode=normalization_mode) else: op = IndependentOperator(openmc.Materials([fuel]), - micro_xs, + [flux], + [micro_xs], chain_file, normalization_mode=normalization_mode) diff --git a/tests/regression_tests/microxs/test.py b/tests/regression_tests/microxs/test.py index 3d67ee8a38d..f4b3649fd5b 100644 --- a/tests/regression_tests/microxs/test.py +++ b/tests/regression_tests/microxs/test.py @@ -4,7 +4,7 @@ import numpy as np import pytest import openmc -from openmc.deplete import MicroXS +from openmc.deplete import MicroXS, get_microxs_and_flux from tests.regression_tests import config @@ -47,13 +47,13 @@ def model(): def test_from_model(model): - fuel = model.materials[0] + domains = model.materials[:1] nuclides = ['U234', 'U235', 'U238', 'U236', 'O16', 'O17', 'I135', 'Xe135', 'Xe136', 'Cs135', 'Gd157', 'Gd156'] - test_xs = MicroXS.from_model(model, fuel, nuclides, chain_file=CHAIN_FILE) + _, test_xs = get_microxs_and_flux(model, domains, nuclides, chain_file=CHAIN_FILE) if config['update']: - test_xs.to_csv('test_reference.csv') + test_xs[0].to_csv('test_reference.csv') ref_xs = MicroXS.from_csv('test_reference.csv') - np.testing.assert_allclose(test_xs, ref_xs, rtol=1e-11) + np.testing.assert_allclose(test_xs[0].data, ref_xs.data, rtol=1e-11) diff --git a/tests/regression_tests/microxs/test_reference.csv b/tests/regression_tests/microxs/test_reference.csv index 0fc41862d66..3e679f98071 100644 --- a/tests/regression_tests/microxs/test_reference.csv +++ b/tests/regression_tests/microxs/test_reference.csv @@ -1,13 +1,25 @@ -nuclide,"(n,gamma)",fission -U234,21.418670317831197,0.5014588470882195 -U235,10.343944102483244,47.46718472611891 -U238,0.8741166723597251,0.10829568455139126 -U236,9.083486784689326,0.3325287927011428 -O16,7.548646353912453e-05,0.0 -O17,0.0004018486221310307,0.0 -I135,6.6912565089429235,0.0 -Xe135,223998.64185667288,0.0 -Xe136,0.022934362666193576,0.0 -Cs135,2.28453952223533,0.0 -Gd157,12582.079620036275,0.0 -Gd156,2.9421127515332417,0.0 +nuclides,reactions,groups,xs +U234,"(n,gamma)",1,21.418670317831076 +U234,fission,1,0.5014588470882162 +U235,"(n,gamma)",1,10.343944102483215 +U235,fission,1,47.46718472611895 +U238,"(n,gamma)",1,0.8741166723597229 +U238,fission,1,0.10829568455139067 +U236,"(n,gamma)",1,9.08348678468935 +U236,fission,1,0.3325287927011424 +O16,"(n,gamma)",1,7.548646353912426e-05 +O16,fission,1,0.0 +O17,"(n,gamma)",1,0.00040184862213103105 +O17,fission,1,0.0 +I135,"(n,gamma)",1,6.691256508942912 +I135,fission,1,0.0 +Xe135,"(n,gamma)",1,223998.6418566729 +Xe135,fission,1,0.0 +Xe136,"(n,gamma)",1,0.022934362666193503 +Xe136,fission,1,0.0 +Cs135,"(n,gamma)",1,2.2845395222353204 +Cs135,fission,1,0.0 +Gd157,"(n,gamma)",1,12582.07962003624 +Gd157,fission,1,0.0 +Gd156,"(n,gamma)",1,2.942112751533234 +Gd156,fission,1,0.0 diff --git a/tests/unit_tests/test_cell.py b/tests/unit_tests/test_cell.py index 1df3df9c55f..888a0ef88f3 100644 --- a/tests/unit_tests/test_cell.py +++ b/tests/unit_tests/test_cell.py @@ -345,3 +345,21 @@ def test_rotation_from_xml(rotation): elem, {s.id: s}, {'void': None}, openmc.Universe ) np.testing.assert_allclose(new_cell.rotation, cell.rotation) + + +def test_plot(run_in_tmpdir): + zcyl = openmc.ZCylinder() + c = openmc.Cell(region=-zcyl) + + # create a universe before the plot + u_before = openmc.Universe() + + # create a plot of the cell + c.plot() + + # create a universe after the plot + u_after = openmc.Universe() + + # ensure that calling the plot method doesn't + # affect the universe ID space + assert u_before.id + 1 == u_after.id diff --git a/tests/unit_tests/test_deplete_decay_products.py b/tests/unit_tests/test_deplete_decay_products.py index 34ca1b067e8..751a96ca6e6 100644 --- a/tests/unit_tests/test_deplete_decay_products.py +++ b/tests/unit_tests/test_deplete_decay_products.py @@ -1,4 +1,5 @@ import openmc.deplete +import numpy as np import pytest @@ -16,11 +17,14 @@ def test_deplete_decay_products(run_in_tmpdir): """) + # Create MicroXS object with no cross sections + micro_xs = openmc.deplete.MicroXS(np.empty((0, 0)), [], []) + # Create depletion operator with no reactions - micro_xs = openmc.deplete.MicroXS() op = openmc.deplete.IndependentOperator.from_nuclides( volume=1.0, nuclides={'Li5': 1.0}, + flux=0.0, micro_xs=micro_xs, chain_file='test_chain.xml', normalization_mode='source-rate' diff --git a/tests/unit_tests/test_deplete_independent_operator.py b/tests/unit_tests/test_deplete_independent_operator.py index f6cc7e2000d..813ac06de80 100644 --- a/tests/unit_tests/test_deplete_independent_operator.py +++ b/tests/unit_tests/test_deplete_independent_operator.py @@ -4,14 +4,13 @@ from pathlib import Path -import pytest from openmc.deplete import IndependentOperator, MicroXS from openmc import Material, Materials -import numpy as np CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml" ONE_GROUP_XS = Path(__file__).parents[1] / "micro_xs_simple.csv" + def test_operator_init(): """The test uses a temporary dummy chain. This file will be removed at the end of the test, and only contains a depletion_chain node.""" @@ -22,15 +21,18 @@ def test_operator_init(): 'U236': 4.5724195495061115e+18, 'O16': 4.639065406771322e+22, 'O17': 1.7588724018066158e+19} + flux = 1.0 micro_xs = MicroXS.from_csv(ONE_GROUP_XS) IndependentOperator.from_nuclides( - volume, nuclides, micro_xs, CHAIN_PATH, nuc_units='atom/cm3') + volume, nuclides, flux, micro_xs, CHAIN_PATH, nuc_units='atom/cm3') fuel = Material(name="uo2") fuel.add_element("U", 1, percent_type="ao", enrichment=4.25) fuel.add_element("O", 2) fuel.set_density("g/cc", 10.4) - fuel.depletable=True + fuel.depletable = True fuel.volume = 1 materials = Materials([fuel]) - IndependentOperator(materials, micro_xs, CHAIN_PATH) + fluxes = [1.0] + micros = [micro_xs] + IndependentOperator(materials, fluxes, micros, CHAIN_PATH) diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py index 557082eda4f..582c584654d 100644 --- a/tests/unit_tests/test_deplete_microxs.py +++ b/tests/unit_tests/test_deplete_microxs.py @@ -43,17 +43,18 @@ def test_from_array(): [0., 0.], [0., 0.1], [0., 0.1]]) + data.shape = (12, 2, 1) - MicroXS.from_array(nuclides, reactions, data) + MicroXS(data, nuclides, reactions) with pytest.raises(ValueError, match=r'Nuclides list of length \d* and ' r'reactions array of length \d* do not ' - r'match dimensions of data array of shape \(\d*\,d*\)'): - MicroXS.from_array(nuclides, reactions, data[:, 0]) + r'match dimensions of data array of shape \(\d*\, \d*\)'): + MicroXS(data[:, 0], nuclides, reactions) def test_csv(): ref_xs = MicroXS.from_csv(ONE_GROUP_XS) ref_xs.to_csv('temp_xs.csv') temp_xs = MicroXS.from_csv('temp_xs.csv') - assert np.all(ref_xs == temp_xs) + assert np.all(ref_xs.data == temp_xs.data) remove('temp_xs.csv') diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 480104154d5..77b40b44fe8 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -56,21 +56,26 @@ def test_cylindrical_mesh_bounding_box(): 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)) + np.testing.assert_array_equal(mesh.lower_left, (-1, -1, 0.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.lower_left, (-1, -1, 0.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)) + np.testing.assert_array_equal(mesh.lower_left, (2, 4, 7.1)) 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.lower_left, (2, 4, 7.1)) np.testing.assert_array_equal(bb.upper_right, (4, 6, 8)) + # changing z grid to contain negative numbers + mesh.z_grid = [-10, 0, 10] + np.testing.assert_array_equal(mesh.lower_left, (2, 4, -3)) + np.testing.assert_array_equal(mesh.upper_right, (4, 6, 17)) + def test_spherical_mesh_bounding_box(): # test with mesh at origin (0, 0, 0) diff --git a/tests/unit_tests/test_mesh_from_domain.py b/tests/unit_tests/test_mesh_from_domain.py index 015c8f88f0d..c75e8fd078c 100644 --- a/tests/unit_tests/test_mesh_from_domain.py +++ b/tests/unit_tests/test_mesh_from_domain.py @@ -18,10 +18,11 @@ def test_reg_mesh_from_cell(): def test_cylindrical_mesh_from_cell(): """Tests a CylindricalMesh can be made from a Cell and the specified - dimensions are propagated through. Cell is not centralized""" + dimensions are propagated through.""" + # Cell is not centralized on Z axis cy_surface = openmc.ZCylinder(r=50) - z_surface_1 = openmc.ZPlane(z0=30) - z_surface_2 = openmc.ZPlane(z0=0) + z_surface_1 = openmc.ZPlane(z0=40) + z_surface_2 = openmc.ZPlane(z0=10) cell = openmc.Cell(region=-cy_surface & -z_surface_1 & +z_surface_2) mesh = openmc.CylindricalMesh.from_domain(domain=cell, dimension=[2, 4, 3]) @@ -29,7 +30,27 @@ def test_cylindrical_mesh_from_cell(): assert np.array_equal(mesh.dimension, (2, 4, 3)) assert np.array_equal(mesh.r_grid, [0., 25., 50.]) assert np.array_equal(mesh.phi_grid, [0., 0.5*np.pi, np.pi, 1.5*np.pi, 2.*np.pi]) - assert np.array_equal(mesh.z_grid, [0., 10., 20., 30.]) + assert np.array_equal(mesh.z_grid, [10., 20., 30., 40.]) + assert np.array_equal(mesh.origin, [0., 0., 10.]) + + # Cell is not centralized on Z or X axis + cy_surface = openmc.ZCylinder(r=50, x0=100) + cell = openmc.Cell(region=-cy_surface & -z_surface_1 & +z_surface_2) + mesh = openmc.CylindricalMesh.from_domain(domain=cell, dimension=[1, 1, 1]) + + assert isinstance(mesh, openmc.CylindricalMesh) + assert np.array_equal(mesh.dimension, (1, 1, 1)) + assert np.array_equal(mesh.r_grid, [0., 150.]) + assert np.array_equal(mesh.origin, [100., 0., 10.]) + + # Cell is not centralized on Z, X or Y axis + cy_surface = openmc.ZCylinder(r=50, x0=100, y0=170) + cell = openmc.Cell(region=-cy_surface & -z_surface_1 & +z_surface_2) + mesh = openmc.CylindricalMesh.from_domain(domain=cell, dimension=[1, 1, 1]) + + assert isinstance(mesh, openmc.CylindricalMesh) + assert np.array_equal(mesh.r_grid, [0., 220.]) + assert np.array_equal(mesh.origin, [100., 170., 10.]) def test_reg_mesh_from_region(): diff --git a/tests/unit_tests/test_no_visible_boundary.py b/tests/unit_tests/test_no_visible_boundary.py new file mode 100644 index 00000000000..8f07b8da9e7 --- /dev/null +++ b/tests/unit_tests/test_no_visible_boundary.py @@ -0,0 +1,30 @@ +import openmc + + +def test_no_visible_boundary(run_in_tmpdir): + copper = openmc.Material() + copper.add_nuclide('Cu63', 1.0) + copper.set_density('g/cm3', 0.3) + air = openmc.Material() + air.add_nuclide('N14', 1.0) + air.set_density('g/cm3', 0.0012) + + # Create a simple model of a neutron source directly impinging on a thin + # disc of copper. Neutrons leaving the back of the disc see no surfaces in + # front of them. + disc = openmc.model.RightCircularCylinder((0., 0., 1.), 0.1, 1.2) + box = openmc.rectangular_prism(width=10, height=10, boundary_type='vacuum') + c1 = openmc.Cell(fill=copper, region=-disc) + c2 = openmc.Cell(fill=air, region=+disc & box) + model = openmc.Model() + model.geometry = openmc.Geometry([c1, c2]) + model.settings.run_mode = 'fixed source' + model.settings.particles = 1000 + model.settings.batches = 5 + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point(), + angle=openmc.stats.Monodirectional((0., 0., 1.)) + ) + + # Run model to ensure it doesn't segfault + model.run()