diff --git a/.gitignore b/.gitignore index 136491a4b8..3260cb4946 100644 --- a/.gitignore +++ b/.gitignore @@ -24,6 +24,7 @@ examples/python/**/*.xml # Documentation builds docs/build docs/source/_images/*.pdf +docs/source/_images/*.aux # Source build build diff --git a/docs/Makefile b/docs/Makefile index 89c71dc074..e84aadd33c 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -17,6 +17,8 @@ ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) sou SVG2PDF = inkscape PDFS = $(patsubst %.svg,%.pdf,$(wildcard $(IMAGEDIR)/*.svg)) +# Tikz to PNG conversion +PNGS = $(patsubst %.tex,%.png,$(wildcard $(IMAGEDIR)/*.tex)) .PHONY: help images clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest @@ -43,8 +45,13 @@ help: %.pdf: %.svg $(SVG2PDF) -f $< -A $@ +%.png: %.tex + pdflatex --interaction=nonstopmode --output-directory=$(IMAGEDIR) $< + pdftoppm -r 120 -singlefile $(patsubst %.tex,%.pdf, $<) $(basename $<) + convert -trim -fuzz 2% -transparent white $(patsubst %.tex,%.ppm,$<) $@ + # Rule to build PDFs -images: $(PDFS) +images: $(PDFS) $(PNGS) clean: -rm -rf $(BUILDDIR)/* diff --git a/docs/requirements-rtd.txt b/docs/requirements-rtd.txt new file mode 100644 index 0000000000..652f880ff8 --- /dev/null +++ b/docs/requirements-rtd.txt @@ -0,0 +1,2 @@ +sphinx-numfig +jupyter diff --git a/docs/source/_images/cmfd_flow.png b/docs/source/_images/cmfd_flow.png new file mode 100644 index 0000000000..078f1a4bf8 Binary files /dev/null and b/docs/source/_images/cmfd_flow.png differ diff --git a/docs/source/methods/cmfd_tikz/cmfd_flow.tikz b/docs/source/_images/cmfd_flow.tex similarity index 76% rename from docs/source/methods/cmfd_tikz/cmfd_flow.tikz rename to docs/source/_images/cmfd_flow.tex index 0611918052..4f673b0d77 100644 --- a/docs/source/methods/cmfd_tikz/cmfd_flow.tikz +++ b/docs/source/_images/cmfd_flow.tex @@ -1,3 +1,12 @@ +\documentclass{standalone} +\usepackage[utf8]{inputenc} +\usepackage{amsmath} +\usepackage{tikz} +\usepackage{pgfplots} +\pgfplotsset{compat=1.11} +\usetikzlibrary{shapes,snakes,shadows,arrows,calc,decorations.markings,patterns,fit,matrix,spy} +\pagestyle{empty} +\begin{document} \begin{tikzpicture} \matrix[every node/.style={draw, thick, minimum width=3cm, minimum height=1cm, align=center}, column sep=2cm, row sep=1cm] (m) { \node[draw, fill=red!40] (start) {Batch $i$ \\ tally NDA}; & \\ @@ -16,4 +25,5 @@ \draw (modify.north) -- (end.south); \end{scope} -\end{tikzpicture} \ No newline at end of file +\end{tikzpicture} +\end{document} diff --git a/docs/source/_images/meshfig.png b/docs/source/_images/meshfig.png new file mode 100644 index 0000000000..c84becf3a0 Binary files /dev/null and b/docs/source/_images/meshfig.png differ diff --git a/docs/source/methods/cmfd_tikz/meshfig.tikz b/docs/source/_images/meshfig.tex similarity index 99% rename from docs/source/methods/cmfd_tikz/meshfig.tikz rename to docs/source/_images/meshfig.tex index 1e5bbfa476..d710107cde 100644 --- a/docs/source/methods/cmfd_tikz/meshfig.tikz +++ b/docs/source/_images/meshfig.tex @@ -1,10 +1,20 @@ +\documentclass[tikz]{standalone} +\usepackage[utf8]{inputenc} +\usepackage{amsmath} +\usepackage{tikz} +\usepackage{pgfplots} +\pgfplotsset{compat=1.11} +\usetikzlibrary{shapes,snakes,shadows,arrows,calc,decorations.markings,patterns,fit,matrix,spy} +\usepackage{fixltx2e} +\pagestyle{empty} +\begin{document} % these dimensions are determined in arrow_dimms.ods \def\scale{1.0} \def\latWidth{0.2808363589*\scale} - + \def\RPVOR{3*\scale} \def\rectW{0.75*\scale} \def\RPVIR{2.8694005485*\scale} @@ -20,7 +30,7 @@ \def\bafMIRy{1.5445999739*\scale} \def\bafMORx{1.8544620609*\scale} \def\bafMORy{1.573625702*\scale} - + \tikzset{Assembly/.style={ inner sep=0pt, text width=\latWidth in, @@ -29,13 +39,13 @@ align=center } } - + \def\tkzRPV{(0,0) circle (\RPVIR) (0,0) circle (\RPVOR)} \def\tkzBarrel{(0,0) circle (\BarrelIR) (0,0) circle (\BarrelOR)} \def\tkzShields{(0,0) circle (\BarrelOR) (0,0) circle (\ShieldOR)} - + \def\tkzBaffCOR{(-\bafCORx, -\bafCORy) rectangle (\bafCORx, \bafCORy)} - \def\tkzBaffCIR{(-\bafCIRx, -\bafCIRy) rectangle (\bafCIRx, \bafCIRy)} + \def\tkzBaffCIR{(-\bafCIRx, -\bafCIRy) rectangle (\bafCIRx, \bafCIRy)} \def\tkzBaffMOR{(-\bafMORx, -\bafMORy) rectangle (\bafMORx, \bafMORy)} \def\tkzBaffMIR{(-\bafMIRx, -\bafMIRy) rectangle (\bafMIRx, \bafMIRy) } \def\tkzBaffleC{ \tkzBaffCIR \tkzBaffCOR } @@ -53,7 +63,7 @@ \begin{tikzpicture}[x=1in,y=1in, xshift=3in] \scalebox{0.6}{ % draw RPV, barrel, and shield panels - + \path[fill=black,even odd rule] \tkzRPV; \path[fill=black,even odd rule] \tkzBarrel; \begin{scope} @@ -61,9 +71,9 @@ \path[fill=black,even odd rule] \tkzShields; \end{scope} - + % draw assembly row/column headers - + \draw[red, thick] ($(-7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {R} -- ($(-7*\latWidth,4*\latWidth)$); \draw[red, thick] ($(-6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {P} -- ($(-6*\latWidth,6*\latWidth)$); \draw[red, thick] ($(-5*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {N} -- ($(-5*\latWidth,7*\latWidth)$); @@ -79,7 +89,7 @@ \draw[red, thick] ($(5*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {C} -- ($(5*\latWidth,7*\latWidth)$); \draw[red, thick] ($(6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {B} -- ($(6*\latWidth,6*\latWidth)$); \draw[red, thick] ($(7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {A} -- ($(7*\latWidth,4*\latWidth)$); - + \begin{scope}[rotate=90] \draw[red, thick] ($(-7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {15} -- ($(-7*\latWidth,4*\latWidth)$); \draw[red, thick] ($(-6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {14} -- ($(-6*\latWidth,6*\latWidth)$); @@ -97,7 +107,7 @@ \draw[red, thick] ($(6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {2} -- ($(6*\latWidth,6*\latWidth)$); \draw[red, thick] ($(7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {1} -- ($(7*\latWidth,4*\latWidth)$); \end{scope} - + % draw fuel assembly nodes \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,8*\latWidth)$) {}; \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,8*\latWidth)$) {}; @@ -116,7 +126,7 @@ \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 6*\latWidth,8*\latWidth)$) {}; \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,8*\latWidth)$) {}; \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,8*\latWidth)$) {}; - + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,7*\latWidth)$) {}; \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,7*\latWidth)$) {}; \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-6*\latWidth,7*\latWidth)$) {}; @@ -581,7 +591,7 @@ \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 6*\latWidth,-7*\latWidth)$) {}; \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,-7*\latWidth)$) {}; \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-7*\latWidth)$) {}; - + \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,-8*\latWidth)$) {}; \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,-8*\latWidth)$) {}; \node [Assembly, fill=\lightgray, opacity=0.3] at ($(-6*\latWidth,-8*\latWidth)$) {}; @@ -601,7 +611,7 @@ \node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-8*\latWidth)$) {}; % draw baffle north/south - + \begin{scope}[even odd rule] \clip[rotate=90] \tkzBaffMClip; \path[fill=black] \tkzBaffleC; @@ -611,9 +621,9 @@ \clip \tkzBaffMClip; \path[fill=black, rotate=90] \tkzBaffleM; \end{scope} - + % draw baffle east/west - + \begin{scope}[rotate=90] \begin{scope}[even odd rule] \clip[rotate=90] \tkzBaffMClip; @@ -626,3 +636,4 @@ \end{scope} \end{scope}} \end{tikzpicture} +\end{document} diff --git a/docs/source/_images/openmc200px.png b/docs/source/_images/openmc200px.png new file mode 100644 index 0000000000..3997c6baa0 Binary files /dev/null and b/docs/source/_images/openmc200px.png differ diff --git a/docs/source/conf.py b/docs/source/conf.py index 118ff2c03b..32154fa60f 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -13,6 +13,26 @@ import sys, os +# Determine if we're on Read the Docs server +on_rtd = os.environ.get('READTHEDOCS', None) == 'True' + +# On Read the Docs, we need to mock a few third-party modules so we don't get +# ImportErrors when building documentation +try: + from unittest.mock import MagicMock +except ImportError: + from mock import Mock as MagicMock + + +class Mock(MagicMock): + @classmethod + def __getattr__(cls, name): + return Mock() + +MOCK_MODULES = ['numpy', 'h5py', 'pandas', 'opencg'] +sys.modules.update((mod_name, Mock()) for mod_name in MOCK_MODULES) + + # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. @@ -26,9 +46,8 @@ sys.path.insert(0, os.path.abspath('../..')) # coming with Sphinx (named 'sphinx.ext.*') or your custom ones. extensions = ['sphinx.ext.autodoc', 'sphinx.ext.napoleon', - 'sphinx.ext.pngmath', + 'sphinx.ext.mathjax', 'sphinx.ext.autosummary', - 'sphinxcontrib.tikz', 'sphinx_numfig', 'notebook_sphinxext'] @@ -105,16 +124,15 @@ pygments_style = 'tango' # The theme to use for HTML and HTML Help pages. Major themes that come with # Sphinx are currently 'default' and 'sphinxdoc'. -html_theme = 'haiku' -#html_theme = 'altered_nature' -#html_theme = 'sphinxdoc' - -# Theme options are theme-specific and customize the look and feel of a theme -# further. For a list of options available for each theme, see the -# documentation. -html_theme_options = {'full_logo': True, - 'linkcolor': '#0c3762', - 'visitedlinkcolor': '#0c3762'} +if on_rtd: + html_theme = 'default' + html_logo = '_images/openmc200px.png' +else: + html_theme = 'haiku' + html_theme_options = {'full_logo': True, + 'linkcolor': '#0c3762', + 'visitedlinkcolor': '#0c3762'} + html_logo = '_images/openmc.png' # Add any paths that contain custom themes here, relative to this directory. #html_theme_path = ["_theme"] @@ -126,10 +144,6 @@ html_title = "OpenMC Documentation" # A shorter title for the navigation bar. Default is the same as html_title. #html_short_title = None -# The name of an image file (relative to this directory) to place at the top -# of the sidebar. -html_logo = '_images/openmc.png' - # The name of an image file (within the static path) to use as favicon of the # docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 # pixels large. diff --git a/docs/source/devguide/docbuild.rst b/docs/source/devguide/docbuild.rst index 4c62ca179f..079a69a2c0 100644 --- a/docs/source/devguide/docbuild.rst +++ b/docs/source/devguide/docbuild.rst @@ -12,15 +12,8 @@ is via pip: sudo pip install sphinx -Additionally, you will also need two Sphinx extensions for TikZ support and -numbering figures. The sphinxcontrib-tikz_ package should be installed directly -from the git repository as such: - -.. code-block:: sh - - sudo pip install https://bitbucket.org/philexander/tikz/get/HEAD.tar.gz - -The Numfig_ package can be installed directly with pip: +Additionally, you will also need a Sphinx extension for numbering figures. The +Numfig_ package can be installed directly with pip: .. code-block:: sh diff --git a/docs/source/methods/cmfd.rst b/docs/source/methods/cmfd.rst index 49cb86c674..4181f1b37e 100644 --- a/docs/source/methods/cmfd.rst +++ b/docs/source/methods/cmfd.rst @@ -110,10 +110,12 @@ illustrated as a flow chart below. After a batch of neutrons is simulated, NDA can take place. Each of the steps described above is described in detail in the following sections. -.. tikz:: Flow chart of NDA process. Note "XS" is used for cross section and - "DC" is used for diffusion coefficient. - :libs: shapes, snakes, shadows, arrows, calc, decorations.markings, patterns, fit, matrix, spy - :include: cmfd_tikz/cmfd_flow.tikz +.. figure:: ../_images/cmfd_flow.png + :align: center + :figclass: align-center + + Flow chart of NDA process. Note "XS" is used for cross section and "DC" is + used for diffusion coefficient. Calculation of Macroscopic Cross Sections ----------------------------------------- @@ -422,9 +424,11 @@ during the MC simulation with incoming and outgoing partial currents. This allows the user to not have to worry about neutrons producing adequate tallies in mesh cells far away from the core. -.. tikz:: Diagram of CMFD acceleration mesh - :libs: shapes, snakes, shadows, arrows, calc, decorations.markings, patterns, fit, matrix, spy - :include: cmfd_tikz/meshfig.tikz +.. figure:: ../_images/meshfig.png + :align: center + :figclass: align-center + + Diagram of CMFD acceleration mesh During an MC simulation, CMFD tallies are accumulated. The basic tallies needed are listed in Table :ref:`tab_tally`. Each tally is performed on a spatial and diff --git a/docs/source/methods/physics.rst b/docs/source/methods/physics.rst index e25057488c..cc486ab67b 100644 --- a/docs/source/methods/physics.rst +++ b/docs/source/methods/physics.rst @@ -1626,7 +1626,7 @@ another. .. _ENDF-6 Format: http://www-nds.iaea.org/ndspub/documents/endf/endf102/endf102.pdf -.. _Monte Carlo Sampler: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721_3rdmcsampler.pdf +.. _Monte Carlo Sampler: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721.pdf .. _LA-UR-14-27694: http://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-UR-14-27694 @@ -1636,4 +1636,4 @@ another. .. _lectures: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-05-4983.pdf -.. _MCNP Manual: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/MCNP5_Manual_Volume_I_LA-UR-03-1987.pdf +.. _MCNP Manual: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-03-1987.pdf diff --git a/docs/source/methods/random_numbers.rst b/docs/source/methods/random_numbers.rst index 3ea61719c5..e0a974f43b 100644 --- a/docs/source/methods/random_numbers.rst +++ b/docs/source/methods/random_numbers.rst @@ -70,5 +70,5 @@ the idea is to determine the new multiplicative and additive constants in Different Sizes and Good Lattice Structures," *Math. Comput.*, **68**, 249 (1999). -.. _Brown: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/anl_rn_arb-strides_1994.pdf +.. _Brown: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/anl-rn-arb-stride.pdf .. _linear congruential generator: http://en.wikipedia.org/wiki/Linear_congruential_generator diff --git a/docs/source/methods/tallies.rst b/docs/source/methods/tallies.rst index fd0812245a..65a0989ac7 100644 --- a/docs/source/methods/tallies.rst +++ b/docs/source/methods/tallies.rst @@ -507,6 +507,6 @@ improve the estimate of the percentile. .. _Cauchy distribution: http://en.wikipedia.org/wiki/Cauchy_distribution -.. _unpublished rational approximation: http://home.online.no/~pjacklam/notes/invnorm/ +.. _unpublished rational approximation: https://web.archive.org/web/20150926021742/http://home.online.no/~pjacklam/notes/invnorm/ .. _MC21: http://www.osti.gov/bridge/servlets/purl/903083-HT5p1o/903083.pdf diff --git a/docs/source/usersguide/output/statepoint.rst b/docs/source/usersguide/output/statepoint.rst index d3c1729af6..15bc79f739 100644 --- a/docs/source/usersguide/output/statepoint.rst +++ b/docs/source/usersguide/output/statepoint.rst @@ -185,10 +185,6 @@ if run_mode == 'k-eigenvalue': Type of the j-th filter. Can be 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', or 'distribcell'. -**/tallies/tally /filter /offset** (*int*) - - Filter offset (used for distribcell filter). - **/tallies/tally /filter /n_bins** (*int*) Number of bins for the j-th filter. diff --git a/docs/source/usersguide/output/summary.rst b/docs/source/usersguide/output/summary.rst index f87f60c4a0..66a5cad154 100644 --- a/docs/source/usersguide/output/summary.rst +++ b/docs/source/usersguide/output/summary.rst @@ -121,6 +121,10 @@ The current revision of the summary file format is 1. Region specification for the cell. +**/geometry/cells/cell /distribcell_index** (*int*) + + Index of this cell in distribcell filter arrays. + **/geometry/surfaces/surface /index** (*int*) Index in surfaces array used internally in OpenMC. @@ -306,6 +310,11 @@ The current revision of the summary file format is 1. than the number of user-specified scores since each score might have multiple scoring bins, e.g., scatter-PN. +**/tallies/tally /moment_orders** (*char[][]*) + + Tallying moment orders for Legendre and spherical harmonic tally expansions + (*e.g.*, 'P2', 'Y1,2', etc.). + **/tallies/tally /score_bins** (*char[][]*) Scoring bins for the tally. diff --git a/docs/sphinxext/notebook_sphinxext.py b/docs/sphinxext/notebook_sphinxext.py index a9d634a351..eab93cc8a5 100644 --- a/docs/sphinxext/notebook_sphinxext.py +++ b/docs/sphinxext/notebook_sphinxext.py @@ -12,10 +12,7 @@ from docutils.parsers.rst import directives, roles, states from docutils.parsers.rst.roles import set_classes from docutils.transforms import misc -try: - from IPython.nbconver.exporters import html -except ImportError: - from IPython.nbconvert import html +from nbconvert import html class Notebook(Directive): diff --git a/openmc/filter.py b/openmc/filter.py index 04935b8edc..d61fb648f6 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -42,8 +42,6 @@ class Filter(object): The number of filter bins mesh : Mesh or None A Mesh object for 'mesh' type filters. - offset : Integral - A value used to index tally bins for 'distribcell' tallies. stride : Integral The number of filter, nuclide and score bins within each of this filter's bins. @@ -57,7 +55,6 @@ class Filter(object): self._num_bins = 0 self._bins = None self._mesh = None - self._offset = -1 self._stride = None if type is not None: @@ -93,7 +90,6 @@ class Filter(object): clone._bins = copy.deepcopy(self.bins, memo) clone._num_bins = self.num_bins clone._mesh = copy.deepcopy(self.mesh, memo) - clone._offset = self.offset clone._stride = self.stride memo[id(self)] = clone @@ -108,7 +104,6 @@ class Filter(object): string = 'Filter\n' string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self.type) string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self.bins) - string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self.offset) return string @property @@ -134,10 +129,6 @@ class Filter(object): def mesh(self): return self._mesh - @property - def offset(self): - return self._offset - @property def stride(self): return self._stride @@ -226,11 +217,6 @@ class Filter(object): self.type = 'mesh' self.bins = self.mesh.id - @offset.setter - def offset(self, offset): - cv.check_type('filter offset', offset, Integral) - self._offset = offset - @stride.setter def stride(self, stride): cv.check_type('filter stride', stride, Integral) @@ -623,7 +609,7 @@ class Filter(object): # If this region is in Cell corresponding to the # distribcell filter bin, store it in dictionary if cell_id == self.bins[0]: - offset = openmc_geometry.get_offset(path, self.offset) + offset = openmc_geometry.get_cell_instance(path) offsets_to_coords[offset] = coords # Each distribcell offset is a DataFrame bin diff --git a/openmc/geometry.py b/openmc/geometry.py index e848e0cddf..fc8e19f071 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -42,10 +42,10 @@ class Geometry(object): self._root_universe = root_universe - def get_offset(self, path, filter_offset): - """Returns the corresponding location in the results array for a given path and - filter number. This is primarily intended to post-processing result when - a distribcell filter is used. + def get_cell_instance(self, path): + """Return the instance number for the final cell in a geometry path. + + The instance is an index into tally distribcell filter arrays. Parameters ---------- @@ -55,24 +55,31 @@ class Geometry(object): lattice passed through. For the case of the lattice, a tuple should be provided to indicate which coordinates in the lattice should be entered. This should be in the form: (lat_id, i_x, i_y, i_z) - filter_offset : int - An integer that specifies which offset map the filter is using Returns ------- - offset : int - Location in the results array for the path and filter + instance : int + Index in tally results array for distribcell filters """ + # Find the distribcell index of the cell. + cells = self.get_all_cells() + if path[-1] in cells: + distribcell_index = cells[path[-1]].distribcell_index + else: + raise RuntimeError('Could not find cell {} specified in a \ + distribcell filter'.format(path[-1])) + # Return memoize'd offset if possible - if (path, filter_offset) in self._offsets: - offset = self._offsets[(path, filter_offset)] + if (path, distribcell_index) in self._offsets: + offset = self._offsets[(path, distribcell_index)] # Begin recursive call to compute offset starting with the base Universe else: - offset = self._root_universe.get_offset(path, filter_offset) - self._offsets[(path, filter_offset)] = offset + offset = self._root_universe.get_cell_instance(path, + distribcell_index) + self._offsets[(path, distribcell_index)] = offset # Return the final offset return offset diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index ef3e90fc41..ca9fa5d669 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -73,6 +73,9 @@ class Library(object): name : str, optional Name of the multi-group cross section library. Used as a label to identify tallies in OpenMC 'tallies.xml' file. + sparse : bool + Whether or not the Library's tallies use SciPy's LIL sparse matrix + format for compressed data storage """ @@ -92,6 +95,7 @@ class Library(object): self._all_mgxs = OrderedDict() self._sp_filename = None self._keff = None + self._sparse = False self.name = name self.openmc_geometry = openmc_geometry @@ -119,6 +123,7 @@ class Library(object): clone._all_mgxs = self.all_mgxs clone._sp_filename = self._sp_filename clone._keff = self._keff + clone._sparse = self.sparse clone._all_mgxs = OrderedDict() for domain in self.domains: @@ -208,6 +213,10 @@ class Library(object): def keff(self): return self._keff + @property + def sparse(self): + return self._sparse + @openmc_geometry.setter def openmc_geometry(self, openmc_geometry): cv.check_type('openmc_geometry', openmc_geometry, openmc.Geometry) @@ -286,6 +295,28 @@ class Library(object): cv.check_type('tally trigger', tally_trigger, openmc.Trigger) self._tally_trigger = tally_trigger + @sparse.setter + def sparse(self, sparse): + """Convert tally data from NumPy arrays to SciPy list of lists (LIL) + sparse matrices, and vice versa. + + This property may be used to reduce the amount of data in memory during + tally data processing. The tally data will be stored as SciPy LIL + matrices internally within the Tally object. All tally data access + properties and methods will return data as a dense NumPy array. + + """ + + cv.check_type('sparse', sparse, bool) + + # Sparsify or densify each MGXS in the Library + for domain in self.domains: + for mgxs_type in self.mgxs_types: + mgxs = self.get_mgxs(domain, mgxs_type) + mgxs.sparse = self.sparse + + self._sparse = sparse + def build_library(self): """Initialize MGXS objects in each domain and for each reaction type in the library. @@ -381,6 +412,7 @@ class Library(object): for mgxs_type in self.mgxs_types: mgxs = self.get_mgxs(domain, mgxs_type) mgxs.load_from_statepoint(statepoint) + mgxs.sparse = self.sparse def get_mgxs(self, domain, mgxs_type): """Return the MGXS object for some domain and reaction rate type. diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 635c822e31..968504a06d 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -47,7 +47,7 @@ _DOMAINS = [openmc.Cell, class MGXS(object): - """An abstract multi-group cross section for some energy group structure + """An abstract multi-group cross section for some energy group structure within some spatial domain. This class can be used for both OpenMC input generation and tally data @@ -92,10 +92,10 @@ class MGXS(object): xs_tally : Tally Derived tally for the multi-group cross section. This attribute is None unless the multi-group cross section has been computed. - num_sumbdomains : Integral + num_subdomains : Integral The number of subdomains is unity for 'material', 'cell' and 'universe' domain types. When the This is equal to the number of cell instances - for 'distribcell' domain types (it is equal to unity prior to loading + for 'distribcell' domain types (it is equal to unity prior to loading tally data from a statepoint file). num_nuclides : Integral The number of nuclides for which the multi-group cross section is @@ -103,6 +103,9 @@ class MGXS(object): nuclides : list of str or 'sum' A list of nuclide string names (e.g., 'U-238', 'O-16') when by_nuclide is True and 'sum' when by_nuclide is False. + sparse : bool + Whether or not the MGXS' tallies use SciPy's LIL sparse matrix format + for compressed data storage """ @@ -121,6 +124,7 @@ class MGXS(object): self._tally_trigger = None self._tallies = None self._xs_tally = None + self._sparse = False self.name = name self.by_nuclide = by_nuclide @@ -146,6 +150,7 @@ class MGXS(object): clone._energy_groups = copy.deepcopy(self.energy_groups, memo) clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo) clone._xs_tally = copy.deepcopy(self.xs_tally, memo) + clone._sparse = self.sparse clone._tallies = OrderedDict() for tally_type, tally in self.tallies.items(): @@ -199,6 +204,10 @@ class MGXS(object): def xs_tally(self): return self._xs_tally + @property + def sparse(self): + return self._sparse + @property def num_subdomains(self): tally = list(self.tallies.values())[0] @@ -249,6 +258,29 @@ class MGXS(object): cv.check_type('tally trigger', tally_trigger, openmc.Trigger) self._tally_trigger = tally_trigger + @sparse.setter + def sparse(self, sparse): + """Convert tally data from NumPy arrays to SciPy list of lists (LIL) + sparse matrices, and vice versa. + + This property may be used to reduce the amount of data in memory during + tally data processing. The tally data will be stored as SciPy LIL + matrices internally within the Tally object. All tally data access + properties and methods will return data as a dense NumPy array. + + """ + + cv.check_type('sparse', sparse, bool) + + # Sparsify or densify the derived MGXS tally and its base tallies + if self.xs_tally: + self.xs_tally.sparse = sparse + + for tally_name in self.tallies: + self.tallies[tally_name].sparse = sparse + + self._sparse = sparse + @staticmethod def get_mgxs(mgxs_type, domain=None, domain_type=None, energy_groups=None, by_nuclide=False, name=''): @@ -503,6 +535,7 @@ class MGXS(object): # Remove NaNs which may have resulted from divide-by-zero operations self._xs_tally._mean = np.nan_to_num(self.xs_tally.mean) self._xs_tally._std_dev = np.nan_to_num(self.xs_tally.std_dev) + self.xs_tally.sparse = self.sparse def load_from_statepoint(self, statepoint): """Extracts tallies in an OpenMC StatePoint with the data needed to @@ -565,6 +598,7 @@ class MGXS(object): estimator=tally.estimator) sp_tally = sp_tally.get_slice(tally.scores, filters, filter_bins, tally.nuclides) + sp_tally.sparse = self.sparse self.tallies[tally_type] = sp_tally # Compute the cross section from the tallies @@ -716,6 +750,7 @@ class MGXS(object): # Clone this MGXS to initialize the condensed version condensed_xs = copy.deepcopy(self) + condensed_xs.sparse = False condensed_xs.energy_groups = coarse_groups # Build energy indices to sum across @@ -754,10 +789,8 @@ class MGXS(object): std_dev = np.sqrt(std_dev) # Reshape condensed data arrays with one dimension for all filters - new_shape = \ - (tally.num_filter_bins, tally.num_nuclides, tally.num_score_bins,) - mean = np.reshape(mean, new_shape) - std_dev = np.reshape(std_dev, new_shape) + mean = np.reshape(mean, tally.shape) + std_dev = np.reshape(std_dev, tally.shape) # Override tally's data with the new condensed data tally._mean = mean @@ -765,6 +798,7 @@ class MGXS(object): # Compute the energy condensed multi-group cross section condensed_xs.compute_xs() + condensed_xs.sparse = self.sparse return condensed_xs def get_subdomain_avg_xs(self, subdomains='all'): @@ -807,8 +841,9 @@ class MGXS(object): # Clone this MGXS to initialize the subdomain-averaged version avg_xs = copy.deepcopy(self) + avg_xs.sparse = False - # If domain is distribcell, make the new domain 'cell' + # If domain is distribcell, make the new domain 'cell' if self.domain_type == 'distribcell': avg_xs.domain_type = 'cell' @@ -836,10 +871,8 @@ class MGXS(object): domain_filter.num_bins = 1 # Reshape averaged data arrays with one dimension for all filters - new_shape = \ - (tally.num_filter_bins, tally.num_nuclides, tally.num_score_bins,) - mean = np.reshape(mean, new_shape) - std_dev = np.reshape(std_dev, new_shape) + mean = np.reshape(mean, tally.shape) + std_dev = np.reshape(std_dev, tally.shape) # Override tally's data with the new condensed data tally._mean = mean @@ -847,7 +880,7 @@ class MGXS(object): # Compute the subdomain-averaged multi-group cross section avg_xs.compute_xs() - + avg_xs.sparse = self.sparse return avg_xs def print_xs(self, subdomains='all', nuclides='all', xs_type='macro'): @@ -1420,8 +1453,8 @@ class AbsorptionXS(MGXS): class CaptureXS(MGXS): """A capture multi-group cross section. - - The Neutron capture reaction rate is defined as the difference between + + The neutron capture reaction rate is defined as the difference between OpenMC's 'absorption' and 'fission' reaction rate score types. This includes not only radiative capture, but all forms of neutron disappearance aside from fission (e.g., MT > 100). @@ -1723,8 +1756,8 @@ class ScatterMatrixXS(MGXS): if self.correction == 'P0': scatter_p1 = self.tallies['scatter-P1'] scatter_p1 = scatter_p1.get_slice(scores=['scatter-P1']) - energy_filter = openmc.Filter(type='energy') - energy_filter.bins = self.energy_groups.group_edges + energy_filter = self.tallies['scatter'].find_filter('energy') + energy_filter = copy.deepcopy(energy_filter) scatter_p1 = scatter_p1.diagonalize_filter(energy_filter) rxn_tally = self.tallies['scatter'] - scatter_p1 else: diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index afa57c78d9..0bda48c160 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -725,11 +725,11 @@ def get_openmc_cell(opencg_cell): else: openmc_cell.fill = get_openmc_material(fill) - if opencg_cell.rotation: + if opencg_cell.rotation is not None: rotation = np.asarray(opencg_cell.rotation, dtype=np.float64) openmc_cell.rotation = rotation - if opencg_cell.translation: + if opencg_cell.translation is not None: translation = np.asarray(opencg_cell.translation, dtype=np.float64) openmc_cell.translation = translation @@ -881,16 +881,30 @@ def get_opencg_lattice(openmc_lattice): universes = openmc_lattice.universes outer = openmc_lattice.outer + # Convert 2D dimension to 3D for OpenCG + if len(dimension) == 2: + new_dimension = np.ones(3, dtype=np.int) + new_dimension[:2] = dimension + dimension = new_dimension + + # Convert 2D pitch to 3D for OpenCG if len(pitch) == 2: - new_pitch = np.ones(3, dtype=np.float64) * np.inf + new_pitch = np.ones(3, dtype=np.float64) * np.finfo(np.float64).max new_pitch[:2] = pitch pitch = new_pitch + # Convert 2D lower left to 3D for OpenCG if len(lower_left) == 2: - new_lower_left = np.ones(3, dtype=np.float64) + new_lower_left = np.ones(3, dtype=np.float64) * np.finfo(np.float64).min new_lower_left[:2] = lower_left lower_left = new_lower_left + # Convert 2D universes array to 3D for OpenCG + if len(universes.shape) == 2: + new_universes = universes.copy() + new_universes.shape = (1,) + universes.shape + universes = new_universes + # Initialize an empty array for the OpenCG nested Universes in this Lattice universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), dtype=opencg.Universe) @@ -905,7 +919,7 @@ def get_opencg_lattice(openmc_lattice): for z in range(dimension[2]): for y in range(dimension[1]): for x in range(dimension[0]): - universe_id = universes[x][dimension[1]-y-1][z].id + universe_id = universes[z][y][x].id universe_array[z][y][x] = unique_universes[universe_id] opencg_lattice = opencg.Lattice(lattice_id, name) @@ -963,7 +977,7 @@ def get_openmc_lattice(opencg_lattice): outer = opencg_lattice.outside # Initialize an empty array for the OpenMC nested Universes in this Lattice - universe_array = np.ndarray(tuple(np.array(dimension)), + universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), dtype=openmc.Universe) # Create OpenMC Universes for each unique nested Universe in this Lattice @@ -977,7 +991,7 @@ def get_openmc_lattice(opencg_lattice): for y in range(dimension[1]): for x in range(dimension[0]): universe_id = universes[z][y][x].id - universe_array[x][y][z] = unique_universes[universe_id] + universe_array[z][y][x] = unique_universes[universe_id] # Reverse y-dimension in array to match ordering in OpenCG universe_array = universe_array[:, ::-1, :] diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 2edf9badd4..3c0759f124 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -9,9 +9,9 @@ if sys.version > '3': class StatePoint(object): - """State information on a simulation at a certain point in time (at the end of a - given batch). Statepoints can be used to analyze tally results as well as - restart a simulation. + """State information on a simulation at a certain point in time (at the end + of a given batch). Statepoints can be used to analyze tally results as well + as restart a simulation. Attributes ---------- @@ -75,6 +75,9 @@ class StatePoint(object): energy of the source site. source_present : bool Indicate whether source sites are present + sparse : bool + Whether or not the tallies uses SciPy's LIL sparse matrix format for + compressed data storage tallies : dict Dictionary whose keys are tally IDs and whose values are Tally objects tallies_present : bool @@ -110,6 +113,7 @@ class StatePoint(object): self._tallies_read = False self._summary = False self._global_tallies = None + self._sparse = False def close(self): self._f.close() @@ -318,6 +322,10 @@ class StatePoint(object): def source_present(self): return self._f['source_present'].value > 0 + @property + def sparse(self): + return self._sparse + @property def tallies(self): if not self._tallies_read: @@ -340,7 +348,8 @@ class StatePoint(object): for tally_key in tally_keys: # Read the Tally size specifications - n_realizations = self._f['{0}{1}/n_realizations'.format(base, tally_key)].value + n_realizations = \ + self._f['{0}{1}/n_realizations'.format(base, tally_key)].value # Create Tally object and assign basic properties tally = openmc.Tally(tally_id=tally_key) @@ -350,7 +359,8 @@ class StatePoint(object): tally.num_realizations = n_realizations # Read the number of Filters - n_filters = self._f['{0}{1}/n_filters'.format(base, tally_key)].value + n_filters = \ + self._f['{0}{1}/n_filters'.format(base, tally_key)].value subbase = '{0}{1}/filter '.format(base, tally_key) @@ -358,10 +368,8 @@ class StatePoint(object): for j in range(1, n_filters+1): # Read the Filter type - filter_type = self._f['{0}{1}/type'.format(subbase, j)].value.decode() - - # Read the Filter offset - offset = self._f['{0}{1}/offset'.format(subbase, j)].value + filter_type = \ + self._f['{0}{1}/type'.format(subbase, j)].value.decode() n_bins = self._f['{0}{1}/n_bins'.format(subbase, j)].value @@ -370,7 +378,6 @@ class StatePoint(object): # Create Filter object filter = openmc.Filter(filter_type, bins) - filter.offset = offset filter.num_bins = n_bins if filter_type == 'mesh': @@ -384,22 +391,18 @@ class StatePoint(object): tally.add_filter(filter) # Read Nuclide bins - nuclide_names = self._f['{0}{1}/nuclides'.format(base, tally_key)].value + nuclide_names = \ + self._f['{0}{1}/nuclides'.format(base, tally_key)].value # Add all Nuclides to the Tally for name in nuclide_names: nuclide = openmc.Nuclide(name.decode().strip()) tally.add_nuclide(nuclide) - # Read score bins - n_score_bins = self._f['{0}{1}/n_score_bins'.format(base, tally_key)].value - - tally.num_score_bins = n_score_bins - scores = self._f['{0}{1}/score_bins'.format( base, tally_key)].value - n_user_scores = self._f['{0}{1}/n_user_score_bins' - .format(base, tally_key)].value + n_score_bins = self._f['{0}{1}/n_score_bins' + .format(base, tally_key)].value # Compute and set the filter strides for i in range(n_filters): @@ -409,7 +412,7 @@ class StatePoint(object): for j in range(i+1, n_filters): filter.stride *= tally.filters[j].num_bins - # Read scattering moment order strings (e.g., P3, Y-1,2, etc.) + # Read scattering moment order strings (e.g., P3, Y1,2, etc.) moments = self._f['{0}{1}/moment_orders'.format( base, tally_key)].value @@ -424,6 +427,7 @@ class StatePoint(object): tally.add_score(score) # Add Tally to the global dictionary of all Tallies + tally.sparse = self.sparse self._tallies[tally_key] = tally self._tallies_read = True @@ -448,6 +452,26 @@ class StatePoint(object): def with_summary(self): return False if self.summary is None else True + @sparse.setter + def sparse(self, sparse): + """Convert tally data from NumPy arrays to SciPy list of lists (LIL) + sparse matrices, and vice versa. + + This property may be used to reduce the amount of data in memory during + tally data processing. The tally data will be stored as SciPy LIL + matrices internally within each Tally object. All tally data access + properties and methods will return data as a dense NumPy array. + + """ + + cv.check_type('sparse', sparse, bool) + self._sparse = sparse + + # Update tally sparsities + if self._tallies_read: + for tally_id in self.tallies: + self.tallies[tally_id].sparse = self.sparse + def get_tally(self, scores=[], filters=[], nuclides=[], name=None, id=None, estimator=None): """Finds and returns a Tally object with certain properties. diff --git a/openmc/summary.py b/openmc/summary.py index aff749200f..80f9fc3976 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -1,4 +1,5 @@ import numpy as np +import re import openmc from openmc.region import Region @@ -274,6 +275,11 @@ class Summary(object): cell.region = Region.from_expression( region, {s.id: s for s in self.surfaces.values()}) + # Get the distribcell index + ind = self._f['geometry/cells'][key]['distribcell_index'].value + if ind != 0: + cell.distribcell_index = ind + # Add the Cell to the global dictionary of all Cells self.cells[index] = cell @@ -332,11 +338,8 @@ class Summary(object): self._f['geometry/lattices'][key]['lower_left'][...] pitch = self._f['geometry/lattices'][key]['pitch'][...] outer = self._f['geometry/lattices'][key]['outer'].value - universe_ids = \ - self._f['geometry/lattices'][key]['universes'][...] - universe_ids = np.swapaxes(universe_ids, 0, 1) - universe_ids = np.swapaxes(universe_ids, 1, 2) + self._f['geometry/lattices'][key]['universes'][...] # Create the Lattice lattice = openmc.RectLattice(lattice_id=lattice_id, name=name) @@ -352,22 +355,22 @@ class Summary(object): universes = \ np.ndarray(tuple(universe_ids.shape), dtype=openmc.Universe) - for x in range(universe_ids.shape[0]): + for z in range(universe_ids.shape[0]): for y in range(universe_ids.shape[1]): - for z in range(universe_ids.shape[2]): - universes[x, y, z] = \ - self.get_universe_by_id(universe_ids[x, y, z]) + for x in range(universe_ids.shape[2]): + universes[z, y, x] = \ + self.get_universe_by_id(universe_ids[z, y, x]) - # Transpose, reverse y-dimension for appropriate ordering - shape = universes.shape - universes = np.transpose(universes, (1, 0, 2)) - universes.shape = shape - universes = universes[:, ::-1, :] + # Use 2D NumPy array to store lattice universes for 2D lattices + if len(dimension) == 2: + universes = np.squeeze(universes) + universes = np.atleast_2d(universes) + + # Set the universes for the lattice lattice.universes = universes if offsets is not None: - offsets = np.swapaxes(offsets, 0, 1) - offsets = np.swapaxes(offsets, 1, 2) + offsets = np.swapaxes(offsets, 0, 2) lattice.offsets = offsets # Add the Lattice to the global dictionary of all Lattices @@ -523,12 +526,18 @@ class Summary(object): # Create Tally object and assign basic properties tally = openmc.Tally(tally_id, tally_name) + # Read scattering moment order strings (e.g., P3, Y1,2, etc.) + moments = self._f['{0}/moment_orders'.format(subbase)].value + # Read score metadata scores = self._f['{0}/score_bins'.format(subbase)].value - for score in scores: - tally.add_score(score.decode()) - num_score_bins = self._f['{0}/n_score_bins'.format(subbase)][...] - tally.num_score_bins = num_score_bins + for j, score in enumerate(scores): + score = score.decode() + + # If this is a moment, use generic moment order + pattern = r'-n$|-pn$|-yn$' + score = re.sub(pattern, '-' + moments[j].decode(), score) + tally.add_score(score) # Read filter metadata num_filters = self._f['{0}/n_filters'.format(subbase)].value diff --git a/openmc/tallies.py b/openmc/tallies.py index 13a219dedd..7ec1ca50cb 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -24,6 +24,13 @@ if sys.version_info[0] >= 3: # "Static" variable for auto-generated Tally IDs AUTO_TALLY_ID = 10000 +# The tally arithmetic product types. The tensor product performs the full +# cross product of the data in two tallies with respect to a specified axis +# (filters, nuclides, or scores). The entrywise product performs the arithmetic +# operation entrywise across the entries in two tallies with respect to a +# specified axis. +_PRODUCT_TYPES = ['tensor', 'entrywise'] + def reset_auto_tally_id(): global AUTO_TALLY_ID @@ -58,16 +65,17 @@ class Tally(object): Type of estimator for the tally triggers : list of openmc.trigger.Trigger List of tally triggers - num_score_bins : Integral + num_scores : Integral Total number of scores, accounting for the fact that a single user-specified score, e.g. scatter-P3 or flux-Y2,2, might have multiple bins - num_scores : Integral - Total number of user-specified scores num_filter_bins : Integral Total number of filter bins accounting for all filters num_bins : Integral Total number of bins for the tally + shape : 3-tuple of Integral + The shape of the tally data array ordered as the number of filter bins, + nuclide bins and score bins num_realizations : Integral Total number of realizations with_summary : bool @@ -81,6 +89,11 @@ class Tally(object): An array containing the sample mean for each bin std_dev : ndarray An array containing the sample standard deviation for each bin + derived : bool + Whether or not the tally is derived from one or more other tallies + sparse : bool + Whether or not the tally uses SciPy's LIL sparse matrix format for + compressed data storage """ @@ -94,7 +107,6 @@ class Tally(object): self._estimator = None self._triggers = [] - self._num_score_bins = 0 self._num_realizations = 0 self._with_summary = False @@ -104,6 +116,7 @@ class Tally(object): self._std_dev = None self._with_batch_statistics = False self._derived = False + self._sparse = False self._sp_filename = None self._results_read = False @@ -117,7 +130,6 @@ class Tally(object): clone.id = self.id clone.name = self.name clone.estimator = self.estimator - clone.num_score_bins = self.num_score_bins clone.num_realizations = self.num_realizations clone._sum = copy.deepcopy(self._sum, memo) clone._sum_sq = copy.deepcopy(self._sum_sq, memo) @@ -126,6 +138,7 @@ class Tally(object): clone._with_summary = self.with_summary clone._with_batch_statistics = self.with_batch_statistics clone._derived = self.derived + clone._sparse = self.sparse clone._sp_filename = self._sp_filename clone._results_read = self._results_read @@ -246,10 +259,6 @@ class Tally(object): def num_scores(self): return len(self._scores) - @property - def num_score_bins(self): - return self._num_score_bins - @property def num_filter_bins(self): num_bins = 1 @@ -263,9 +272,13 @@ class Tally(object): def num_bins(self): num_bins = self.num_filter_bins num_bins *= self.num_nuclides - num_bins *= self.num_score_bins + num_bins *= self.num_scores return num_bins + @property + def shape(self): + return (self.num_filter_bins, self.num_nuclides, self.num_scores) + @property def estimator(self): return self._estimator @@ -299,29 +312,33 @@ class Tally(object): sum = data['sum'] sum_sq = data['sum_sq'] - # Define a routine to convert 0 to 1 - def nonzero(val): - return 1 if not val else val - # Reshape the results arrays - new_shape = (nonzero(self.num_filter_bins), - nonzero(self.num_nuclides), - nonzero(self.num_score_bins)) - - sum = np.reshape(sum, new_shape) - sum_sq = np.reshape(sum_sq, new_shape) + sum = np.reshape(sum, self.shape) + sum_sq = np.reshape(sum_sq, self.shape) # Set the data for this Tally self._sum = sum self._sum_sq = sum_sq + # Convert NumPy arrays to SciPy sparse LIL matrices + if self.sparse: + import scipy.sparse as sps + + self._sum = \ + sps.lil_matrix(self._sum.flatten(), self._sum.shape) + self._sum_sq = \ + sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape) + # Indicate that Tally results have been read self._results_read = True # Close the HDF5 statepoint file f.close() - return self._sum + if self.sparse: + return np.reshape(self._sum.toarray(), self.shape) + else: + return self._sum @property def sum_sq(self): @@ -332,7 +349,10 @@ class Tally(object): # Force reading of sum and sum_sq self.sum - return self._sum_sq + if self.sparse: + return np.reshape(self._sum_sq.toarray(), self.shape) + else: + return self._sum_sq @property def mean(self): @@ -341,7 +361,18 @@ class Tally(object): return None self._mean = self.sum / self.num_realizations - return self._mean + + # Convert NumPy array to SciPy sparse LIL matrix + if self.sparse: + import scipy.sparse as sps + + self._mean = \ + sps.lil_matrix(self._mean.flatten(), self._mean.shape) + + if self.sparse: + return np.reshape(self._mean.toarray(), self.shape) + else: + return self._mean @property def std_dev(self): @@ -354,8 +385,20 @@ class Tally(object): self._std_dev = np.zeros_like(self.mean) self._std_dev[nonzero] = np.sqrt((self.sum_sq[nonzero]/n - self.mean[nonzero]**2)/(n - 1)) + + # Convert NumPy array to SciPy sparse LIL matrix + if self.sparse: + import scipy.sparse as sps + + self._std_dev = \ + sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape) + self.with_batch_statistics = True - return self._std_dev + + if self.sparse: + return np.reshape(self._std_dev.toarray(), self.shape) + else: + return self._std_dev @property def with_batch_statistics(self): @@ -365,6 +408,10 @@ class Tally(object): def derived(self): return self._derived + @property + def sparse(self): + return self._sparse + @estimator.setter def estimator(self, estimator): cv.check_value('estimator', estimator, @@ -452,9 +499,13 @@ class Tally(object): 'not a string'.format(score, self.id) raise ValueError(msg) - # If the score is already in the Tally, don't add it again + # If the score is already in the Tally, raise an error if score in self.scores: - return + msg = 'Unable to add a duplicate score {0} to Tally ID="{1}" ' \ + 'since duplicate scores are not supported in the OpenMC ' \ + 'Python API'.format(score, self.id) + raise ValueError(msg) + # Normal score strings if isinstance(score, basestring): self._scores.append(score.strip()) @@ -462,10 +513,6 @@ class Tally(object): else: self._scores.append(score) - @num_score_bins.setter - def num_score_bins(self, num_score_bins): - self._num_score_bins = num_score_bins - @num_realizations.setter def num_realizations(self, num_realizations): cv.check_type('number of realizations', num_realizations, Integral) @@ -492,6 +539,51 @@ class Tally(object): cv.check_type('sum_sq', sum_sq, Iterable) self._sum_sq = sum_sq + @sparse.setter + def sparse(self, sparse): + """Convert tally data from NumPy arrays to SciPy list of lists (LIL) + sparse matrices, and vice versa. + + This property may be used to reduce the amount of data in memory during + tally data processing. The tally data will be stored as SciPy LIL + matrices internally within the Tally object. All tally data access + properties and methods will return data as a dense NumPy array. + + """ + + cv.check_type('sparse', sparse, bool) + + # Convert NumPy arrays to SciPy sparse LIL matrices + if sparse and not self.sparse: + import scipy.sparse as sps + + if self._sum is not None: + self._sum = \ + sps.lil_matrix(self._sum.flatten(), self._sum.shape) + if self._sum_sq is not None: + self._sum_sq = \ + sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape) + if self._mean is not None: + self._mean = \ + sps.lil_matrix(self._mean.flatten(), self._mean.shape) + if self._std_dev is not None: + self._std_dev = \ + sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape) + + self._sparse = True + + # Convert SciPy sparse LIL matrices to NumPy arrays + elif not sparse and self.sparse: + if self._sum is not None: + self._sum = np.reshape(self._sum.toarray(), self.shape) + if self._sum_sq is not None: + self._sum_sq = np.reshape(self._sum_sq.toarray(), self.shape) + if self._mean is not None: + self._mean = np.reshape(self._mean.toarray(), self.shape) + if self._std_dev is not None: + self._std_dev = np.reshape(self._std_dev.toarray(), self.shape) + self._sparse = False + def remove_score(self, score): """Remove a score from the tally @@ -619,7 +711,8 @@ class Tally(object): """ if not self.can_merge(tally): - msg = 'Unable to merge tally ID="{0}" with "{1}"'.format(tally.id, self.id) + msg = 'Unable to merge tally ID="{0}" with ' + \ + '"{1}"'.format(tally.id, self.id) raise ValueError(msg) # Create deep copy of tally to return as merged tally @@ -636,9 +729,10 @@ class Tally(object): merged_tally.filters[i] = merged_filter break - # Add scores from second tally to merged tally + # Add unique scores from second tally to merged tally for score in tally.scores: - merged_tally.add_score(score) + if score not in merged_tally.scores: + merged_tally.add_score(score) # Add triggers from second tally to merged tally for trigger in tally.triggers: @@ -1004,7 +1098,12 @@ class Tally(object): """ - cv.check_iterable_type('scores', scores, basestring) + for score in scores: + if not isinstance(score, (basestring, CrossScore)): + msg = 'Unable to get score indices for score "{0}" in Tally ' \ + 'ID="{1}" since it is not a string or CrossScore'\ + .format(score, self.id) + raise ValueError(msg) # Determine the score indices from any of the requested scores if scores: @@ -1062,22 +1161,19 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. ValueError is also thrown - if the input parameters do not correspond to the Tally's attributes, + When this method is called before the Tally is populated with data, + or the input parameters do not correspond to the Tally's attributes, e.g., if the score(s) do not match those in the Tally. """ - # Ensure that StatePoint.read_results() was called first + # Ensure that the tally has data if (value == 'mean' and self.mean is None) or \ (value == 'std_dev' and self.std_dev is None) or \ (value == 'rel_err' and self.mean is None) or \ (value == 'sum' and self.sum is None) or \ (value == 'sum_sq' and self.sum_sq is None): - msg = 'The Tally ID="{0}" has no data to return. Call the ' \ - 'StatePoint.read_results() method before using ' \ - 'Tally.get_values(...)'.format(self.id) + msg = 'The Tally ID="{0}" has no data to return'.format(self.id) raise ValueError(msg) # Get filter, nuclide and score indices @@ -1144,17 +1240,14 @@ class Tally(object): ------ KeyError When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. ImportError When Pandas can not be found on the caller's system """ - # Ensure that StatePoint.read_results() was called first + # Ensure that the tally has data if self.mean is None or self.std_dev is None: - msg = 'The Tally ID="{0}" has no data to return. Call the ' \ - 'StatePoint.read_results() method before using ' \ - 'Tally.get_pandas_dataframe(...)'.format(self.id) + msg = 'The Tally ID="{0}" has no data to return'.format(self.id) raise KeyError(msg) # If using Summary, ensure StatePoint.link_with_summary(...) was called @@ -1200,12 +1293,12 @@ class Tally(object): # Tile the nuclide bins into a DataFrame column nuclides = np.repeat(nuclides, len(self.scores)) tile_factor = data_size / len(nuclides) - df['nuclide'] = np.tile(nuclides, tile_factor) + df['nuclide'] = np.tile(nuclides, int(tile_factor)) # Include column for scores if user requested it if scores: tile_factor = data_size / len(self.scores) - df['score'] = np.tile(self.scores, tile_factor) + df['score'] = np.tile(self.scores, int(tile_factor)) # Append columns with mean, std. dev. for each tally bin df['mean'] = self.mean.ravel() @@ -1275,7 +1368,7 @@ class Tally(object): for filter in self.filters: new_shape += (filter.num_bins, ) new_shape += (self.num_nuclides,) - new_shape += (self.num_score_bins,) + new_shape += (self.num_scores,) # Reshape the data with one dimension for each filter data = np.reshape(data, new_shape) @@ -1300,16 +1393,13 @@ class Tally(object): Raises ------ KeyError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ - # Ensure that StatePoint.read_results() was called first + # Ensure that the tally has data if self._sum is None or self._sum_sq is None and not self.derived: - msg = 'The Tally ID="{0}" has no data to export. Call the ' \ - 'StatePoint.read_results() method before using ' \ - 'Tally.export_results(...)'.format(self.id) + msg = 'The Tally ID="{0}" has no data to export'.format(self.id) raise KeyError(msg) if not isinstance(filename, basestring): @@ -1425,34 +1515,83 @@ class Tally(object): # Pickle the Tally results to a file pickle.dump(tally_results, open(filename, 'wb')) - def _outer_product(self, other, binary_op): + def hybrid_product(self, other, binary_op, filter_product=None, + nuclide_product=None, score_product=None): """Combines filters, scores and nuclides with another tally. - This is a helper method for the tally arithmetic methods. The filters, - scores and nuclides from both tallies are enumerated into all possible - combinations and expressed as CrossFilter, CrossScore and - CrossNuclide objects in the new derived tally. + This is a helper method for the tally arithmetic operator overloaded + methods. It is called a "hybrid product" because it performs a + combination of tensor (or Kronecker) and entrywise (or Hadamard) + products. The filters from both tallies are combined using an entrywise + (or Hadamard) product on matching filters. By default, if all nuclides + are identical in the two tallies, the entrywise product is performed + across nuclides; else the tensor product is performed. By default, if + all scores are identical in the two tallies, the entrywise product is + performed across scores; else the tensor product is performed. Users + can also call the method explicitly and specify the desired product. Parameters ---------- other : Tally - The tally on the right hand side of the outer product + The tally on the right hand side of the hybrid product binary_op : {'+', '-', '*', '/', '^'} - The binary operation in the outer product + The binary operation in the hybrid product + filter_product : {'tensor', 'entrywise' or None} + The type of product (tensor or entrywise) to be performed between + filter data. The default is the entrywise product. Currently only + the entrywise product is supported since a tally cannot contain + two of the same filter. + nuclide_product : {'tensor', 'entrywise' or None} + The type of product (tensor or entrywise) to be performed between + nuclide data. The default is the entrywise product if all nuclides + between the two tallies are the same; otherwise the default is + the tensor product. + score_product : {'tensor', 'entrywise' or None} + The type of product (tensor or entrywise) to be performed between + score data. The default is the entrywise product if all scores + between the two tallies are the same; otherwise the default is + the tensor product. Returns ------- Tally - A new Tally that is the outer product with this one. + A new Tally that is the hybrid product with this one. Raises ------ ValueError When this method is called before the other tally is populated - with data by the StatePoint.read_results() method. + with data. """ + # Set default value for filter product if it was not set + if filter_product is None: + filter_product = 'entrywise' + elif filter_product == 'tensor': + msg = 'Unable to perform Tally arithmetic with a tensor product' \ + 'for the filter data as this is not currently supported.' + raise ValueError(msg) + + # Set default value for nuclide product if it was not set + if nuclide_product is None: + if self.nuclides == other.nuclides: + nuclide_product = 'entrywise' + else: + nuclide_product = 'tensor' + + # Set default value for score product if it was not set + if score_product is None: + if self.scores == other.scores: + score_product = 'entrywise' + else: + score_product = 'tensor' + + # Check product types + cv.check_value('filter product', filter_product, _PRODUCT_TYPES) + cv.check_value('nuclide product', nuclide_product, _PRODUCT_TYPES) + cv.check_value('score product', score_product, _PRODUCT_TYPES) + # Check that results have been read if not other.derived and other.sum is None: msg = 'Unable to use tally arithmetic with Tally ID="{0}" ' \ @@ -1468,32 +1607,23 @@ class Tally(object): new_name = '({0} {1} {2})'.format(self.name, binary_op, other.name) new_tally.name = new_name + # Query the mean and std dev so the tally data is read in from file + # if it has not already been read in. + self.mean, self.std_dev, other.mean, other.std_dev + # Create copies of self and other tallies to rearrange for tally # arithmetic self_copy = copy.deepcopy(self) other_copy = copy.deepcopy(other) - # Find any shared filters between the two tallies - filter_intersect = [] - for filter in self_copy.filters: - if filter in other_copy.filters: - filter_intersect.append(filter) + self_copy.sparse = False + other_copy.sparse = False - # Align the shared filters in successive order - for i, filter in enumerate(filter_intersect): - self_index = self_copy.filters.index(filter) - other_index = other_copy.filters.index(filter) - - # If necessary, swap self filter - if self_index != i: - self_copy.swap_filters(filter, self_copy.filters[i], inplace=True) - - # If necessary, swap other filter - if other_index != i: - other_copy.swap_filters(filter, other_copy.filters[i], inplace=True) - - data = self_copy._align_tally_data(other_copy) + # Align the tally data based on desired hybrid product + data = self_copy._align_tally_data(other_copy, filter_product, + nuclide_product, score_product) + # Perform tally arithmetic operation if binary_op == '+': new_tally._mean = data['self']['mean'] + data['other']['mean'] new_tally._std_dev = np.sqrt(data['self']['std. dev.']**2 + @@ -1523,6 +1653,13 @@ class Tally(object): new_tally._std_dev = np.abs(new_tally.mean) * \ np.sqrt(first_term**2 + second_term**2) + # Convert any infs and nans to zero + new_tally._mean[np.isinf(new_tally._mean)] = 0 + new_tally._mean = np.nan_to_num(new_tally._mean) + new_tally._std_dev[np.isinf(new_tally._std_dev)] = 0 + new_tally._std_dev = np.nan_to_num(new_tally._std_dev) + + # Set tally attributes if self_copy.estimator == other_copy.estimator: new_tally.estimator = self_copy.estimator if self_copy.with_summary and other_copy.with_summary: @@ -1530,77 +1667,52 @@ class Tally(object): if self_copy.num_realizations == other_copy.num_realizations: new_tally.num_realizations = self_copy.num_realizations - # If filters are identical, simply reuse them in derived tally - if self_copy.filters == other_copy.filters: + # Add filters to the new tally + if filter_product == 'entrywise': for self_filter in self_copy.filters: new_tally.add_filter(self_filter) - - # Generate filter "outer products" for non-identical filters else: + all_filters = [self_copy.filters, other_copy.filters] + for self_filter, other_filter in itertools.product(*all_filters): + new_filter = CrossFilter(self_filter, other_filter, binary_op) + new_tally.add_filter(new_filter) - # Find the common longest sequence of shared filters - match = 0 - for self_filter, other_filter in zip(self_copy.filters, other_copy.filters): - if self_filter == other_filter: - match += 1 - else: - break + # Add nuclides to the new tally + if nuclide_product == 'entrywise': + for self_nuclide in self_copy.nuclides: + new_tally.add_nuclide(self_nuclide) + else: + all_nuclides = [self_copy.nuclides, other_copy.nuclides] + for self_nuclide, other_nuclide in itertools.product(*all_nuclides): + new_nuclide = \ + CrossNuclide(self_nuclide, other_nuclide, binary_op) + new_tally.add_nuclide(new_nuclide) - match_filters = self_copy.filters[:match] - cross_filters = [self_copy.filters[match:], other_copy.filters[match:]] - - # Simply reuse shared filters in derived tally - for filter in match_filters: - new_tally.add_filter(filter) - - # Use cross filters to combine non-shared filters in derived tally - if len(self_copy.filters) != match and len(other_copy.filters) == match: - for filter in cross_filters[0]: - new_tally.add_filter(filter) - elif len(self_copy.filters) == match and len(other_copy.filters) != match: - for filter in cross_filters[1]: - new_tally.add_filter(filter) - else: - for self_filter, other_filter in itertools.product(*cross_filters): - new_filter = CrossFilter(self_filter, other_filter, binary_op) - new_tally.add_filter(new_filter) - - # Generate score "outer products" - if self_copy.scores == other_copy.scores: - new_tally.num_score_bins = self_copy.num_score_bins + # Add scores to the new tally + if score_product == 'entrywise': for self_score in self_copy.scores: new_tally.add_score(self_score) else: - new_tally.num_score_bins = self_copy.num_score_bins * other_copy.num_score_bins all_scores = [self_copy.scores, other_copy.scores] for self_score, other_score in itertools.product(*all_scores): new_score = CrossScore(self_score, other_score, binary_op) new_tally.add_score(new_score) - # Generate nuclide "outer products" - if self_copy.nuclides == other_copy.nuclides: - for self_nuclide in self_copy.nuclides: - new_tally.nuclides.append(self_nuclide) - else: - all_nuclides = [self_copy.nuclides, other_copy.nuclides] - for self_nuclide, other_nuclide in itertools.product(*all_nuclides): - new_nuclide = CrossNuclide(self_nuclide, other_nuclide, binary_op) - new_tally.add_nuclide(new_nuclide) - # Correct each Filter's stride - stride = new_tally.num_nuclides * new_tally.num_score_bins + stride = new_tally.num_nuclides * new_tally.num_scores for filter in reversed(new_tally.filters): filter.stride = stride stride *= filter.num_bins return new_tally - def _align_tally_data(self, other): + def _align_tally_data(self, other, filter_product, nuclide_product, + score_product): """Aligns data from two tallies for tally arithmetic. This is a helper method to construct a dict of dicts of the "aligned" data arrays from each tally for tally arithmetic. The method analyzes - the filters, scores and nuclides in both tally's and determines how to + the filters, scores and nuclides in both tallies and determines how to appropriately align the data for vectorized arithmetic. For example, if the two tallies have different filters, this method will use NumPy 'tile' and 'repeat' operations to the new data arrays such that all @@ -1611,6 +1723,15 @@ class Tally(object): ---------- other : Tally The tally to outer product with this tally + filter_product : {'entrywise'} + The type of product to be performed between filter data. Currently, + only the entrywise product is supported for the filter product. + nuclide_product : {'tensor', 'entrywise'} + The type of product (tensor or entrywise) to be performed between + nuclide data. + score_product : {'tensor', 'entrywise'} + The type of product (tensor or entrywise) to be performed between + score data. Returns ------- @@ -1620,113 +1741,147 @@ class Tally(object): """ - self_mean = copy.deepcopy(self.mean) - self_std_dev = copy.deepcopy(self.std_dev) - other_mean = copy.deepcopy(other.mean) - other_std_dev = copy.deepcopy(other.std_dev) + # Get the set of filters that each tally is missing + other_missing_filters = \ + set(self.filters).difference(set(other.filters)) + self_missing_filters = \ + set(other.filters).difference(set(self.filters)) - if self.filters != other.filters: + # Add filters present in self but not in other to other + for filter in other_missing_filters: + filter = copy.deepcopy(filter) + other._mean = np.repeat(other.mean, filter.num_bins, axis=0) + other._std_dev = np.repeat(other.std_dev, filter.num_bins, axis=0) + other.add_filter(filter) - # Determine the number of paired combinations of filter bins - # between the two tallies and repeat arrays along filter axes - diff1 = list(set(self.filters).difference(set(other.filters))) - diff2 = list(set(other.filters).difference(set(self.filters))) + # Add filters present in other but not in self to self + for filter in self_missing_filters: + filter = copy.deepcopy(filter) + self._mean = np.repeat(self.mean, filter.num_bins, axis=0) + self._std_dev = np.repeat(self.std_dev, filter.num_bins, axis=0) + self.add_filter(filter) - # Determine the factors by which each tally operands' data arrays - # must be tiled or repeated for the tally outer product - other_tile_factor = 1 - self_repeat_factor = 1 - for filter in diff1: - other_tile_factor *= filter.num_bins - for filter in diff2: - self_repeat_factor *= filter.num_bins + # Align other filters with self filters + for i, filter in enumerate(self.filters): + other_index = other.filters.index(filter) - # Tile / repeat the tally data for the tally outer product - self_shape = list(self_mean.shape) - other_shape = list(other_mean.shape) - self_shape[0] *= self_repeat_factor - self_mean = np.repeat(self_mean, self_repeat_factor) - self_std_dev = np.repeat(self_std_dev, self_repeat_factor) + # If necessary, swap other filter + if other_index != i: + other._swap_filters(filter, other.filters[i]) - if self_repeat_factor == 1: - other_shape[0] *= other_tile_factor - other_mean = np.repeat(other_mean, other_tile_factor, axis=0) - other_std_dev = np.repeat(other_std_dev, other_tile_factor, - axis=0) - else: - other_mean = np.tile(other_mean, (other_tile_factor, 1, 1)) - other_std_dev = np.tile(other_std_dev, (other_tile_factor, 1, 1)) + # Repeat and tile the data by nuclide in preparation for performing + # the tensor product across nuclides. + if nuclide_product == 'tensor': + self._mean = \ + np.repeat(self.mean, other.num_nuclides, axis=1) + self._std_dev = \ + np.repeat(self.std_dev, other.num_nuclides, axis=1) + other._mean = \ + np.tile(other.mean, (1, self.num_nuclides, 1)) + other._std_dev = \ + np.tile(other.std_dev, (1, self.num_nuclides, 1)) - # NumPy repeat and tile routines return 1D flattened arrays - # Reshape arrays as 3D with filters, nuclides and scores axes - self_mean.shape = tuple(self_shape) - self_std_dev.shape = tuple(self_shape) - other_mean.shape = tuple(other_shape) - other_std_dev.shape = tuple(other_shape) + # Add nuclides to each tally such that each tally contains the complete + # set of nuclides necessary to perform an entrywise product. New + # nuclides added to a tally will have all their scores set to zero. + else: - if self.nuclides != other.nuclides: + # Get the set of nuclides that each tally is missing + other_missing_nuclides = \ + set(self.nuclides).difference(set(other.nuclides)) + self_missing_nuclides = \ + set(other.nuclides).difference(set(self.nuclides)) - # Determine the number of paired combinations of nuclides - # between the two tallies and repeat arrays along nuclide axes - self_repeat_factor = other.num_nuclides - other_tile_factor = self.num_nuclides + # Add nuclides present in self but not in other to other + for nuclide in other_missing_nuclides: + other._mean = \ + np.insert(other.mean, other.num_nuclides, 0, axis=1) + other._std_dev = \ + np.insert(other.std_dev, other.num_nuclides, 0, axis=1) + other.add_nuclide(nuclide) - # Tile / repeat the tally data for the tally outer product - self_shape = list(self_mean.shape) - other_shape = list(other_mean.shape) - self_shape[1] *= self_repeat_factor - other_shape[1] *= other_tile_factor - self_mean = np.repeat(self_mean, self_repeat_factor, axis=1) - other_mean = np.tile(other_mean, (1, other_tile_factor, 1)) - self_std_dev = np.repeat(self_std_dev, self_repeat_factor, axis=1) - other_std_dev = np.tile(other_std_dev, (1, other_tile_factor, 1)) + # Add nuclides present in other but not in self to self + for nuclide in self_missing_nuclides: + self._mean = \ + np.insert(self.mean, self.num_nuclides, 0, axis=1) + self._std_dev = \ + np.insert(self.std_dev, self.num_nuclides, 0, axis=1) + self.add_nuclide(nuclide) - # NumPy repeat and tile routines return 1D flattened arrays - # Reshape arrays as 3D with filters, nuclides and scores axes - self_mean.shape = tuple(self_shape) - self_std_dev.shape = tuple(self_shape) - other_mean.shape = tuple(other_shape) - other_std_dev.shape = tuple(other_shape) + # Align other nuclides with self nuclides + for i, nuclide in enumerate(self.nuclides): + other_index = other.get_nuclide_index(nuclide) - if self.scores != other.scores: + # If necessary, swap other nuclide + if other_index != i: + other._swap_nuclides(nuclide, other.nuclides[i]) - # Determine the number of paired combinations of score bins - # between the two tallies and repeat arrays along score axes - self_repeat_factor = other.num_score_bins - other_tile_factor = self.num_score_bins + # Repeat and tile the data by score in preparation for performing + # the tensor product across scores. + if score_product == 'tensor': + self._mean = np.repeat(self.mean, other.num_scores, axis=2) + self._std_dev = np.repeat(self.std_dev, other.num_scores, axis=2) + other._mean = np.tile(other.mean, (1, 1, self.num_scores)) + other._std_dev = np.tile(other.std_dev, (1, 1, self.num_scores)) - # Tile / repeat the tally data for the tally outer product - self_shape = list(self_mean.shape) - other_shape = list(other_mean.shape) - self_shape[2] *= self_repeat_factor - other_shape[2] *= other_tile_factor - self_mean = np.repeat(self_mean, self_repeat_factor, axis=2) - other_mean = np.tile(other_mean, (1, 1, other_tile_factor)) - self_std_dev = np.repeat(self_std_dev, self_repeat_factor, axis=2) - other_std_dev = np.tile(other_std_dev, (1, 1, other_tile_factor)) + # Add scores to each tally such that each tally contains the complete set + # of scores necessary to perform an entrywise product. New scores added + # to a tally will be set to zero. + else: - # NumPy repeat and tile routines return 1D flattened arrays - # Reshape arrays as 3D with filters, nuclides and scores axes - self_mean.shape = tuple(self_shape) - self_std_dev.shape = tuple(self_shape) - other_mean.shape = tuple(other_shape) - other_std_dev.shape = tuple(other_shape) + # Get the set of scores that each tally is missing + other_missing_scores = \ + set(self.scores).difference(set(other.scores)) + self_missing_scores = \ + set(other.scores).difference(set(self.scores)) + + # Add scores present in self but not in other to other + for score in other_missing_scores: + other._mean = np.insert(other.mean, other.num_scores, 0, axis=2) + other._std_dev = np.insert(other.std_dev, other.num_scores, 0, axis=2) + other.add_score(score) + + # Add scores present in other but not in self to self + for score in self_missing_scores: + self._mean = np.insert(self.mean, self.num_scores, 0, axis=2) + self._std_dev = np.insert(self.std_dev, self.num_scores, 0, axis=2) + self.add_score(score) + + # Align other scores with self scores + for i, score in enumerate(self.scores): + other_index = other.scores.index(score) + + # If necessary, swap other score + if other_index != i: + other._swap_scores(score, other.scores[i]) + + # Correct the stride for other filters + stride = other.num_nuclides * other.num_scores + for filter in reversed(other.filters): + filter.stride = stride + stride *= filter.num_bins + + # Correct the stride for self filters + stride = self.num_nuclides * self.num_scores + for filter in reversed(self.filters): + filter.stride = stride + stride *= filter.num_bins data = {} data['self'] = {} data['other'] = {} - data['self']['mean'] = self_mean - data['other']['mean'] = other_mean - data['self']['std. dev.'] = self_std_dev - data['other']['std. dev.'] = other_std_dev + data['self']['mean'] = self.mean + data['other']['mean'] = other.mean + data['self']['std. dev.'] = self.std_dev + data['other']['std. dev.'] = other.std_dev return data - def swap_filters(self, filter1, filter2, inplace=False): + def _swap_filters(self, filter1, filter2): """Reverse the ordering of two filters in this tally This is a helper method for tally arithmetic which helps align the data - in two tallies with shared filters. This method copies this tally and - reverses the order of the two filters. + in two tallies with shared filters. This method reverses the order of + the two filters in place. Parameters ---------- @@ -1736,21 +1891,11 @@ class Tally(object): filter2 : Filter The filter to swap with filter1 - inplace : bool, optional - Whether to perform operation inplace or return new tally with the - filters swapped. - - Returns - ------- - swap_tally - If inplace is false, a copy of this tally with the filters swapped. - Otherwise, nothing is returned. - Raises ------ ValueError If this is a derived tally or this method is called before the tally - is populated with data by the StatePoint.read_results() method. + is populated with data. """ @@ -1763,6 +1908,7 @@ class Tally(object): cv.check_type('filter1', filter1, Filter) cv.check_type('filter2', filter2, Filter) + # Check that the filters exist in the tally and are not the same if filter1 == filter2: msg = 'Unable to swap a filter with itself' raise ValueError(msg) @@ -1775,25 +1921,15 @@ class Tally(object): 'does not contain such a filter'.format(filter2.type, self.id) raise ValueError(msg) - # Create a copy of the tally that preserves the original data formatting - # throughout swapping process - tally_copy = copy.deepcopy(self) - - # Set the swap tally - if inplace: - swap_tally = self - else: - swap_tally = copy.deepcopy(self) - # Swap the filters in the copied version of this Tally - filter1_index = swap_tally.filters.index(filter1) - filter2_index = swap_tally.filters.index(filter2) - swap_tally.filters[filter1_index] = filter2 - swap_tally.filters[filter2_index] = filter1 + filter1_index = self.filters.index(filter1) + filter2_index = self.filters.index(filter2) + self.filters[filter1_index] = filter2 + self.filters[filter2_index] = filter1 # Update the strides for each of the filters - stride = swap_tally.num_nuclides * swap_tally.num_score_bins - for filter in reversed(swap_tally.filters): + stride = self.num_nuclides * self.num_scores + for filter in reversed(self.filters): filter.stride = stride stride *= filter.num_bins @@ -1804,49 +1940,167 @@ class Tally(object): else: filter1_bins = [(filter1.get_bin(i)) for i in range(filter1.num_bins)] - if filter1.type == 'distribcell': + if filter2.type == 'distribcell': filter2_bins = np.arange(filter2.num_bins) else: filter2_bins = [filter2.get_bin(i) for i in range(filter2.num_bins)] - # Adjust the sum data array to relect the new filter order - if swap_tally.sum is not None: - for bin1, bin2 in itertools.product(filter1_bins, filter2_bins): - filter_bins = [(bin1,), (bin2,)] - data = tally_copy.get_values( - filters=filters, filter_bins=filter_bins, value='sum') - indices = swap_tally.get_filter_indices(filters, filter_bins) - swap_tally.sum[indices, :, :] = data - - # Adjust the sum_sq data array to relect the new filter order - if swap_tally.sum_sq is not None: - for bin1, bin2 in itertools.product(filter1_bins, filter2_bins): - filter_bins = [(bin1,), (bin2,)] - data = tally_copy.get_values( - filters=filters, filter_bins=filter_bins, value='sum_sq') - indices = swap_tally.get_filter_indices(filters, filter_bins) - swap_tally.sum_sq[indices, :, :] = data - # Adjust the mean data array to relect the new filter order - if swap_tally.mean is not None: + if self.mean is not None: for bin1, bin2 in itertools.product(filter1_bins, filter2_bins): filter_bins = [(bin1,), (bin2,)] - data = tally_copy.get_values( + data = self.get_values( filters=filters, filter_bins=filter_bins, value='mean') - indices = swap_tally.get_filter_indices(filters, filter_bins) - swap_tally._mean[indices, :, :] = data + indices = self.get_filter_indices(filters, filter_bins) + self.mean[indices, :, :] = data # Adjust the std_dev data array to relect the new filter order - if swap_tally.std_dev is not None: + if self.std_dev is not None: for bin1, bin2 in itertools.product(filter1_bins, filter2_bins): filter_bins = [(bin1,), (bin2,)] - data = tally_copy.get_values( + data = self.get_values( filters=filters, filter_bins=filter_bins, value='std_dev') - indices = swap_tally.get_filter_indices(filters, filter_bins) - swap_tally._std_dev[indices, :, :] = data + indices = self.get_filter_indices(filters, filter_bins) + self.std_dev[indices, :, :] = data - if not inplace: - return swap_tally + def _swap_nuclides(self, nuclide1, nuclide2): + """Reverse the ordering of two nuclides in this tally + + This is a helper method for tally arithmetic which helps align the data + in two tallies with shared nuclides. This method reverses the order of + the two nuclides in place. + + Parameters + ---------- + nuclide1 : Nuclide + The nuclide to swap with nuclide2 + + nuclide2 : Nuclide + The nuclide to swap with nuclide1 + + Raises + ------ + ValueError + If this is a derived tally or this method is called before the tally + is populated with data. + + """ + + # Check that results have been read + if not self.derived and self.sum is None: + msg = 'Unable to use tally arithmetic with Tally ID="{0}" ' \ + 'since it does not contain any results.'.format(self.id) + raise ValueError(msg) + + cv.check_type('nuclide1', nuclide1, Nuclide) + cv.check_type('nuclide2', nuclide2, Nuclide) + + # Check that the nuclides exist in the tally and are not the same + if nuclide1 == nuclide2: + msg = 'Unable to swap a nuclide with itself' + raise ValueError(msg) + elif nuclide1 not in self.nuclides: + msg = 'Unable to swap nuclide1 "{0}" in Tally ID="{1}" since it ' \ + 'does not contain such a nuclide'\ + .format(nuclide1.name, self.id) + raise ValueError(msg) + elif nuclide2 not in self.nuclides: + msg = 'Unable to swap "{0}" nuclide2 in Tally ID="{1}" since it ' \ + 'does not contain such a nuclide'\ + .format(nuclide2.name, self.id) + raise ValueError(msg) + + # Swap the nuclides in the Tally + nuclide1_index = self.get_nuclide_index(nuclide1) + nuclide2_index = self.get_nuclide_index(nuclide2) + self.nuclides[nuclide1_index] = nuclide2 + self.nuclides[nuclide2_index] = nuclide1 + + # Adjust the mean data array to relect the new nuclide order + if self.mean is not None: + nuclide1_mean = self.mean[:, nuclide1_index, :].copy() + nuclide2_mean = self.mean[:, nuclide2_index, :].copy() + self.mean[:, nuclide2_index, :] = nuclide1_mean + self.mean[:, nuclide1_index, :] = nuclide2_mean + + # Adjust the std_dev data array to relect the new nuclide order + if self.std_dev is not None: + nuclide1_std_dev = self.std_dev[:, nuclide1_index, :].copy() + nuclide2_std_dev = self.std_dev[:, nuclide2_index, :].copy() + self.std_dev[:, nuclide2_index, :] = nuclide1_std_dev + self.std_dev[:, nuclide1_index, :] = nuclide2_std_dev + + def _swap_scores(self, score1, score2): + """Reverse the ordering of two scores in this tally + + This is a helper method for tally arithmetic which helps align the data + in two tallies with shared scores. This method reverses the order + of the two scores in place. + + Parameters + ---------- + score1 : str or CrossScore + The score to swap with score2 + + score2 : str or CrossScore + The score to swap with score1 + + Raises + ------ + ValueError + If this is a derived tally or this method is called before the tally + is populated with data. + + """ + + # Check that results have been read + if not self.derived and self.sum is None: + msg = 'Unable to use tally arithmetic with Tally ID="{0}" ' \ + 'since it does not contain any results.'.format(self.id) + raise ValueError(msg) + + # Check that the scores are valid + if not isinstance(score1, (basestring, CrossScore)): + msg = 'Unable to swap score1 "{0}" in Tally ID="{1}" since it is ' \ + 'not a string or CrossScore'.format(score1, self.id) + raise ValueError(msg) + elif not isinstance(score2, (basestring, CrossScore)): + msg = 'Unable to swap score2 "{0}" in Tally ID="{1}" since it is ' \ + 'not a string or CrossScore'.format(score2, self.id) + raise ValueError(msg) + + # Check that the scores exist in the tally and are not the same + if score1 == score2: + msg = 'Unable to swap a score with itself' + raise ValueError(msg) + elif score1 not in self.scores: + msg = 'Unable to swap score1 "{0}" in Tally ID="{1}" since it ' \ + 'does not contain such a score'.format(score1, self.id) + raise ValueError(msg) + elif score2 not in self.scores: + msg = 'Unable to swap score2 "{0}" in Tally ID="{1}" since it ' \ + 'does not contain such a score'.format(score2, self.id) + raise ValueError(msg) + + # Swap the scores in the Tally + score1_index = self.get_score_index(score1) + score2_index = self.get_score_index(score2) + self.scores[score1_index] = score2 + self.scores[score2_index] = score1 + + # Adjust the mean data array to relect the new nuclide order + if self.mean is not None: + score1_mean = self.mean[:, :, score1_index].copy() + score2_mean = self.mean[:, :, score2_index].copy() + self.mean[:, :, score2_index] = score1_mean + self.mean[:, :, score1_index] = score2_mean + + # Adjust the std_dev data array to relect the new nuclide order + if self.std_dev is not None: + score1_std_dev = self.std_dev[:, :, score1_index].copy() + score2_std_dev = self.std_dev[:, :, score2_index].copy() + self.std_dev[:, :, score2_index] = score1_std_dev + self.std_dev[:, :, score1_index] = score2_std_dev def __add__(self, other): """Adds this tally to another tally or scalar value. @@ -1879,8 +2133,7 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ @@ -1891,7 +2144,11 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._outer_product(other, binary_op='+') + new_tally = self.hybrid_product(other, binary_op='+') + + # If both tally operands were sparse, sparsify the new tally + if self.sparse and other.sparse: + new_tally.sparse = True elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -1903,7 +2160,6 @@ class Tally(object): new_tally.estimator = self.estimator new_tally.with_summary = self.with_summary new_tally.num_realization = self.num_realizations - new_tally.num_score_bins = self.num_score_bins for filter in self.filters: new_tally.add_filter(filter) @@ -1912,6 +2168,9 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) + # If this tally operand is sparse, sparsify the new tally + new_tally.sparse = self.sparse + else: msg = 'Unable to add "{0}" to Tally ID="{1}"'.format(other, self.id) raise ValueError(msg) @@ -1949,8 +2208,7 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ @@ -1961,7 +2219,11 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._outer_product(other, binary_op='-') + new_tally = self.hybrid_product(other, binary_op='-') + + # If both tally operands were sparse, sparsify the new tally + if self.sparse and other.sparse: + new_tally.sparse = True elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -1972,7 +2234,6 @@ class Tally(object): new_tally.estimator = self.estimator new_tally.with_summary = self.with_summary new_tally.num_realization = self.num_realizations - new_tally.num_score_bins = self.num_score_bins for filter in self.filters: new_tally.add_filter(filter) @@ -1981,6 +2242,9 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) + # If this tally operand is sparse, sparsify the new tally + new_tally.sparse = self.sparse + else: msg = 'Unable to subtract "{0}" from Tally ' \ 'ID="{1}"'.format(other, self.id) @@ -2019,8 +2283,7 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ @@ -2031,7 +2294,11 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._outer_product(other, binary_op='*') + new_tally = self.hybrid_product(other, binary_op='*') + + # If original tally operands were sparse, sparsify the new tally + if self.sparse and other.sparse: + new_tally.sparse = True elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -2042,7 +2309,6 @@ class Tally(object): new_tally.estimator = self.estimator new_tally.with_summary = self.with_summary new_tally.num_realization = self.num_realizations - new_tally.num_score_bins = self.num_score_bins for filter in self.filters: new_tally.add_filter(filter) @@ -2051,6 +2317,9 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) + # If this tally operand is sparse, sparsify the new tally + new_tally.sparse = self.sparse + else: msg = 'Unable to multiply Tally ID="{0}" ' \ 'by "{1}"'.format(self.id, other) @@ -2089,8 +2358,7 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ @@ -2101,7 +2369,11 @@ class Tally(object): raise ValueError(msg) if isinstance(other, Tally): - new_tally = self._outer_product(other, binary_op='/') + new_tally = self.hybrid_product(other, binary_op='/') + + # If original tally operands were sparse, sparsify the new tally + if self.sparse and other.sparse: + new_tally.sparse = True elif isinstance(other, Real): new_tally = Tally(name='derived') @@ -2112,7 +2384,6 @@ class Tally(object): new_tally.estimator = self.estimator new_tally.with_summary = self.with_summary new_tally.num_realization = self.num_realizations - new_tally.num_score_bins = self.num_score_bins for filter in self.filters: new_tally.add_filter(filter) @@ -2121,6 +2392,9 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) + # If this tally operand is sparse, sparsify the new tally + new_tally.sparse = self.sparse + else: msg = 'Unable to divide Tally ID="{0}" ' \ 'by "{1}"'.format(self.id, other) @@ -2162,8 +2436,7 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ @@ -2174,7 +2447,11 @@ class Tally(object): raise ValueError(msg) if isinstance(power, Tally): - new_tally = self._outer_product(power, binary_op='^') + new_tally = self.hybrid_product(power, binary_op='^') + + # If original tally operands were sparse, sparsify the new tally + if self.sparse and other.sparse: + new_tally.sparse = True elif isinstance(power, Real): new_tally = Tally(name='derived') @@ -2186,7 +2463,6 @@ class Tally(object): new_tally.estimator = self.estimator new_tally.with_summary = self.with_summary new_tally.num_realization = self.num_realizations - new_tally.num_score_bins = self.num_score_bins for filter in self.filters: new_tally.add_filter(filter) @@ -2195,6 +2471,9 @@ class Tally(object): for score in self.scores: new_tally.add_score(score) + # If original tally was sparse, sparsify the exponentiated tally + new_tally.sparse = self.sparse + else: msg = 'Unable to raise Tally ID="{0}" to ' \ 'power "{1}"'.format(self.id, power) @@ -2345,18 +2624,18 @@ class Tally(object): Raises ------ ValueError - When this method is called before the Tally is populated with data - by the StatePoint.read_results() method. + When this method is called before the Tally is populated with data. """ - # Ensure that StatePoint.read_results() was called first + # Ensure that the tally has data if not self.derived and self.sum is None: msg = 'Unable to use tally arithmetic with Tally ID="{0}" ' \ 'since it does not contain any results.'.format(self.id) raise ValueError(msg) new_tally = copy.deepcopy(self) + new_tally.sparse = False if self.sum is not None: new_sum = self.get_values(scores, filters, filter_bins, @@ -2372,7 +2651,7 @@ class Tally(object): new_tally._mean = new_mean if self.std_dev is not None: new_std_dev = self.get_values(scores, filters, filter_bins, - nuclides, 'std_dev') + nuclides, 'std_dev') new_tally._std_dev = new_std_dev # SCORES @@ -2388,7 +2667,6 @@ class Tally(object): # Loop over indices in reverse to remove excluded scores for score_index in reversed(score_indices): new_tally.remove_score(self.scores[score_index]) - new_tally.num_score_bins -= 1 # NUCLIDES if nuclides: @@ -2428,11 +2706,13 @@ class Tally(object): filter.num_bins = len(filter_bins[i]) # Correct each Filter's stride - stride = new_tally.num_nuclides * new_tally.num_score_bins + stride = new_tally.num_nuclides * new_tally.num_scores for filter in reversed(new_tally.filters): filter.stride = stride stride *= filter.num_bins + # If original tally was sparse, sparsify the sliced tally + new_tally.sparse = self.sparse return new_tally def summation(self, scores=[], filter_type=None, @@ -2535,6 +2815,8 @@ class Tally(object): filters[i] = CrossFilter(filters[i-1], filters[i], '+') tally_sum.add_filter(filters[-1]) + # If original tally was sparse, sparsify the tally summation + tally_sum.sparse = self.sparse return tally_sum def diagonalize_filter(self, new_filter): @@ -2571,12 +2853,6 @@ class Tally(object): new_tally = copy.deepcopy(self) new_tally.add_filter(new_filter) - # Determine the shape of data in the new diagonalized Tally - num_filter_bins = new_tally.num_filter_bins - num_nuclides = new_tally.num_nuclides - num_score_bins = new_tally.num_score_bins - new_shape = (num_filter_bins, num_nuclides, num_score_bins) - # Determine "base" indices along the new "diagonal", and the factor # by which the "base" indices should be repeated to account for all # other filter bins in the diagonalized tally @@ -2592,24 +2868,26 @@ class Tally(object): # Inject this Tally's data along the diagonal of the diagonalized Tally if self.sum is not None: - new_tally._sum = np.zeros(new_shape, dtype=np.float64) + new_tally._sum = np.zeros(new_tally.shape, dtype=np.float64) new_tally._sum[diag_indices, :, :] = self.sum if self.sum_sq is not None: - new_tally._sum_sq = np.zeros(new_shape, dtype=np.float64) + new_tally._sum_sq = np.zeros(new_tally.shape, dtype=np.float64) new_tally._sum_sq[diag_indices, :, :] = self.sum_sq if self.mean is not None: - new_tally._mean = np.zeros(new_shape, dtype=np.float64) + new_tally._mean = np.zeros(new_tally.shape, dtype=np.float64) new_tally._mean[diag_indices, :, :] = self.mean if self.std_dev is not None: - new_tally._std_dev = np.zeros(new_shape, dtype=np.float64) + new_tally._std_dev = np.zeros(new_tally.shape, dtype=np.float64) new_tally._std_dev[diag_indices, :, :] = self.std_dev # Correct each Filter's stride - stride = new_tally.num_nuclides * new_tally.num_score_bins + stride = new_tally.num_nuclides * new_tally.num_scores for filter in reversed(new_tally.filters): filter.stride = stride stride *= filter.num_bins + # If original tally was sparse, sparsify the diagonalized tally + new_tally.sparse = self.sparse return new_tally diff --git a/openmc/universe.py b/openmc/universe.py index 98367c381b..8417f27c71 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -63,6 +63,8 @@ class Cell(object): that is used to translate (shift) the universe. offsets : ndarray Array of offsets used for distributed cell searches + distribcell_index : int + Index of this cell in distribcell arrays """ @@ -76,6 +78,7 @@ class Cell(object): self._rotation = None self._translation = None self._offsets = None + self._distribcell_index = None def __eq__(self, other): if not isinstance(other, Cell): @@ -122,6 +125,8 @@ class Cell(object): string += '{0: <16}{1}{2}\n'.format('\tTranslation', '=\t', self._translation) string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self._offsets) + string += '{0: <16}{1}{2}\n'.format('\tDistribcell index', '=\t', + self._distribcell_index) return string @@ -164,6 +169,10 @@ class Cell(object): def offsets(self): return self._offsets + @property + def distribcell_index(self): + return self._distribcell_index + @id.setter def id(self, cell_id): if cell_id is None: @@ -231,6 +240,11 @@ class Cell(object): cv.check_type('cell region', region, Region) self._region = region + @distribcell_index.setter + def distribcell_index(self, ind): + cv.check_type('distribcell index', ind, Integral) + self._distribcell_index = ind + def add_surface(self, surface, halfspace): """Add a half-space to the list of half-spaces whose intersection defines the cell. @@ -271,7 +285,7 @@ class Cell(object): else: self.region = Intersection(self.region, region) - def get_offset(self, path, filter_offset): + def get_cell_instance(self, path, distribcell_index): # Get the current element and remove it from the list cell_id = path[0] path = path[1:] @@ -282,12 +296,12 @@ class Cell(object): # If the Cell is filled by a Universe elif self._type == 'fill': - offset = self._offsets[filter_offset-1] - offset += self._fill.get_offset(path, filter_offset) + offset = self.offsets[distribcell_index-1] + offset += self.fill.get_cell_instance(path, distribcell_index) # If the Cell is filled by a Lattice else: - offset = self._fill.get_offset(path, filter_offset) + offset = self.fill.get_cell_instance(path, distribcell_index) return offset @@ -591,7 +605,7 @@ class Universe(object): self._cells.clear() - def get_offset(self, path, filter_offset): + def get_cell_instance(self, path, distribcell_index): # Get the current element and remove it from the list path = path[1:] @@ -599,7 +613,7 @@ class Universe(object): cell_id = path[0] # Make a recursive call to the Cell within this Universe - offset = self._cells[cell_id].get_offset(path, filter_offset) + offset = self.cells[cell_id].get_cell_instance(path, distribcell_index) # Return the offset computed at all nested Universe levels return offset @@ -807,7 +821,7 @@ class Lattice(object): def universes(self, universes): cv.check_iterable_type('lattice universes', universes, Universe, min_depth=2, max_depth=3) - self._universes = universes + self._universes = np.asarray(universes) def get_unique_universes(self): """Determine all unique universes in the lattice @@ -1059,21 +1073,22 @@ class RectLattice(Lattice): cv.check_greater_than('lattice pitch', dim, 0.0) self._pitch = pitch - def get_offset(self, path, filter_offset): + def get_cell_instance(self, path, distribcell_index): # Get the current element and remove it from the list i = path[0] path = path[1:] # For 2D Lattices if len(self._dimension) == 2: - offset = self._offsets[i[1]-1, i[2]-1, 0, filter_offset-1] - offset += self._universes[i[1]][i[2]].get_offset(path, filter_offset) + offset = self._offsets[i[1]-1, i[2]-1, 0, distribcell_index-1] + offset += self._universes[i[1]-1][i[2]-1].get_cell_instance(path, + distribcell_index) # For 3D Lattices else: - offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, filter_offset-1] - offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_offset(path, - filter_offset) + offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, distribcell_index-1] + offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_cell_instance( + path, distribcell_index) return offset @@ -1119,7 +1134,7 @@ class RectLattice(Lattice): for z in range(self._dimension[2]): for y in range(self._dimension[1]): for x in range(self._dimension[0]): - universe = self._universes[x][y][z] + universe = self._universes[z][y][x] # Append Universe ID to the Lattice XML subelement universe_ids += '{0} '.format(universe._id) @@ -1137,7 +1152,7 @@ class RectLattice(Lattice): else: for y in range(self._dimension[1]): for x in range(self._dimension[0]): - universe = self._universes[x][y] + universe = self._universes[y][x] # Append Universe ID to Lattice XML subelement universe_ids += '{0} '.format(universe._id) diff --git a/setup.py b/setup.py index 907d80e31b..d401f5fbaf 100644 --- a/setup.py +++ b/setup.py @@ -37,6 +37,7 @@ if have_setuptools: # Optional dependencies 'extras_require': { 'pandas': ['pandas'], + 'sparse' : ['scipy'], 'vtk': ['vtk', 'silomesh'], 'validate': ['lxml'] }}) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index a58135e88a..2cc9078dc5 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -48,7 +48,7 @@ module ace_header integer :: MT ! ENDF MT value real(8) :: Q_value ! Reaction Q value integer :: multiplicity ! Number of secondary particles released - type(Tab1), allocatable :: multiplicity_E ! Energy-dependent neutron yield + type(Tab1), pointer :: multiplicity_E => null() ! Energy-dependent neutron yield integer :: threshold ! Energy grid index of threshold logical :: scatter_in_cm ! scattering system in center-of-mass? logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity @@ -104,6 +104,8 @@ module ace_header class(Reaction), intent(inout) :: this ! The Reaction object to clear + if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E) + if (associated(this % edist)) then call this % edist % clear() deallocate(this % edist) diff --git a/src/constants.F90 b/src/constants.F90 index e08fdf8b61..81418dcaf6 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -14,7 +14,7 @@ module constants integer, parameter :: REVISION_STATEPOINT = 14 integer, parameter :: REVISION_PARTICLE_RESTART = 1 integer, parameter :: REVISION_TRACK = 1 - integer, parameter :: REVISION_SUMMARY = 1 + integer, parameter :: REVISION_SUMMARY = 2 ! ============================================================================ ! ADJUSTABLE PARAMETERS @@ -45,10 +45,11 @@ module constants ! Maximum number of words in a single line, length of line, and length of ! single word - integer, parameter :: MAX_WORDS = 500 - integer, parameter :: MAX_LINE_LEN = 250 - integer, parameter :: MAX_WORD_LEN = 150 - integer, parameter :: MAX_FILE_LEN = 255 + integer, parameter :: MAX_WORDS = 500 + integer, parameter :: MAX_LINE_LEN = 250 + integer, parameter :: MAX_WORD_LEN = 150 + integer, parameter :: MAX_FILE_LEN = 255 + integer, parameter :: REGION_SPEC_LEN = 1000 ! Maximum number of external source spatial resamples to encounter before an ! error is thrown. diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 6ca13740ba..3dad3ba395 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -131,10 +131,13 @@ module geometry_header integer, allocatable :: offset (:) ! Distribcell offset for tally ! counter integer, allocatable :: region(:) ! Definition of spatial region as - ! Boolean expression of half-spaces + ! Boolean expression of half-spaces integer, allocatable :: rpn(:) ! Reverse Polish notation for region - ! expression - logical :: simple ! Is the region simple (intersections only) + ! expression + logical :: simple ! Is the region simple (intersections + ! only) + integer :: distribcell_index ! Index corresponding to this cell in + ! distribcell arrays ! Rotation matrix and translation vector real(8), allocatable :: translation(:) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index fc7a462e6a..01e50d983a 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -1898,7 +1898,7 @@ contains logical :: mpio integer :: hdf5_err - integer :: driver + integer(HID_T) :: driver integer(HID_T) :: file_id integer(HID_T) :: fapl_id diff --git a/src/initialize.F90 b/src/initialize.F90 index 7aa351b498..820367f821 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1073,12 +1073,12 @@ contains do i = 1, n_tallies t => tallies(i) - do j = 1, t%n_filters - filter => t%filters(j) + do j = 1, t % n_filters + filter => t % filters(j) - if (filter%type == FILTER_DISTRIBCELL) then - if (.not. cell_list%contains(filter%int_bins(1))) then - call cell_list%add(filter%int_bins(1)) + if (filter % type == FILTER_DISTRIBCELL) then + if (.not. cell_list % contains(filter % int_bins(1))) then + call cell_list % add(filter % int_bins(1)) end if end if @@ -1089,8 +1089,8 @@ contains ! to determine the number of offset tables to allocate do i = 1, n_universes univ => universes(i) - do j = 1, univ%n_cells - if (cell_list%contains(univ%cells(j))) then + do j = 1, univ % n_cells + if (cell_list % contains(univ % cells(j))) then n_maps = n_maps + 1 end if end do @@ -1109,32 +1109,14 @@ contains found(:,:) = .false. k = 1 + ! Search through universes for distributed cells and assign each one a + ! unique distribcell array index. do i = 1, n_universes univ => universes(i) - - do j = 1, univ%n_cells - - if (cell_list%contains(univ%cells(j))) then - - ! Loop over all tallies - do l = 1, n_tallies - t => tallies(l) - - do m = 1, t%n_filters - filter => t%filters(m) - - ! Loop over only distribcell filters - ! If filter points to cell we just found, set offset index - if (filter%type == FILTER_DISTRIBCELL) then - if (filter%int_bins(1) == univ%cells(j)) then - filter%offset = k - end if - end if - - end do - end do - - univ_list(k) = univ%id + do j = 1, univ % n_cells + if (cell_list % contains(univ % cells(j))) then + cells(univ % cells(j)) % distribcell_index = k + univ_list(k) = univ % id k = k + 1 end if end do @@ -1142,26 +1124,26 @@ contains ! Allocate the offset tables for lattices do i = 1, n_lattices - lat => lattices(i)%obj + lat => lattices(i) % obj select type(lat) type is (RectLattice) - allocate(lat%offset(n_maps, lat%n_cells(1), lat%n_cells(2), & - lat%n_cells(3))) + allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), & + lat % n_cells(3))) type is (HexLattice) - allocate(lat%offset(n_maps, 2 * lat%n_rings - 1, & - 2 * lat%n_rings - 1, lat%n_axial)) + allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, & + 2 * lat % n_rings - 1, lat % n_axial)) end select - lat%offset(:, :, :, :) = 0 + lat % offset(:, :, :, :) = 0 end do ! Allocate offset table for fill cells do i = 1, n_cells - if (cells(i)%material == NONE) then - allocate(cells(i)%offset(n_maps)) + if (cells(i) % material == NONE) then + allocate(cells(i) % offset(n_maps)) end if end do diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 866ed7b452..a6f836ac61 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1037,7 +1037,7 @@ contains logical :: boundary_exists character(MAX_LINE_LEN) :: filename character(MAX_WORD_LEN) :: word - character(1000) :: region_spec + character(REGION_SPEC_LEN) :: region_spec type(Cell), pointer :: c class(Surface), pointer :: s class(Lattice), pointer :: lat @@ -1092,8 +1092,9 @@ contains do i = 1, n_cells c => cells(i) - ! Initialize the number of cell instances - this is a base case for distribcells + ! Initialize distribcell instances and distribcell index c % instances = 0 + c % distribcell_index = NONE ! Get pointer to i-th cell node call get_list_item(node_cell_list, i, node_cell) diff --git a/src/output.F90 b/src/output.F90 index 3267dbd58c..5379555ff0 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1186,8 +1186,7 @@ contains label = '' univ => universes(BASE_UNIVERSE) offset = 0 - call find_offset(t % filters(i_filter) % offset, & - t % filters(i_filter) % int_bins(1), & + call find_offset(t % filters(i_filter) % int_bins(1), & univ, bin-1, offset, label) case (FILTER_SURFACE) i = t % filters(i_filter) % int_bins(bin) @@ -1221,15 +1220,15 @@ contains ! with the given offset !=============================================================================== - recursive subroutine find_offset(map, goal, univ, final, offset, path) + recursive subroutine find_offset(goal, univ, final, offset, path) - integer, intent(in) :: map ! Index in maps vector - integer, intent(in) :: goal ! The target cell ID + integer, intent(in) :: goal ! The target cell index type(Universe), intent(in) :: univ ! Universe to begin search integer, intent(in) :: final ! Target offset integer, intent(inout) :: offset ! Current offset character(*), intent(inout) :: path ! Path to offset + integer :: map ! Index in maps vector integer :: i, j ! Index over cells integer :: k, l, m ! Indices in lattice integer :: old_k, old_l, old_m ! Previous indices in lattice @@ -1244,6 +1243,9 @@ contains type(Universe), pointer :: next_univ ! Next universe to loop through class(Lattice), pointer :: lat ! Pointer to current lattice + ! Get the distribcell index for this cell + map = cells(goal) % distribcell_index + n = univ % n_cells ! Write to the geometry stack @@ -1255,17 +1257,13 @@ contains ! Look through all cells in this universe do i = 1, n - - cell_index = univ % cells(i) - c => cells(cell_index) - - ! If the cell ID matches the goal and the offset matches final, - ! write to the geometry stack - if (cell_dict % get_key(c % id) == goal .AND. offset == final) then - path = trim(path) // "->" // to_str(c%id) + ! If the cell matches the goal and the offset matches final, write to the + ! geometry stack + if (univ % cells(i) == goal .AND. offset == final) then + c => cells(univ % cells(i)) + path = trim(path) // "->" // to_str(c % id) return end if - end do ! Find the fill cell or lattice cell that we need to enter @@ -1345,7 +1343,7 @@ contains offset = c % offset(map) + offset next_univ => universes(c % fill) - call find_offset(map, goal, next_univ, final, offset, path) + call find_offset(goal, next_univ, final, offset, path) return ! ==================================================================== @@ -1385,7 +1383,7 @@ contains path = trim(path) // "(" // trim(to_str(k)) // & "," // trim(to_str(l)) // "," // & trim(to_str(m)) // ")" - call find_offset(map, goal, next_univ, final, offset, path) + call find_offset(goal, next_univ, final, offset, path) return else old_m = m @@ -1401,7 +1399,7 @@ contains path = trim(path) // "(" // trim(to_str(old_k)) // & "," // trim(to_str(old_l)) // "," // & trim(to_str(old_m)) // ")" - call find_offset(map, goal, next_univ, final, offset, path) + call find_offset(goal, next_univ, final, offset, path) return end if @@ -1446,8 +1444,7 @@ contains trim(to_str(k - lat % n_rings)) // "," // & trim(to_str(l - lat % n_rings)) // "," // & trim(to_str(m)) // ")" - call find_offset(map, goal, next_univ, final, offset, & - path) + call find_offset(goal, next_univ, final, offset, path) return else old_m = m @@ -1464,7 +1461,7 @@ contains trim(to_str(old_k - lat % n_rings)) // "," // & trim(to_str(old_l - lat % n_rings)) // "," // & trim(to_str(old_m)) // ")" - call find_offset(map, goal, next_univ, final, offset, path) + call find_offset(goal, next_univ, final, offset, path) return end if diff --git a/src/state_point.F90 b/src/state_point.F90 index 0a818ed297..4b58cb7704 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -274,7 +274,6 @@ contains call write_dataset(filter_group, "type", "delayedgroup") end select - call write_dataset(filter_group, "offset", tally%filters(j)%offset) call write_dataset(filter_group, "n_bins", tally%filters(j)%n_bins) if (tally % filters(j) % type == FILTER_ENERGYIN .or. & tally % filters(j) % type == FILTER_ENERGYOUT .or. & diff --git a/src/summary.F90 b/src/summary.F90 index cf6c39e8ff..feeb6e4903 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -120,7 +120,7 @@ contains integer(HID_T) :: universes_group, univ_group integer(HID_T) :: lattices_group, lattice_group real(8), allocatable :: coeffs(:) - character(MAX_LINE_LEN) :: region_spec + character(REGION_SPEC_LEN) :: region_spec type(Cell), pointer :: c class(Surface), pointer :: s type(Universe), pointer :: u @@ -203,6 +203,8 @@ contains end do call write_dataset(cell_group, "region", adjustl(region_spec)) + call write_dataset(cell_group, "distribcell_index", c % distribcell_index) + call close_group(cell_group) end do CELL_LOOP @@ -362,16 +364,22 @@ contains call write_dataset(lattice_group, "type", "rectangular") ! Write lattice dimensions, lower left corner, and pitch - call write_dataset(lattice_group, "dimension", lat%n_cells) - call write_dataset(lattice_group, "lower_left", lat%lower_left) + if (lat % is_3d) then + call write_dataset(lattice_group, "dimension", lat % n_cells) + call write_dataset(lattice_group, "lower_left", lat % lower_left) + else + call write_dataset(lattice_group, "dimension", lat % n_cells(1:2)) + call write_dataset(lattice_group, "lower_left", lat % lower_left) + end if ! Write lattice universes. allocate(lattice_universes(lat%n_cells(1), lat%n_cells(2), & &lat%n_cells(3))) do j = 1, lat%n_cells(1) - do k = 1, lat%n_cells(2) + do k = 0, lat%n_cells(2) - 1 do m = 1, lat%n_cells(3) - lattice_universes(j,k,m) = universes(lat%universes(j,k,m))%id + lattice_universes(j, k+1, m) = & + universes(lat%universes(j, lat%n_cells(2) - k, m))%id end do end do end do @@ -495,8 +503,10 @@ contains subroutine write_tallies(file_id) integer(HID_T), intent(in) :: file_id - integer :: i, j + integer :: i, j, k integer :: i_list, i_xs + integer :: n_order ! loop index for moment orders + integer :: nm_order ! loop index for Ynm moment orders integer(HID_T) :: tallies_group integer(HID_T) :: mesh_group integer(HID_T) :: tally_group @@ -551,7 +561,6 @@ contains filter_group = create_group(tally_group, "filter " // trim(to_str(j))) ! Write number of bins for this filter - call write_dataset(filter_group, "offset", t%filters(j)%offset) call write_dataset(filter_group, "n_bins", t%filters(j)%n_bins) ! Write filter bins @@ -675,6 +684,37 @@ contains deallocate(str_array) + ! Write explicit moment order strings for each score bin + k = 1 + allocate(str_array(t%n_score_bins)) + MOMENT_LOOP: do j = 1, t%n_user_score_bins + select case(t%score_bins(k)) + case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) + str_array(k) = 'P' // trim(to_str(t%moment_order(k))) + k = k + 1 + case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN) + do n_order = 0, t%moment_order(k) + str_array(k) = 'P' // trim(to_str(n_order)) + k = k + 1 + end do + case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, & + SCORE_TOTAL_YN) + do n_order = 0, t%moment_order(k) + do nm_order = -n_order, n_order + str_array(k) = 'Y' // trim(to_str(n_order)) // ',' // & + trim(to_str(nm_order)) + k = k + 1 + end do + end do + case default + str_array(k) = '' + k = k + 1 + end select + end do MOMENT_LOOP + + call write_dataset(tally_group, "moment_orders", str_array) + deallocate(str_array) + call close_group(tally_group) end do TALLY_METADATA diff --git a/src/tally.F90 b/src/tally.F90 index 85a603419b..f16a87000a 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -2276,6 +2276,7 @@ contains integer :: j integer :: n ! number of bins for single filter integer :: offset ! offset for distribcell + integer :: distribcell_index ! index in distribcell arrays real(8) :: E ! particle energy real(8) :: theta, phi ! Polar and Azimuthal Angles, respectively type(TallyObject), pointer :: t @@ -2320,12 +2321,14 @@ contains case (FILTER_DISTRIBCELL) ! determine next distribcell bin + distribcell_index = cells(t % filters(i) % int_bins(1)) & + % distribcell_index matching_bins(i) = NO_BIN_FOUND offset = 0 do j = 1, p % n_coord if (cells(p % coord(j) % cell) % type == CELL_FILL) then offset = offset + cells(p % coord(j) % cell) % & - offset(t % filters(i) % offset) + offset(distribcell_index) elseif(cells(p % coord(j) % cell) % type == CELL_LATTICE) then if (lattices(p % coord(j + 1) % lattice) % obj & % are_valid_indices([& @@ -2333,7 +2336,7 @@ contains p % coord(j + 1) % lattice_y, & p % coord(j + 1) % lattice_z])) then offset = offset + lattices(p % coord(j + 1) % lattice) % obj % & - offset(t % filters(i) % offset, & + offset(distribcell_index, & p % coord(j + 1) % lattice_x, & p % coord(j + 1) % lattice_y, & p % coord(j + 1) % lattice_z) diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 182c1e967d..8db9ab0e1b 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -55,7 +55,6 @@ module tally_header type TallyFilter integer :: type = NONE integer :: n_bins = 0 - integer :: offset = 0 ! Only used for distribcell filters integer, allocatable :: int_bins(:) real(8), allocatable :: real_bins(:) ! Only used for energy filters end type TallyFilter diff --git a/tests/test_many_scores/results_true.dat b/tests/test_many_scores/results_true.dat index 0d5dd6e32e..ae5430e30a 100644 --- a/tests/test_many_scores/results_true.dat +++ b/tests/test_many_scores/results_true.dat @@ -11,18 +11,6 @@ tally 1: 2.483728E+01 5.102293E-01 8.710841E-02 -8.628000E+00 -2.481430E+01 -9.329009E-01 -2.902534E-01 -5.102293E-01 -8.710841E-02 -5.087118E-01 -8.657086E-02 -8.632000E+00 -2.483728E+01 -9.328366E-01 -2.902108E-01 5.087118E-01 8.657086E-02 9.212024E+00 diff --git a/tests/test_many_scores/tallies.xml b/tests/test_many_scores/tallies.xml index 2df5597d04..b8a154f1ac 100644 --- a/tests/test_many_scores/tallies.xml +++ b/tests/test_many_scores/tallies.xml @@ -4,9 +4,9 @@ - flux total scatter nu-scatter scatter-2 scatter-p2 nu-scatter-2 - nu-scatter-p2 transport n1n absorption nu-fission kappa-fission - flux-y2 total-y2 scatter-y2 nu-scatter-y2 events delayed-nu-fission + flux total scatter nu-scatter scatter-2 nu-scatter-2 transport n1n + absorption nu-fission kappa-fission flux-y2 total-y2 scatter-y2 + nu-scatter-y2 events delayed-nu-fission diff --git a/tests/test_score_MT/inputs_true.dat b/tests/test_score_MT/inputs_true.dat index 3789a69cba..b361c28cf2 100644 --- a/tests/test_score_MT/inputs_true.dat +++ b/tests/test_score_MT/inputs_true.dat @@ -1 +1 @@ -63295b9d510370e65e63a3db627d47d286f5479e53c8eabeda9a5cb25ffe35becb636835aadad11691e34c23292fe11b5f688daee76d76ceac4b5dfd2f9ede4c \ No newline at end of file +7270299e4a4dde19d250825b0124fa36b60df847f046e058ccdfc74d4690beaaaaf387e050f596cb3445caf793d6a390ddb2d19650a3dccd4eb60056ae3b1477 \ No newline at end of file diff --git a/tests/test_score_MT/results_true.dat b/tests/test_score_MT/results_true.dat index 44656963f0..04b9c5ca50 100644 --- a/tests/test_score_MT/results_true.dat +++ b/tests/test_score_MT/results_true.dat @@ -7,10 +7,6 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -0.000000E+00 -0.000000E+00 -1.538090E-03 -1.381456E-06 1.538090E-03 1.381456E-06 3.974412E-01 @@ -19,16 +15,12 @@ tally 1: 3.796357E-01 5.252455E-06 2.290531E-11 -5.252455E-06 -2.290531E-11 3.359792E-02 2.267331E-04 2.459115E-02 1.233078E-04 0.000000E+00 0.000000E+00 -0.000000E+00 -0.000000E+00 7.004005E-05 1.748739E-09 8.104946E-02 @@ -42,18 +34,12 @@ tally 2: 0.000000E+00 0.000000E+00 0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 2.000000E-01 9.000000E-03 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -0.000000E+00 -0.000000E+00 2.000000E-02 2.000000E-04 0.000000E+00 @@ -64,8 +50,6 @@ tally 2: 0.000000E+00 0.000000E+00 0.000000E+00 -0.000000E+00 -0.000000E+00 tally 3: 0.000000E+00 0.000000E+00 @@ -73,10 +57,6 @@ tally 3: 0.000000E+00 0.000000E+00 0.000000E+00 -0.000000E+00 -0.000000E+00 -1.408027E-03 -1.849607E-06 1.408027E-03 1.849607E-06 3.946436E-01 @@ -85,16 +65,12 @@ tally 3: 3.382059E-01 2.179241E-05 4.749090E-10 -2.179241E-05 -4.749090E-10 3.146594E-02 2.159308E-04 4.278352E-02 4.094478E-04 0.000000E+00 0.000000E+00 -0.000000E+00 -0.000000E+00 1.736514E-04 1.245171E-08 7.989710E-02 diff --git a/tests/test_score_MT/test_score_MT.py b/tests/test_score_MT/test_score_MT.py index dae86d418c..fb8aa7cfb4 100644 --- a/tests/test_score_MT/test_score_MT.py +++ b/tests/test_score_MT/test_score_MT.py @@ -12,7 +12,6 @@ class ScoreMTTestHarness(PyAPITestHarness): filt = openmc.Filter(type='cell', bins=(10, 21, 22, 23)) tallies = [openmc.Tally(tally_id=i) for i in range(1, 4)] [t.add_filter(filt) for t in tallies] - [t.add_score('n2n') for t in tallies] [t.add_score('16') for t in tallies] [t.add_score('51') for t in tallies] [t.add_score('102') for t in tallies] diff --git a/tests/test_tally_arithmetic/inputs_true.dat b/tests/test_tally_arithmetic/inputs_true.dat new file mode 100644 index 0000000000..8e8838131a --- /dev/null +++ b/tests/test_tally_arithmetic/inputs_true.dat @@ -0,0 +1 @@ +df6318b76cd37a29ef9dd09da48a70c191ed07c1a2ceb75eb502ab35086c31250872406f22c2587edfcee16f67dd95c81db012ccd728710ed2509084438aea56 \ No newline at end of file diff --git a/tests/test_tally_arithmetic/results_true.dat b/tests/test_tally_arithmetic/results_true.dat new file mode 100644 index 0000000000..ded2efa665 --- /dev/null +++ b/tests/test_tally_arithmetic/results_true.dat @@ -0,0 +1,134 @@ +[[[ 8.90240785e-05 8.31464209e-11 4.41507090e-05 4.12357364e-11] + [ 7.29618709e-05 5.98879928e-05 3.61847984e-05 2.97009235e-05] + [ 4.69276830e-05 4.38293656e-11 2.35903380e-05 2.20328275e-11] + [ 3.84607356e-05 3.15690405e-05 1.93340411e-05 1.58696166e-05]] + + [[ 8.90240785e-05 8.31464209e-11 4.41507090e-05 4.12357364e-11] + [ 7.29618709e-05 5.98879928e-05 3.61847984e-05 2.97009235e-05] + [ 4.69276830e-05 4.38293656e-11 2.35903380e-05 2.20328275e-11] + [ 3.84607356e-05 3.15690405e-05 1.93340411e-05 1.58696166e-05]] + + [[ 8.90240785e-05 8.31464209e-11 4.41507090e-05 4.12357364e-11] + [ 7.29618709e-05 5.98879928e-05 3.61847984e-05 2.97009235e-05] + [ 4.69276830e-05 4.38293656e-11 2.35903380e-05 2.20328275e-11] + [ 3.84607356e-05 3.15690405e-05 1.93340411e-05 1.58696166e-05]] + + ..., + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]]][[[ 8.90240785e-05 8.31464209e-11 4.41507090e-05 4.12357364e-11] + [ 7.29618709e-05 5.98879928e-05 3.61847984e-05 2.97009235e-05] + [ 4.69276830e-05 4.38293656e-11 2.35903380e-05 2.20328275e-11] + [ 3.84607356e-05 3.15690405e-05 1.93340411e-05 1.58696166e-05]] + + [[ 8.90240785e-05 8.31464209e-11 4.41507090e-05 4.12357364e-11] + [ 7.29618709e-05 5.98879928e-05 3.61847984e-05 2.97009235e-05] + [ 4.69276830e-05 4.38293656e-11 2.35903380e-05 2.20328275e-11] + [ 3.84607356e-05 3.15690405e-05 1.93340411e-05 1.58696166e-05]] + + [[ 8.90240785e-05 8.31464209e-11 4.41507090e-05 4.12357364e-11] + [ 7.29618709e-05 5.98879928e-05 3.61847984e-05 2.97009235e-05] + [ 4.69276830e-05 4.38293656e-11 2.35903380e-05 2.20328275e-11] + [ 3.84607356e-05 3.15690405e-05 1.93340411e-05 1.58696166e-05]] + + ..., + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]]][[[ 7.29618709e-05 5.98879928e-05 3.61847984e-05 2.97009235e-05] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 7.29618709e-05 5.98879928e-05 3.61847984e-05 2.97009235e-05] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 7.29618709e-05 5.98879928e-05 3.61847984e-05 2.97009235e-05] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + ..., + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]]][[[ 0.00000000e+00 4.41507090e-05 0.00000000e+00] + [ 0.00000000e+00 3.61847984e-05 0.00000000e+00] + [ 0.00000000e+00 2.35903380e-05 0.00000000e+00] + [ 0.00000000e+00 1.93340411e-05 0.00000000e+00]] + + [[ 0.00000000e+00 4.41507090e-05 0.00000000e+00] + [ 0.00000000e+00 3.61847984e-05 0.00000000e+00] + [ 0.00000000e+00 2.35903380e-05 0.00000000e+00] + [ 0.00000000e+00 1.93340411e-05 0.00000000e+00]] + + [[ 0.00000000e+00 4.41507090e-05 0.00000000e+00] + [ 0.00000000e+00 3.61847984e-05 0.00000000e+00] + [ 0.00000000e+00 2.35903380e-05 0.00000000e+00] + [ 0.00000000e+00 1.93340411e-05 0.00000000e+00]] + + ..., + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]][[[ 0.00000000e+00 3.61847984e-05 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 3.61847984e-05 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 3.61847984e-05 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + ..., + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]] + + [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] + [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]] \ No newline at end of file diff --git a/tests/test_tally_arithmetic/test_tally_arithmetic.py b/tests/test_tally_arithmetic/test_tally_arithmetic.py new file mode 100644 index 0000000000..6643fd5a22 --- /dev/null +++ b/tests/test_tally_arithmetic/test_tally_arithmetic.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python + +import os +import sys +import glob +import hashlib +sys.path.insert(0, os.pardir) +from testing_harness import PyAPITestHarness +import openmc + + +class TallyArithmeticTestHarness(PyAPITestHarness): + def _build_inputs(self): + + # The summary.h5 file needs to be created to read in the tallies + self._input_set.settings.output = {'summary': True} + + # Initialize the tallies file + tallies_file = openmc.TalliesFile() + + # Initialize the nuclides + u235 = openmc.Nuclide('U-235') + u238 = openmc.Nuclide('U-238') + pu239 = openmc.Nuclide('Pu-239') + + # Initialize Mesh + mesh = openmc.Mesh(mesh_id=1) + mesh.type = 'regular' + mesh.dimension = [2, 2, 2] + mesh.lower_left = [-160.0, -160.0, -183.0] + mesh.upper_right = [160.0, 160.0, 183.0] + + # Initialize the filters + energy_filter = openmc.Filter(type='energy', bins=(0.0, 0.253e-6, + 1.0e-3, 1.0, 20.0)) + material_filter = openmc.Filter(type='material', bins=(1, 3)) + distrib_filter = openmc.Filter(type='distribcell', bins=(60)) + mesh_filter = openmc.Filter(type='mesh') + mesh_filter.mesh = mesh + + # Initialized the tallies + tally = openmc.Tally(name='tally 1') + tally.add_filter(material_filter) + tally.add_filter(energy_filter) + tally.add_filter(distrib_filter) + tally.add_score('nu-fission') + tally.add_score('total') + tally.add_nuclide(u235) + tally.add_nuclide(pu239) + tallies_file.add_tally(tally) + + tally = openmc.Tally(name='tally 2') + tally.add_filter(energy_filter) + tally.add_filter(mesh_filter) + tally.add_score('total') + tally.add_score('fission') + tally.add_nuclide(u238) + tally.add_nuclide(u235) + tallies_file.add_tally(tally) + tallies_file.add_mesh(mesh) + + # Export tallies to file + self._input_set.tallies = tallies_file + super(TallyArithmeticTestHarness, self)._build_inputs() + + def _get_results(self, hash_output=False): + """Digest info in the statepoint and return as a string.""" + + # Read the statepoint file. + statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0] + sp = openmc.StatePoint(statepoint) + + # Read the summary file. + summary = glob.glob(os.path.join(os.getcwd(), 'summary.h5'))[0] + su = openmc.Summary(summary) + sp.link_with_summary(su) + + # Load the tallies + tally_1 = sp.get_tally(name='tally 1') + tally_2 = sp.get_tally(name='tally 2') + + # Perform all the tally arithmetic operations and output results + outstr = '' + tally_3 = tally_1 * tally_2 + outstr += str(tally_3.mean) + + tally_3 = tally_1.hybrid_product(tally_2, '*', 'entrywise', 'tensor', + 'tensor') + outstr += str(tally_3.mean) + + tally_3 = tally_1.hybrid_product(tally_2, '*', 'entrywise', 'entrywise', + 'tensor') + outstr += str(tally_3.mean) + + tally_3 = tally_1.hybrid_product(tally_2, '*', 'entrywise', 'tensor', + 'entrywise') + outstr += str(tally_3.mean) + + tally_3 = tally_1.hybrid_product(tally_2, '*', 'entrywise', 'entrywise', + 'entrywise') + outstr += str(tally_3.mean) + + print(outstr) + + # Hash the results if necessary + if hash_output: + sha512 = hashlib.sha512() + sha512.update(outstr.encode('utf-8')) + outstr = sha512.hexdigest() + + return outstr + + def _cleanup(self): + super(TallyArithmeticTestHarness, self)._cleanup() + f = os.path.join(os.getcwd(), 'tallies.xml') + if os.path.exists(f): os.remove(f) + +if __name__ == '__main__': + harness = TallyArithmeticTestHarness('statepoint.10.h5', True) + harness.main()