From eca70ba3df84335c777c1f97a2f669fd6058ef04 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 3 Nov 2016 16:49:10 -0400 Subject: [PATCH 01/19] Fix PyAPI plot pixels bug; add aspect ratio --- openmc/universe.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index 545b05d5a1..fc71d9a2ec 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -153,7 +153,8 @@ class Universe(object): return [] def plot(self, center=(0., 0., 0.), width=(1., 1.), pixels=(200, 200), - basis='xy', color_by='cell', colors=None, filename=None, seed=None): + basis='xy', color_by='cell', colors=None, filename=None, seed=None, + aspect='auto'): """Display a slice plot of the universe. Parameters @@ -171,9 +172,9 @@ class Universe(object): colors : dict Assigns colors to specific materials or cells. Keys are instances of - :class:`Cell` or :class:`Material` and values are RGB 3-tuples or RGBA - 4-tuples. Red, green, blue, and alpha should all be floats in the - range [0.0, 1.0], for example: + :class:`Cell` or :class:`Material` and values are RGB 3-tuples or + RGBA 4-tuples. Red, green, blue, and alpha should all be floats in + the range [0.0, 1.0], for example: .. code-block:: python @@ -188,6 +189,9 @@ class Universe(object): Hashable object which is used to seed the random number generator used to select colors. If None, the generator is seeded from the current time. + aspect : 'auto' or Real or None + This argument is passed directly to matplotlib.pyplot.imshow. + 'auto' makes the image aspect ratio match the physical units. """ import matplotlib.pyplot as plt @@ -211,7 +215,8 @@ class Universe(object): y_min = center[1] - 0.5*width[1] y_max = center[1] + 0.5*width[1] elif basis == 'yz': - # The x-axis will correspond to physical y and the y-axis will correspond to physical z + # The x-axis will correspond to physical y and the y-axis will + # correspond to physical z x_min = center[1] - 0.5*width[0] x_max = center[1] + 0.5*width[0] y_min = center[2] - 0.5*width[1] @@ -229,8 +234,9 @@ class Universe(object): y_coords = np.linspace(y_max, y_min, pixels[1], endpoint=False) - \ 0.5*(y_max - y_min)/pixels[1] - # Search for locations and assign colors - img = np.zeros(pixels + (4,)) # Use RGBA form + # Initialize output image in RGBA format. Flip the pixels from + # traditional (x, y) to (y, x) used in graphics. + img = np.zeros((pixels[1], pixels[0], 4)) for i, x in enumerate(x_coords): for j, y in enumerate(y_coords): if basis == 'xy': @@ -257,7 +263,7 @@ class Universe(object): img[j, i, :] = colors[obj] # Display image - plt.imshow(img, extent=(x_min, x_max, y_min, y_max)) + plt.imshow(img, extent=(x_min, x_max, y_min, y_max), aspect=aspect) # Show or save the plot if filename is None: From 11040fb5f94cee76f551ef3681e0d2839c6d450a Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 7 Nov 2016 18:59:50 -0500 Subject: [PATCH 02/19] Added piecewise and combination functions --- openmc/data/function.py | 113 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 113 insertions(+) diff --git a/openmc/data/function.py b/openmc/data/function.py index e9ecf82c58..e87f0a9c42 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -397,6 +397,62 @@ class Polynomial(np.polynomial.Polynomial, Function1D): return cls(dataset.value) +class Combination(EqualityMixin): + """Combination of multiple functions with a user-defined operator + + This class allows you to create a callable object which represents the + combination of other callable objects by way of a series of user-defined + operators connecting each of the callable objects. + + Parameters + ---------- + functions : Iterable of Callable + Functions which are to be added together + operations : Iterable of numpy.ufunc + Operations to perform between functions; note that the standard order + of operations (i.e., PEMDAS) will not be followed. + + + Attributes + ---------- + functions : Iterable of Callable + Functions which are to be added together + operations : Iterable of numpy.ufunc + Operations to perform between functions + + """ + + def __init__(self, functions, operations): + self.functions = functions + self.operations = operations + + def __call__(self, x): + ans = np.zeros_like(x) + for i, operation in enumerate(self.operations): + ans = operation(ans, self.functions[i](x)) + return ans + + @property + def functions(self): + return self._functions + + @functions.setter + def functions(self, functions): + cv.check_type('functions', functions, Iterable, Callable) + self._functions = functions + + @property + def operations(self): + return self._operations + + @operations.setter + def operations(self, operations): + cv.check_type('operations', operations, Iterable, np.ufunc) + length = len(self.functions) + cv.check_length('operations', operations, length, length_max=length) + self._operations = operations + + class Sum(EqualityMixin): """Sum of multiple functions. @@ -432,6 +488,63 @@ class Sum(EqualityMixin): self._functions = functions +class Piecewise(EqualityMixin): + """Piecewise composition of multiple functions. + + This class allows you to create a callable object which is composed + of multiple other callable objects, each applying to a specific interval + + Parameters + ---------- + functions : Iterable of Callable + Functions which are to be combined in a piecewise fashion + breakpoints : Iterable of float + The breakpoints between each function. The functions + in the functions attribute must be provided in order of + increasing intervals. + + Attributes + ---------- + functions : Iterable of Callable + Functions which are to be combined in a piecewise fashion + breakpoints : Iterable of float + The breakpoints between each function + + """ + + def __init__(self, functions, breakpoints): + self.functions = functions + self.breakpoints = breakpoints + + def __call__(self, x): + i = np.searchsorted(self.breakpoints, x) + if isinstance(x, Iterable): + ans = np.empty_like(x) + for j in range(len(i)): + ans[j] = self.functions[i[j]](x[j]) + return ans + else: + return self.functions[i](x) + + @property + def functions(self): + return self._functions + + @property + def breakpoints(self): + return self._breakpoints + + @functions.setter + def functions(self, functions): + cv.check_type('functions', functions, Iterable, Callable) + self._functions = functions + + @breakpoints.setter + def breakpoints(self, breakpoints): + cv.check_iterable_type('breakpoints', breakpoints, Real) + self._breakpoints = breakpoints + + class ResonancesWithBackground(EqualityMixin): """Cross section in resolved resonance region. From fb512fb1ef1aee744b778c2247ba817f254d674e Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 7 Nov 2016 19:01:35 -0500 Subject: [PATCH 03/19] Added DataLibrary from_xml function --- openmc/data/library.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/openmc/data/library.py b/openmc/data/library.py index 0485af0641..3e01c24643 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -22,6 +22,12 @@ class DataLibrary(EqualityMixin): def __init__(self): self.libraries = [] + def __getitem__(self, material): + for library in self.libraries: + if material in library['materials']: + return library + return None + def register_file(self, filename): """Register a file with the data library. @@ -78,3 +84,33 @@ class DataLibrary(EqualityMixin): tree = ET.ElementTree(root) tree.write(path, xml_declaration=True, encoding='utf-8', method='xml') + + @classmethod + def from_xml(cls, path): + """Read cross section data library from an XML file. + + Parameters + ---------- + path : str + Path to XML file to read. + + """ + + data = cls() + + tree = ET.parse(path) + root = tree.getroot() + if root.find('directory') is not None: + directory = root.find('directory').text + else: + directory = os.path.dirname(path) + + for lib_element in root.findall('library'): + filename = os.path.join(directory, lib_element.attrib['path']) + filetype = lib_element.attrib['type'] + materials = lib_element.attrib['materials'].split() + library = {'path': filename, 'type': filetype, + 'materials': materials} + data.libraries.append(library) + + return data From df1d82708ed41134d57c3dc13828e76a5e45ea4b Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 8 Nov 2016 13:16:20 -0500 Subject: [PATCH 04/19] Make Universe.plot pass kwargs to matplotlib --- openmc/universe.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index fc71d9a2ec..260b9faf4c 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -154,7 +154,7 @@ class Universe(object): def plot(self, center=(0., 0., 0.), width=(1., 1.), pixels=(200, 200), basis='xy', color_by='cell', colors=None, filename=None, seed=None, - aspect='auto'): + **kwargs): """Display a slice plot of the universe. Parameters @@ -189,9 +189,11 @@ class Universe(object): Hashable object which is used to seed the random number generator used to select colors. If None, the generator is seeded from the current time. - aspect : 'auto' or Real or None - This argument is passed directly to matplotlib.pyplot.imshow. - 'auto' makes the image aspect ratio match the physical units. + + Keyword arguments + ----------------- + All keyword arguments are passed to matplotlib.pyplot.imshow. See the + matplotlib documentation for available kwargs. """ import matplotlib.pyplot as plt @@ -263,7 +265,7 @@ class Universe(object): img[j, i, :] = colors[obj] # Display image - plt.imshow(img, extent=(x_min, x_max, y_min, y_max), aspect=aspect) + plt.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs) # Show or save the plot if filename is None: From daa6b9a10b9bbb523d93995d3e0bb9fd28810589 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 8 Nov 2016 18:20:08 -0500 Subject: [PATCH 05/19] Link Universe.plot docstring to matplotlib docs --- docs/source/conf.py | 3 ++- openmc/universe.py | 8 +++----- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 75e621bc17..672217c84f 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -249,5 +249,6 @@ napoleon_use_ivar = True intersphinx_mapping = { 'python': ('https://docs.python.org/3', None), 'numpy': ('http://docs.scipy.org/doc/numpy/', None), - 'pandas': ('http://pandas.pydata.org/pandas-docs/stable/', None) + 'pandas': ('http://pandas.pydata.org/pandas-docs/stable/', None), + 'matplotlib': ('http://matplotlib.org/', None) } diff --git a/openmc/universe.py b/openmc/universe.py index 260b9faf4c..b2e7dbb920 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -189,11 +189,9 @@ class Universe(object): Hashable object which is used to seed the random number generator used to select colors. If None, the generator is seeded from the current time. - - Keyword arguments - ----------------- - All keyword arguments are passed to matplotlib.pyplot.imshow. See the - matplotlib documentation for available kwargs. + **kwargs + All keyword arguments are passed to + :func:`matplotlib.pyplot.imshow`. """ import matplotlib.pyplot as plt From cdae8f70f76f49c90792927170d11e2abd339010 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Tue, 8 Nov 2016 18:45:17 -0500 Subject: [PATCH 06/19] resolved @paulromano comments --- openmc/data/function.py | 12 ++++++------ openmc/data/library.py | 19 +++++++++++++++++-- 2 files changed, 23 insertions(+), 8 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index e87f0a9c42..4a1f9d7a38 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -428,8 +428,8 @@ class Combination(EqualityMixin): def __call__(self, x): ans = np.zeros_like(x) - for i, operation in enumerate(self.operations): - ans = operation(ans, self.functions[i](x)) + for op, func in zip(self.operations, self.functions): + ans = op(ans, func(x)) return ans @property @@ -488,7 +488,7 @@ class Sum(EqualityMixin): self._functions = functions -class Piecewise(EqualityMixin): +class Regions1D(EqualityMixin): """Piecewise composition of multiple functions. This class allows you to create a callable object which is composed @@ -499,9 +499,9 @@ class Piecewise(EqualityMixin): functions : Iterable of Callable Functions which are to be combined in a piecewise fashion breakpoints : Iterable of float - The breakpoints between each function. The functions - in the functions attribute must be provided in order of - increasing intervals. + The values of the dependent variable that define the domain of + each function. The *i*th and *(i+1)*th values are the limits of the + domain of the *i*th function. Values must be monotonically increasing. Attributes ---------- diff --git a/openmc/data/library.py b/openmc/data/library.py index 3e01c24643..0363de4fc5 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -22,9 +22,19 @@ class DataLibrary(EqualityMixin): def __init__(self): self.libraries = [] - def __getitem__(self, material): + def get_by_materials(self, value): + """Access a library entry by passing a value of the 'materials' + attribute of the library entry. + + Returns + ------- + library : dict + Dictionary summarizing cross section data from a single file; + the dictionary has keys 'path', 'type', and 'materials'. + + """ for library in self.libraries: - if material in library['materials']: + if value in library['materials']: return library return None @@ -94,6 +104,11 @@ class DataLibrary(EqualityMixin): path : str Path to XML file to read. + Returns + ------- + data : openmc.data.DataLibrary + openmc.data.DataLibrary object initialized from the provided XML + """ data = cls() From deb4de90a25928ff59985e12f5006dc172708dbf Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Tue, 8 Nov 2016 18:52:39 -0500 Subject: [PATCH 07/19] Whoops, missed one --- openmc/data/function.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index 4a1f9d7a38..b1cdb09de0 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -427,9 +427,9 @@ class Combination(EqualityMixin): self.operations = operations def __call__(self, x): - ans = np.zeros_like(x) - for op, func in zip(self.operations, self.functions): - ans = op(ans, func(x)) + ans = self.functions[0](x) + for i, operation in enumerate(self.operations): + ans = operation(ans, self.functions[i + 1](x)) return ans @property @@ -448,7 +448,7 @@ class Combination(EqualityMixin): @operations.setter def operations(self, operations): cv.check_type('operations', operations, Iterable, np.ufunc) - length = len(self.functions) + length = len(self.functions) - 1 cv.check_length('operations', operations, length, length_max=length) self._operations = operations From ac8909fa29ad8507ef99dcca38a22eea91422426 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Tue, 8 Nov 2016 20:15:19 -0500 Subject: [PATCH 08/19] minor comment change --- openmc/data/function.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index b1cdb09de0..c8fa15b768 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -410,7 +410,9 @@ class Combination(EqualityMixin): Functions which are to be added together operations : Iterable of numpy.ufunc Operations to perform between functions; note that the standard order - of operations (i.e., PEMDAS) will not be followed. + of operations will not be followed, but can be simulated by + combinations of Combination objects. The operations parameter must have + a length one less than the number of functions. Attributes From db1cf06a4d96e8bd13ae734c35103145002382aa Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 9 Nov 2016 12:56:03 -0600 Subject: [PATCH 09/19] Fix argument in h5sget_simple_extend_ndims_f --- src/hdf5_interface.F90 | 2 +- src/mgxs_header.F90 | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index c8892c53c4..2a0711954f 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -2589,7 +2589,7 @@ contains subroutine get_ndims(obj_id, ndims) integer(HID_T), intent(in) :: obj_id - integer(HID_T), intent(out) :: ndims + integer, intent(out) :: ndims integer :: hdf5_err integer :: type diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index fc5e8cb3e1..17104d5088 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -430,7 +430,8 @@ module mgxs_header ! in that conversion character(MAX_LINE_LEN) :: temp_str - integer(HID_T) :: xsdata, xsdata_grp, scatt_grp, ndims + integer(HID_T) :: xsdata, xsdata_grp, scatt_grp + integer :: ndims integer(HSIZE_T) :: dims(2) real(8), allocatable :: temp_arr(:), temp_2d(:, :) real(8), allocatable :: temp_beta(:, :) @@ -1113,7 +1114,8 @@ module mgxs_header ! in that conversion character(MAX_LINE_LEN) :: temp_str - integer(HID_T) :: xsdata, xsdata_grp, scatt_grp, ndims + integer(HID_T) :: xsdata, xsdata_grp, scatt_grp + integer :: ndims integer(HSIZE_T) :: dims(4) integer, allocatable :: int_arr(:) real(8), allocatable :: temp_1d(:), temp_3d(:, :, :) From 93597752f27d3b4e0255b911fe6e8423ceb54957 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Wed, 9 Nov 2016 18:19:13 -0500 Subject: [PATCH 10/19] Resolved @paulromano comments --- openmc/data/function.py | 9 ++++++--- openmc/data/library.py | 9 ++++----- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index c8fa15b768..0515a57c0c 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -407,7 +407,7 @@ class Combination(EqualityMixin): Parameters ---------- functions : Iterable of Callable - Functions which are to be added together + Functions to combine according to operations operations : Iterable of numpy.ufunc Operations to perform between functions; note that the standard order of operations will not be followed, but can be simulated by @@ -418,9 +418,12 @@ class Combination(EqualityMixin): Attributes ---------- functions : Iterable of Callable - Functions which are to be added together + Functions to combine according to operations operations : Iterable of numpy.ufunc - Operations to perform between functions + Operations to perform between functions; note that the standard order + of operations will not be followed, but can be simulated by + combinations of Combination objects. The operations parameter must have + a length one less than the number of functions. """ diff --git a/openmc/data/library.py b/openmc/data/library.py index 0363de4fc5..9f6d931467 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -22,13 +22,12 @@ class DataLibrary(EqualityMixin): def __init__(self): self.libraries = [] - def get_by_materials(self, value): - """Access a library entry by passing a value of the 'materials' - attribute of the library entry. + def get_by_material(self, value): + """Return the library dictionary containing a given material. Returns ------- - library : dict + library : dict or None Dictionary summarizing cross section data from a single file; the dictionary has keys 'path', 'type', and 'materials'. @@ -107,7 +106,7 @@ class DataLibrary(EqualityMixin): Returns ------- data : openmc.data.DataLibrary - openmc.data.DataLibrary object initialized from the provided XML + Data library object initialized from the provided XML """ From 80bf6c752839c083482d99975fb18cb2f5a2c0e9 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Wed, 9 Nov 2016 22:20:29 -0500 Subject: [PATCH 11/19] add get_molar_mass() method to material and fixed two other issues --- examples/python/pincell/build-xml.py | 2 +- openmc/data/data.py | 23 +++++++++++++++++++++++ openmc/element.py | 2 +- openmc/material.py | 28 ++++++++++++++++++++++++++++ 4 files changed, 53 insertions(+), 2 deletions(-) diff --git a/examples/python/pincell/build-xml.py b/examples/python/pincell/build-xml.py index 979614d635..2a550b9d73 100644 --- a/examples/python/pincell/build-xml.py +++ b/examples/python/pincell/build-xml.py @@ -51,7 +51,7 @@ o = openmc.Element('O') # Instantiate some Materials and register the appropriate Nuclides uo2 = openmc.Material(material_id=1, name='UO2 fuel at 2.4% wt enrichment') uo2.set_density('g/cm3', 10.29769) -uo2.add_element(u, 1., enrichment=0.024) +uo2.add_element(u, 1., enrichment=2.4) uo2.add_element(o, 2.) helium = openmc.Material(material_id=2, name='Helium for gap') diff --git a/openmc/data/data.py b/openmc/data/data.py index 1f759d3d4d..5dee5aca73 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -1,5 +1,6 @@ import itertools import os +import re # Isotopic abundances from M. Berglund and M. E. Wieser, "Isotopic compositions @@ -149,6 +150,28 @@ def atomic_mass(isotope): """ if not _ATOMIC_MASS: + + # If the isotope represents all natural isotopes of the element + # (e.g. C0), set atomic mass manually using the values from Atomic + # weights of the elements 2013 (IUPAC Technical Report) + # (doi:10.1515/pac-2015-0305). In cases where an atomic mass range is + # given (e.g. C), the average value is used. + if int(re.split(r"(\d+)", isotope)[1]) == 0: + if re.split(r"(\d+)", isotope)[0] == 'C': + _ATOMIC_MASS[isotope.lower()] = 12.0106 + elif re.split(r"(\d+)", isotope)[0] == 'Zn': + _ATOMIC_MASS[isotope.lower()] = 65.38 + elif re.split(r"(\d+)", isotope)[0] == 'Pt': + _ATOMIC_MASS[isotope.lower()] = 195.084 + elif re.split(r"(\d+)", isotope)[0] == 'Os': + _ATOMIC_MASS[isotope.lower()] = 190.23 + elif re.split(r"(\d+)", isotope)[0] == 'Tl': + _ATOMIC_MASS[isotope.lower()] = 204.3835 + else: + msg = 'Could not get the atomic mass for zero isotope'\ + ' {0}.'.format(isotope) + raise ValueError(msg) + # Load data from AME2012 file mass_file = os.path.join(os.path.dirname(__file__), 'mass.mas12') with open(mass_file, 'r') as ame: diff --git a/openmc/element.py b/openmc/element.py index d56c3b3507..ce0b5e4386 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -233,7 +233,7 @@ class Element(object): # Compute the ratio of the nuclide atomic masses to the element # atomic mass - if percent_type == 'wo': + if percent_type == 'wo' and len(abundances) > 1: # Compute the element atomic mass element_am = 0. diff --git a/openmc/material.py b/openmc/material.py index c265c90694..abd5467b81 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -4,6 +4,7 @@ from numbers import Real, Integral import warnings from xml.etree import ElementTree as ET import sys +import re from six import string_types @@ -581,6 +582,33 @@ class Material(object): return nuclides + def get_molar_mass(self): + """Returns the molar mass of the material + + Returns + ------- + molar_mass : float + The molar mass of the material + + """ + + # Get a list of all the nuclides, with elements expanded + nuclide_densities = self.get_nuclide_densities() + + # Using the sum of specified atomic or weight amounts as a basis, sum + # the mass and moles of the material + mass = 0. + moles = 0. + for nuc,vals in nuclide_densities.items(): + if vals[2] == 'ao': + mass += vals[1] * openmc.data.atomic_mass(nuc) + moles += vals[1] + else: + moles += vals[1] / openmc.data.atomic_mass(nuc) + mass += vals[1] + + return (mass / moles) + def _get_nuclide_xml(self, nuclide, distrib=False): xml_element = ET.Element("nuclide") xml_element.set("name", nuclide[0].name) From 9e9e401aaad77bad96630b5cfd8a06c4607cb332 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Wed, 9 Nov 2016 22:23:59 -0500 Subject: [PATCH 12/19] removed import of re in material.py --- openmc/material.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index abd5467b81..8003cc1929 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -4,7 +4,6 @@ from numbers import Real, Integral import warnings from xml.etree import ElementTree as ET import sys -import re from six import string_types @@ -607,7 +606,9 @@ class Material(object): moles += vals[1] / openmc.data.atomic_mass(nuc) mass += vals[1] - return (mass / moles) + # Compute and return the molar mass + molar_mass = mass / moles + return molar_mass def _get_nuclide_xml(self, nuclide, distrib=False): xml_element = ET.Element("nuclide") From e16b0388253b5507a93fe5675127357ce338672f Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Wed, 9 Nov 2016 22:29:22 -0500 Subject: [PATCH 13/19] removed modification to element.expand() --- openmc/element.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/element.py b/openmc/element.py index ce0b5e4386..d56c3b3507 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -233,7 +233,7 @@ class Element(object): # Compute the ratio of the nuclide atomic masses to the element # atomic mass - if percent_type == 'wo' and len(abundances) > 1: + if percent_type == 'wo': # Compute the element atomic mass element_am = 0. From 2208800981ec4a5e74026e400726ffaa9bfc7a33 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Wed, 9 Nov 2016 23:00:08 -0500 Subject: [PATCH 14/19] modified data.py to set the atomic masses for all zero isotopes --- openmc/data/data.py | 23 ++++++----------------- 1 file changed, 6 insertions(+), 17 deletions(-) diff --git a/openmc/data/data.py b/openmc/data/data.py index 5dee5aca73..2d1e471d26 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -1,6 +1,5 @@ import itertools import os -import re # Isotopic abundances from M. Berglund and M. E. Wieser, "Isotopic compositions @@ -151,26 +150,16 @@ def atomic_mass(isotope): """ if not _ATOMIC_MASS: - # If the isotope represents all natural isotopes of the element + # For the isotope represented by all natural isotopes of the element # (e.g. C0), set atomic mass manually using the values from Atomic # weights of the elements 2013 (IUPAC Technical Report) # (doi:10.1515/pac-2015-0305). In cases where an atomic mass range is # given (e.g. C), the average value is used. - if int(re.split(r"(\d+)", isotope)[1]) == 0: - if re.split(r"(\d+)", isotope)[0] == 'C': - _ATOMIC_MASS[isotope.lower()] = 12.0106 - elif re.split(r"(\d+)", isotope)[0] == 'Zn': - _ATOMIC_MASS[isotope.lower()] = 65.38 - elif re.split(r"(\d+)", isotope)[0] == 'Pt': - _ATOMIC_MASS[isotope.lower()] = 195.084 - elif re.split(r"(\d+)", isotope)[0] == 'Os': - _ATOMIC_MASS[isotope.lower()] = 190.23 - elif re.split(r"(\d+)", isotope)[0] == 'Tl': - _ATOMIC_MASS[isotope.lower()] = 204.3835 - else: - msg = 'Could not get the atomic mass for zero isotope'\ - ' {0}.'.format(isotope) - raise ValueError(msg) + _ATOMIC_MASS['c0'] = 12.0106 + _ATOMIC_MASS['zn0'] = 65.38 + _ATOMIC_MASS['pt0'] = 195.084 + _ATOMIC_MASS['os0'] = 190.23 + _ATOMIC_MASS['tl0'] = 204.3835 # Load data from AME2012 file mass_file = os.path.join(os.path.dirname(__file__), 'mass.mas12') From e375620e78627728b0d27094a83b3d9d7d1019b0 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Wed, 9 Nov 2016 23:02:22 -0500 Subject: [PATCH 15/19] fixed typo in data.py --- openmc/data/data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/data/data.py b/openmc/data/data.py index 2d1e471d26..e3a6d53601 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -150,7 +150,7 @@ def atomic_mass(isotope): """ if not _ATOMIC_MASS: - # For the isotope represented by all natural isotopes of the element + # For the isotopes representing all natural isotopes of their element # (e.g. C0), set atomic mass manually using the values from Atomic # weights of the elements 2013 (IUPAC Technical Report) # (doi:10.1515/pac-2015-0305). In cases where an atomic mass range is From 8394970acb28f2bbfd81795ea8cad5d989931680 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Thu, 10 Nov 2016 15:42:22 -0500 Subject: [PATCH 16/19] addressed PR comments --- openmc/data/data.py | 24 ++++++++++--------- openmc/material.py | 56 ++++++++++++++++++++++----------------------- 2 files changed, 40 insertions(+), 40 deletions(-) diff --git a/openmc/data/data.py b/openmc/data/data.py index e3a6d53601..a3e7158126 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -1,5 +1,6 @@ import itertools import os +import re # Isotopic abundances from M. Berglund and M. E. Wieser, "Isotopic compositions @@ -150,17 +151,6 @@ def atomic_mass(isotope): """ if not _ATOMIC_MASS: - # For the isotopes representing all natural isotopes of their element - # (e.g. C0), set atomic mass manually using the values from Atomic - # weights of the elements 2013 (IUPAC Technical Report) - # (doi:10.1515/pac-2015-0305). In cases where an atomic mass range is - # given (e.g. C), the average value is used. - _ATOMIC_MASS['c0'] = 12.0106 - _ATOMIC_MASS['zn0'] = 65.38 - _ATOMIC_MASS['pt0'] = 195.084 - _ATOMIC_MASS['os0'] = 190.23 - _ATOMIC_MASS['tl0'] = 204.3835 - # Load data from AME2012 file mass_file = os.path.join(os.path.dirname(__file__), 'mass.mas12') with open(mass_file, 'r') as ame: @@ -171,6 +161,18 @@ def atomic_mass(isotope): line[100:106] + '.' + line[107:112]) _ATOMIC_MASS[name.lower()] = mass + # For isotopes found in some libraries that represent all natural + # isotopes of their element (e.g. C0), calculate the atomic mass as + # the sum of the atomic mass times the natural abudance of the isotopes + # that make up the element. + for element in ['C', 'Zn', 'Pt', 'Os', 'Tl']: + isotope_zero = element.lower() + '0' + _ATOMIC_MASS[isotope_zero] = 0. + for iso, abundance in NATURAL_ABUNDANCE.items(): + if re.match(r'{}\d+'.format(element), iso): + _ATOMIC_MASS[isotope_zero] += abundance * \ + _ATOMIC_MASS[iso.lower()] + # Get rid of metastable information if '_' in isotope: isotope = isotope[:isotope.find('_')] diff --git a/openmc/material.py b/openmc/material.py index 8003cc1929..4b382c5e34 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -64,6 +64,12 @@ class Material(object): List in which each item is a 3-tuple consisting of an :class:`openmc.Nuclide` instance, the percent density, and the percent type ('ao' or 'wo'). + molar_mass : float + The molar mass of the material computed in units of grams per mole of + nuclides in a material. This entails that the molar mass does not depend + on the magnitude of the sum of atomic amounts of elements and nuclides + in the material. For instance, the molar mass of UO2 would be + ~90 g/mol. """ @@ -194,6 +200,27 @@ class Material(object): def distrib_otf_file(self): return self._distrib_otf_file + @property + def molar_mass(self): + + # Get a list of all the nuclides, with elements expanded + nuclide_densities = self.get_nuclide_densities() + + # Using the sum of specified atomic or weight amounts as a basis, sum + # the mass and moles of the material + mass = 0. + moles = 0. + for nuc, vals in nuclide_densities.items(): + if vals[2] == 'ao': + mass += vals[1] * openmc.data.atomic_mass(nuc) + moles += vals[1] + else: + moles += vals[1] / openmc.data.atomic_mass(nuc) + mass += vals[1] + + # Compute and return the molar mass + return mass / moles + @id.setter def id(self, material_id): @@ -581,35 +608,6 @@ class Material(object): return nuclides - def get_molar_mass(self): - """Returns the molar mass of the material - - Returns - ------- - molar_mass : float - The molar mass of the material - - """ - - # Get a list of all the nuclides, with elements expanded - nuclide_densities = self.get_nuclide_densities() - - # Using the sum of specified atomic or weight amounts as a basis, sum - # the mass and moles of the material - mass = 0. - moles = 0. - for nuc,vals in nuclide_densities.items(): - if vals[2] == 'ao': - mass += vals[1] * openmc.data.atomic_mass(nuc) - moles += vals[1] - else: - moles += vals[1] / openmc.data.atomic_mass(nuc) - mass += vals[1] - - # Compute and return the molar mass - molar_mass = mass / moles - return molar_mass - def _get_nuclide_xml(self, nuclide, distrib=False): xml_element = ET.Element("nuclide") xml_element.set("name", nuclide[0].name) From 4dd663dd781560802f8e4374a2ddf769a71894db Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Thu, 10 Nov 2016 17:19:39 -0500 Subject: [PATCH 17/19] changed molar_mass to average_molar_mass --- openmc/material.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 4b382c5e34..70a070ddd6 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -64,12 +64,10 @@ class Material(object): List in which each item is a 3-tuple consisting of an :class:`openmc.Nuclide` instance, the percent density, and the percent type ('ao' or 'wo'). - molar_mass : float - The molar mass of the material computed in units of grams per mole of - nuclides in a material. This entails that the molar mass does not depend - on the magnitude of the sum of atomic amounts of elements and nuclides - in the material. For instance, the molar mass of UO2 would be - ~90 g/mol. + average_molar_mass : float + The average molar mass of nuclides in the material in units of grams per + mol. For example, UO2 with 3 nuclides will have an average molar mass + of 270 / 3 = 90 g / mol. """ @@ -201,7 +199,7 @@ class Material(object): return self._distrib_otf_file @property - def molar_mass(self): + def average_molar_mass(self): # Get a list of all the nuclides, with elements expanded nuclide_densities = self.get_nuclide_densities() From fb009c6b5e1d5fe5a29eb9707839fef613ba3d3b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 14 Nov 2016 11:26:08 -0600 Subject: [PATCH 18/19] Use multidimensional array for tally results rather than derived type --- docs/source/io_formats/statepoint.rst | 10 +- openmc/particle_restart.py | 8 +- openmc/statepoint.py | 13 +- openmc/summary.py | 10 +- openmc/tallies.py | 13 +- openmc/volume.py | 3 +- src/cmfd_data.F90 | 39 +++-- src/cmfd_execute.F90 | 12 +- src/constants.F90 | 6 + src/eigenvalue.F90 | 14 +- src/finalize.F90 | 4 +- src/global.F90 | 16 +- src/hdf5_interface.F90 | 130 ----------------- src/initialize.F90 | 53 +------ src/input_xml.F90 | 1 - src/math.F90 | 3 +- src/output.F90 | 95 ++++++------ src/physics_mg.F90 | 8 +- src/simulation.F90 | 31 ++-- src/state_point.F90 | 50 +++---- src/tally.F90 | 201 ++++++++++++-------------- src/tally_header.F90 | 91 ++++++++++-- src/tally_initialize.F90 | 6 +- src/trigger.F90 | 16 +- 24 files changed, 331 insertions(+), 502 deletions(-) diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index e8b8f8b90b..7551f4293f 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -224,12 +224,12 @@ if run_mode == 'k-eigenvalue': Tallying moment orders for Legendre and spherical harmonic tally expansions (*e.g.*, 'P2', 'Y1,2', etc.). -**/tallies/tally /results** (Compound type) +**/tallies/tally /results** (*double[][][2]*) - Accumulated sum and sum-of-squares for each bin of the i-th tally. This is a - two-dimensional array, the first dimension of which represents combinations - of filter bins and the second dimensions of which represents scoring - bins. Each element of the array has fields 'sum' and 'sum_sq'. + Accumulated sum and sum-of-squares for each bin of the i-th tally. The first + dimension represents combinations of filter bins, the second dimensions + represents scoring bins, and the third dimension has two entries for the sum + and the sum-of-squares. **/source_present** (*int*) diff --git a/openmc/particle_restart.py b/openmc/particle_restart.py index a0410298b6..cb13265860 100644 --- a/openmc/particle_restart.py +++ b/openmc/particle_restart.py @@ -1,3 +1,5 @@ +import h5py + class Particle(object): """Information used to restart a specific particle that caused a simulation to fail. @@ -33,12 +35,6 @@ class Particle(object): """ def __init__(self, filename): - import h5py - if h5py.__version__ == '2.6.0': - raise ImportError("h5py 2.6.0 has a known bug which makes it " - "incompatible with OpenMC's HDF5 files. " - "Please switch to a different version.") - self._f = h5py.File(filename, 'r') # Ensure filetype and revision are correct diff --git a/openmc/statepoint.py b/openmc/statepoint.py index aaa5c0b840..e05b18ebed 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -5,6 +5,7 @@ import warnings import glob import numpy as np +import h5py import openmc import openmc.checkvalue as cv @@ -104,12 +105,6 @@ class StatePoint(object): """ def __init__(self, filename, autolink=True): - import h5py - if h5py.__version__ == '2.6.0': - raise ImportError("h5py 2.6.0 has a known bug which makes it " - "incompatible with OpenMC's HDF5 files. " - "Please switch to a different version.") - self._f = h5py.File(filename, 'r') # Ensure filetype and revision are correct @@ -209,13 +204,13 @@ class StatePoint(object): def global_tallies(self): if self._global_tallies is None: data = self._f['global_tallies'].value - gt = np.zeros_like(data, dtype=[ + gt = np.zeros(data.shape[0], dtype=[ ('name', 'a14'), ('sum', 'f8'), ('sum_sq', 'f8'), ('mean', 'f8'), ('std_dev', 'f8')]) gt['name'] = ['k-collision', 'k-absorption', 'k-tracklength', 'leakage'] - gt['sum'] = data['sum'] - gt['sum_sq'] = data['sum_sq'] + gt['sum'] = data[:,1] + gt['sum_sq'] = data[:,2] # Calculate mean and sample standard deviation of mean n = self.n_realizations diff --git a/openmc/summary.py b/openmc/summary.py index fd1689eb28..37509ec373 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -2,6 +2,7 @@ from collections import Iterable import re import numpy as np +import h5py import openmc from openmc.region import Region @@ -23,15 +24,6 @@ class Summary(object): """ def __init__(self, filename): - # A user may not have h5py, but they can still use the rest of the - # Python API so we'll only try to import h5py if the user actually inits - # a Summary object. - import h5py - if h5py.__version__ == '2.6.0': - raise ImportError("h5py 2.6.0 has a known bug which makes it " - "incompatible with OpenMC's HDF5 files. " - "Please switch to a different version.") - openmc.reset_auto_ids() if not filename.endswith(('.h5', '.hdf5')): diff --git a/openmc/tallies.py b/openmc/tallies.py index 6047131450..28da08badf 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -13,6 +13,7 @@ from xml.etree import ElementTree as ET from six import string_types import numpy as np +import h5py import openmc import openmc.checkvalue as cv @@ -269,20 +270,14 @@ class Tally(object): return None if not self._results_read: - import h5py - if h5py.__version__ == '2.6.0': - raise ImportError("h5py 2.6.0 has a known bug which makes it " - "incompatible with OpenMC's HDF5 files. " - "Please switch to a different version.") - # Open the HDF5 statepoint file f = h5py.File(self._sp_filename, 'r') # Extract Tally data from the file data = f['tallies/tally {0}/results'.format( self.id)].value - sum = data['sum'] - sum_sq = data['sum_sq'] + sum = data[:,:,0] + sum_sq = data[:,:,1] # Reshape the results arrays sum = np.reshape(sum, self.shape) @@ -1726,8 +1721,6 @@ class Tally(object): # HDF5 binary file if format == 'hdf5': - import h5py - filename = directory + '/' + filename + '.h5' if append: diff --git a/openmc/volume.py b/openmc/volume.py index 2a2a92f61e..b1391c687c 100644 --- a/openmc/volume.py +++ b/openmc/volume.py @@ -5,6 +5,7 @@ from warnings import warn import numpy as np import pandas as pd +import h5py import openmc import openmc.checkvalue as cv @@ -187,8 +188,6 @@ class VolumeCalculation(object): Results of the stochastic volume calculation """ - import h5py - with h5py.File(filename, 'r') as f: domain_type = f.attrs['domain_type'].decode() samples = f.attrs['samples'] diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index e1bf1f9c35..e5ed5df70e 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -164,7 +164,7 @@ contains * t%stride) + 1 ! Get flux - flux = t % results(1,score_index) % sum + flux = t % results(RESULT_SUM,1,score_index) cmfd % flux(h,i,j,k) = flux ! Detect zero flux, abort if located @@ -175,10 +175,10 @@ contains end if ! Get total rr and convert to total xs - cmfd % totalxs(h,i,j,k) = t % results(2,score_index) % sum / flux + cmfd % totalxs(h,i,j,k) = t % results(RESULT_SUM,2,score_index) / flux ! Get p1 scatter rr and convert to p1 scatter xs - cmfd % p1scattxs(h,i,j,k) = t % results(3,score_index) % sum / flux + cmfd % p1scattxs(h,i,j,k) = t % results(RESULT_SUM,3,score_index) / flux ! Calculate diffusion coefficient cmfd % diffcof(h,i,j,k) = ONE/(3.0_8*(cmfd % totalxs(h,i,j,k) - & @@ -211,19 +211,18 @@ contains * t%stride) + 1 ! Get scattering - cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /& + cmfd % scattxs(h,g,i,j,k) = t % results(RESULT_SUM,1,score_index) /& cmfd % flux(h,i,j,k) ! Get nu-fission - cmfd % nfissxs(h,g,i,j,k) = t % results(2,score_index) % sum /& + cmfd % nfissxs(h,g,i,j,k) = t % results(RESULT_SUM,2,score_index) /& cmfd % flux(h,i,j,k) ! Bank source cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + & - t % results(2,score_index) % sum + t % results(RESULT_SUM,2,score_index) cmfd % keff_bal = cmfd % keff_bal + & - t % results(2,score_index) % sum / & - dble(t % n_realizations) + t % results(RESULT_SUM,2,score_index) / t % n_realizations end do INGROUP @@ -243,67 +242,67 @@ contains matching_bins(i_filter_surf) = OUT_LEFT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(1,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = IN_LEFT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(2,h,i,j,k) = t % results(RESULT_SUM,1,score_index) ! Right surface matching_bins(i_filter_surf) = IN_RIGHT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(3,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = OUT_RIGHT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(4,h,i,j,k) = t % results(RESULT_SUM,1,score_index) ! Back surface matching_bins(i_filter_surf) = OUT_BACK score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(5,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = IN_BACK score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(6,h,i,j,k) = t % results(RESULT_SUM,1,score_index) ! Front surface matching_bins(i_filter_surf) = IN_FRONT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(7,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = OUT_FRONT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(8,h,i,j,k) = t % results(RESULT_SUM,1,score_index) ! Bottom surface matching_bins(i_filter_surf) = OUT_BOTTOM score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(9,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = IN_BOTTOM score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(10,h,i,j,k) = t % results(RESULT_SUM,1,score_index) ! Top surface matching_bins(i_filter_surf) = IN_TOP score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(11,h,i,j,k) = t % results(RESULT_SUM,1,score_index) matching_bins(i_filter_surf) = OUT_TOP score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 - cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current(12,h,i,j,k) = t % results(RESULT_SUM,1,score_index) end if TALLY diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index d2631254b1..a5d92002b8 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -365,22 +365,18 @@ contains subroutine cmfd_tally_reset() - use global, only: n_cmfd_tallies, cmfd_tallies + use global, only: cmfd_tallies use output, only: write_message - use tally, only: reset_result integer :: i ! loop counter ! Print message call write_message("CMFD tallies reset", 7) - ! Begin loop around CMFD tallies - do i = 1, n_cmfd_tallies - - ! Reset that tally + ! Reset CMFD tallies + do i = 1, size(cmfd_tallies) cmfd_tallies(i) % n_realizations = 0 - call reset_result(cmfd_tallies(i) % results) - + cmfd_tallies(i) % results(:,:,:) = ZERO end do end subroutine cmfd_tally_reset diff --git a/src/constants.F90 b/src/constants.F90 index ac157d7ad5..f2a9fb6f1a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -269,6 +269,12 @@ module constants ! ============================================================================ ! TALLY-RELATED CONSTANTS + ! Tally result entries + integer, parameter :: & + RESULT_VALUE = 1, & + RESULT_SUM = 2, & + RESULT_SUM_SQ = 3 + ! Tally type integer, parameter :: & TALLY_VOLUME = 1, & diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 713bbc351a..05c8190a65 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -374,7 +374,7 @@ contains subroutine calculate_generation_keff() ! Get keff for this generation by subtracting off the starting value - keff_generation = global_tallies(K_TRACKLENGTH) % value - keff_generation + keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation #ifdef MPI ! Combine values across all processors @@ -466,14 +466,14 @@ contains k_combined = ZERO ! Copy estimates of k-effective and its variance (not variance of the mean) - kv(1) = global_tallies(K_COLLISION) % sum / n - kv(2) = global_tallies(K_ABSORPTION) % sum / n - kv(3) = global_tallies(K_TRACKLENGTH) % sum / n - cov(1,1) = (global_tallies(K_COLLISION) % sum_sq - & + kv(1) = global_tallies(RESULT_SUM, K_COLLISION) / n + kv(2) = global_tallies(RESULT_SUM, K_ABSORPTION) / n + kv(3) = global_tallies(RESULT_SUM, K_TRACKLENGTH) / n + cov(1,1) = (global_tallies(RESULT_SUM_SQ, K_COLLISION) - & n * kv(1) * kv(1)) / (n - 1) - cov(2,2) = (global_tallies(K_ABSORPTION) % sum_sq - & + cov(2,2) = (global_tallies(RESULT_SUM_SQ, K_ABSORPTION) - & n * kv(2) * kv(2)) / (n - 1) - cov(3,3) = (global_tallies(K_TRACKLENGTH) % sum_sq - & + cov(3,3) = (global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH) - & n * kv(3) * kv(3)) / (n - 1) ! Calculate covariances based on sums with Bessel's correction diff --git a/src/finalize.F90 b/src/finalize.F90 index 86100d1950..83a5ac5fae 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -9,7 +9,7 @@ module finalize use message_passing #endif - use hdf5_interface, only: hdf5_bank_t, hdf5_tallyresult_t + use hdf5_interface, only: hdf5_bank_t use hdf5, only: h5tclose_f, h5close_f implicit none @@ -53,7 +53,6 @@ contains call free_memory() ! Release compound datatypes - call h5tclose_f(hdf5_tallyresult_t, hdf5_err) call h5tclose_f(hdf5_bank_t, hdf5_err) ! Close FORTRAN interface. @@ -62,7 +61,6 @@ contains #ifdef MPI ! Free all MPI types call MPI_TYPE_FREE(MPI_BANK, mpi_err) - call MPI_TYPE_FREE(MPI_TALLYRESULT, mpi_err) ! If MPI is in use and enabled, terminate it call MPI_FINALIZE(mpi_err) diff --git a/src/global.F90 b/src/global.F90 index 8111fd5144..eb4b33d7d6 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -1,5 +1,11 @@ module global + use, intrinsic :: ISO_C_BINDING + +#ifdef MPIF08 + use mpi_f08 +#endif + use bank_header, only: Bank use cmfd_header use constants @@ -14,15 +20,11 @@ module global use set_header, only: SetInt use surface_header, only: SurfaceContainer use source_header, only: SourceDistribution - use tally_header, only: TallyObject, TallyResult + use tally_header, only: TallyObject use trigger_header, only: KTrigger use timer_header, only: Timer use volume_header, only: VolumeCalculation -#ifdef MPIF08 - use mpi_f08 -#endif - implicit none ! ============================================================================ @@ -164,7 +166,7 @@ module global ! 3) track-length estimate of k-eff ! 4) leakage fraction - type(TallyResult), allocatable, target :: global_tallies(:) + real(C_DOUBLE), allocatable, target :: global_tallies(:,:) ! It is possible to protect accumulate operations on global tallies by using ! an atomic update. However, when multiple threads accumulate to the same @@ -272,10 +274,8 @@ module global integer :: mpi_err ! MPI error code #ifdef MPIF08 type(MPI_Datatype) :: MPI_BANK - type(MPI_Datatype) :: MPI_TALLYRESULT #else integer :: MPI_BANK ! MPI datatype for fission bank - integer :: MPI_TALLYRESULT ! MPI datatype for TallyResult #endif #ifdef _OPENMP diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 2a0711954f..5e14605ae4 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -16,7 +16,6 @@ module hdf5_interface use h5lt use error, only: fatal_error - use tally_header, only: TallyResult #ifdef PHDF5 use message_passing, only: MPI_COMM_WORLD, MPI_INFO_NULL #endif @@ -24,7 +23,6 @@ module hdf5_interface implicit none private - integer(HID_T), public :: hdf5_tallyresult_t ! Compound type for TallyResult integer(HID_T), public :: hdf5_bank_t ! Compound type for Bank integer(HID_T), public :: hdf5_integer8_t ! type for integer(8) @@ -42,8 +40,6 @@ module hdf5_interface module procedure write_long module procedure write_string module procedure write_string_1D - module procedure write_tally_result_1D - module procedure write_tally_result_2D end interface write_dataset interface read_dataset @@ -60,8 +56,6 @@ module hdf5_interface module procedure read_long module procedure read_string module procedure read_string_1D - module procedure read_tally_result_1D - module procedure read_tally_result_2D module procedure read_complex_2D end interface read_dataset @@ -2063,130 +2057,6 @@ contains call h5ltset_attribute_string_f(group_id, var, attr_type, attr_str, hdf5_err) end subroutine write_attribute_string -!=============================================================================== -! WRITE_TALLY_RESULT writes an OpenMC TallyResult type -!=============================================================================== - - subroutine write_tally_result_1D(group_id, name, buffer) - integer(HID_T), intent(in) :: group_id - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(in), target :: buffer(:) ! data to write - - integer(HSIZE_T) :: dims(1) - - dims(:) = shape(buffer) - call write_tally_result_1D_explicit(group_id, dims, name, buffer) - end subroutine write_tally_result_1D - - subroutine write_tally_result_1D_explicit(group_id, dims, name, buffer) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(1) - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(in), target :: buffer(dims(1)) - - integer :: hdf5_err - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr - - call h5screate_simple_f(1, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), hdf5_tallyresult_t, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - call h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_tally_result_1D_explicit - - subroutine write_tally_result_2D(group_id, name, buffer) - integer(HID_T), intent(in) :: group_id - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(in), target :: buffer(:,:) ! data to write - - integer(HSIZE_T) :: dims(2) - - dims(:) = shape(buffer) - call write_tally_result_2D_explicit(group_id, dims, name, buffer) - end subroutine write_tally_result_2D - - subroutine write_tally_result_2D_explicit(group_id, dims, name, buffer) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(2) - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(in), target :: buffer(dims(1),dims(2)) - - integer :: hdf5_err - integer(HID_T) :: dset ! data set handle - integer(HID_T) :: dspace ! data or file space handle - type(c_ptr) :: f_ptr - - call h5screate_simple_f(2, dims, dspace, hdf5_err) - call h5dcreate_f(group_id, trim(name), hdf5_tallyresult_t, & - dspace, dset, hdf5_err) - f_ptr = c_loc(buffer) - call h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end subroutine write_tally_result_2D_explicit - -!=============================================================================== -! READ_TALLY_RESULT reads OpenMC TallyResult data -!=============================================================================== - - subroutine read_tally_result_1D(group_id, name, buffer) - integer(HID_T), intent(in) :: group_id - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(inout), target :: buffer(:) ! read data here - - integer(HSIZE_T) :: dims(1) - - dims(:) = shape(buffer) - call read_tally_result_1D_explicit(group_id, dims, name, buffer) - end subroutine read_tally_result_1D - - subroutine read_tally_result_1D_explicit(group_id, dims, name, buffer) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(1) - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(inout), target :: buffer(dims(1)) - - integer :: hdf5_err - integer(HID_T) :: dset ! data set handle - type(c_ptr) :: f_ptr - - call h5dopen_f(group_id, trim(name), dset, hdf5_err) - f_ptr = c_loc(buffer) - call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - call h5dclose_f(dset, hdf5_err) - end subroutine read_tally_result_1D_explicit - - subroutine read_tally_result_2D(group_id, name, buffer) - integer(HID_T), intent(in) :: group_id - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(inout), target :: buffer(:,:) - - integer(HSIZE_T) :: dims(2) - - dims(:) = shape(buffer) - call read_tally_result_2D_explicit(group_id, dims, name, buffer) - end subroutine read_tally_result_2D - - subroutine read_tally_result_2D_explicit(group_id, dims, name, buffer) - integer(HID_T), intent(in) :: group_id - integer(HSIZE_T), intent(in) :: dims(2) - character(*), intent(in) :: name ! name of data - type(TallyResult), intent(inout), target :: buffer(dims(1),dims(2)) - - integer :: hdf5_err - integer(HID_T) :: dset ! data set handle - type(c_ptr) :: f_ptr - - call h5dopen_f(group_id, trim(name), dset, hdf5_err) - f_ptr = c_loc(buffer) - call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - call h5dclose_f(dset, hdf5_err) - end subroutine read_tally_result_2D_explicit - subroutine read_attribute_double(buffer, obj_id, name) real(8), intent(inout), target :: buffer integer(HID_T), intent(in) :: obj_id diff --git a/src/initialize.F90 b/src/initialize.F90 index 4fe6f8769c..d51f9ef8eb 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -12,7 +12,7 @@ module initialize &BASE_UNIVERSE use global use hdf5_interface, only: file_open, read_dataset, file_close, hdf5_bank_t,& - hdf5_tallyresult_t, hdf5_integer8_t + hdf5_integer8_t use input_xml, only: read_input_xml, cells_in_univ_dict, read_plots_xml use material_header, only: Material use mgxs_data, only: read_mgxs, create_macro_xs @@ -22,7 +22,7 @@ module initialize use state_point, only: load_state_point use string, only: to_str, starts_with, ends_with, str_to_int use summary, only: write_summary - use tally_header, only: TallyObject, TallyResult + use tally_header, only: TallyObject use tally_initialize,only: configure_tallies use tally_filter use tally, only: init_tally_routines @@ -169,21 +169,11 @@ contains integer :: bank_blocks(5) ! Count for each datatype #ifdef MPIF08 type(MPI_Datatype) :: bank_types(5) - type(MPI_Datatype) :: result_types(1) - type(MPI_Datatype) :: temp_type #else integer :: bank_types(5) ! Datatypes - integer :: result_types(1) ! Datatypes - integer :: temp_type ! temporary derived type #endif integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements - integer :: result_blocks(1) ! Count for each datatype - integer(MPI_ADDRESS_KIND) :: result_disp(1) ! Displacements - integer(MPI_ADDRESS_KIND) :: result_base_disp ! Base displacement - integer(MPI_ADDRESS_KIND) :: lower_bound ! Lower bound for TallyResult - integer(MPI_ADDRESS_KIND) :: extent ! Extent for TallyResult type(Bank) :: b - type(TallyResult) :: tr ! Indicate that MPI is turned on mpi_enabled = .true. @@ -221,36 +211,6 @@ contains call MPI_TYPE_CREATE_STRUCT(5, bank_blocks, bank_disp, & bank_types, MPI_BANK, mpi_err) call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) - - ! ========================================================================== - ! CREATE MPI_TALLYRESULT TYPE - - ! Determine displacements for MPI_BANK type - call MPI_GET_ADDRESS(tr%value, result_base_disp, mpi_err) - call MPI_GET_ADDRESS(tr%sum, result_disp(1), mpi_err) - - ! Adjust displacements - result_disp = result_disp - result_base_disp - - ! Define temporary type for TallyResult - result_blocks = (/ 2 /) - result_types = (/ MPI_REAL8 /) - call MPI_TYPE_CREATE_STRUCT(1, result_blocks, result_disp, result_types, & - temp_type, mpi_err) - - ! Adjust lower-bound and extent of type for tally score - lower_bound = 0 - extent = result_disp(1) + 16 - call MPI_TYPE_CREATE_RESIZED(temp_type, lower_bound, extent, & - MPI_TALLYRESULT, mpi_err) - - ! Commit derived type for tally scores - call MPI_TYPE_COMMIT(MPI_TALLYRESULT, mpi_err) - - ! Free temporary MPI type - call MPI_TYPE_FREE(temp_type, mpi_err) - - end subroutine initialize_mpi #endif !=============================================================================== @@ -259,7 +219,6 @@ contains subroutine hdf5_initialize() - type(TallyResult), target :: tmp(2) ! temporary TallyResult type(Bank), target :: tmpb(2) ! temporary Bank integer :: hdf5_err integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals @@ -268,14 +227,6 @@ contains ! Initialize FORTRAN interface. call h5open_f(hdf5_err) - ! Create the compound datatype for TallyResult - call h5tcreate_f(H5T_COMPOUND_F, h5offsetof(c_loc(tmp(1)), & - c_loc(tmp(2))), hdf5_tallyresult_t, hdf5_err) - call h5tinsert_f(hdf5_tallyresult_t, "sum", h5offsetof(c_loc(tmp(1)), & - c_loc(tmp(1)%sum)), H5T_NATIVE_DOUBLE, hdf5_err) - call h5tinsert_f(hdf5_tallyresult_t, "sum_sq", h5offsetof(c_loc(tmp(1)), & - c_loc(tmp(1)%sum_sq)), H5T_NATIVE_DOUBLE, hdf5_err) - ! Create compound type for xyz and uvw call h5tarray_create_f(H5T_NATIVE_DOUBLE, 1, dims, coordinates_t, hdf5_err) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index ff530c618d..53cb11f985 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2039,7 +2039,6 @@ contains type(Library), allocatable :: libraries(:) type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b) - character(MAX_LINE_LEN) :: temp_str real(8), allocatable :: material_temps(:) logical :: file_exists character(MAX_FILE_LEN) :: env_variable diff --git a/src/math.F90 b/src/math.F90 index 3ab8195aa0..11ace6976a 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -1,8 +1,9 @@ module math + use, intrinsic :: ISO_C_BINDING + use constants use random_lcg, only: prn - use ISO_C_BINDING implicit none diff --git a/src/output.F90 b/src/output.F90 index 8127725a68..6a8ba23b58 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -611,7 +611,7 @@ contains t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1) ! Adjust sum_sq - global_tallies(:) % sum_sq = t_value * global_tallies(:) % sum_sq + global_tallies(RESULT_SUM_SQ,:) = t_value * global_tallies(RESULT_SUM_SQ,:) ! Adjust combined estimator if (n_realizations > 3) then @@ -623,26 +623,26 @@ contains ! write global tallies if (n_realizations > 1) then if (run_mode == MODE_EIGENVALUE) then - write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) & - % sum, global_tallies(K_COLLISION) % sum_sq - write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) & - % sum, global_tallies(K_TRACKLENGTH) % sum_sq - write(ou,102) "k-effective (Absorption)", global_tallies(K_ABSORPTION) & - % sum, global_tallies(K_ABSORPTION) % sum_sq + write(ou,102) "k-effective (Collision)", global_tallies(RESULT_SUM, & + K_COLLISION), global_tallies(RESULT_SUM_SQ, K_COLLISION) + write(ou,102) "k-effective (Track-length)", global_tallies(RESULT_SUM, & + K_TRACKLENGTH), global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH) + write(ou,102) "k-effective (Absorption)", global_tallies(RESULT_SUM, & + K_ABSORPTION), global_tallies(RESULT_SUM_SQ, K_ABSORPTION) if (n_realizations > 3) write(ou,102) "Combined k-effective", k_combined end if - write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, & - global_tallies(LEAKAGE) % sum_sq + write(ou,102) "Leakage Fraction", global_tallies(RESULT_SUM, LEAKAGE), & + global_tallies(RESULT_SUM_SQ, LEAKAGE) else if (master) call warning("Could not compute uncertainties -- only one & &active batch simulated!") if (run_mode == MODE_EIGENVALUE) then - write(ou,103) "k-effective (Collision)", global_tallies(K_COLLISION) % sum - write(ou,103) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum - write(ou,103) "k-effective (Absorption)", global_tallies(K_ABSORPTION) % sum + write(ou,103) "k-effective (Collision)", global_tallies(RESULT_SUM, K_COLLISION) + write(ou,103) "k-effective (Track-length)", global_tallies(RESULT_SUM, K_TRACKLENGTH) + write(ou,103) "k-effective (Absorption)", global_tallies(RESULT_SUM, K_ABSORPTION) end if - write(ou,103) "Leakage Fraction", global_tallies(LEAKAGE) % sum + write(ou,103) "Leakage Fraction", global_tallies(RESULT_SUM, LEAKAGE) end if write(ou,*) @@ -765,7 +765,7 @@ contains end if ! Multiply uncertainty by t-value - t % results % sum_sq = t_value * t % results % sum_sq + t % results(RESULT_SUM_SQ,:,:) = t_value * t % results(RESULT_SUM_SQ,:,:) end if ! Write header block @@ -876,8 +876,8 @@ contains score_names(abs(t % score_bins(k))) write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,score_index,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index))) case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN) score_index = score_index - 1 do n_order = 0, t % moment_order(k) @@ -886,9 +886,8 @@ contains score_names(abs(t % score_bins(k))) write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) & - % sum_sq)) + to_str(t % results(RESULT_SUM,score_index,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index))) end do k = k + t % moment_order(k) case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, & @@ -902,9 +901,9 @@ contains // score_names(abs(t % score_bins(k))) write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index)& - % sum_sq)) + to_str(t % results(RESULT_SUM,score_index,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,score_index,& + filter_index))) end do end do k = k + (t % moment_order(k) + 1)**2 - 1 @@ -916,8 +915,8 @@ contains end if write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,score_index,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index))) end select end do indent = indent - 2 @@ -1020,16 +1019,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current on Left", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_LEFT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current on Left", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) ! Right Surface matching_bins(i_filter_surf) = OUT_RIGHT @@ -1037,16 +1036,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current on Right", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current on Right", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) if (n_dim >= 2) then @@ -1056,16 +1055,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current on Back", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_BACK filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current on Back", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) ! Front Surface matching_bins(i_filter_surf) = OUT_FRONT @@ -1073,16 +1072,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Net Current on Front", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Net Current on Front", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) end if if (n_dim == 3) then @@ -1092,16 +1091,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current on Bottom", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_BOTTOM filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current on Bottom", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) ! Top Surface matching_bins(i_filter_surf) = OUT_TOP @@ -1109,16 +1108,16 @@ contains * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current on Top", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) matching_bins(i_filter_surf) = IN_TOP filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current on Top", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) + to_str(t % results(RESULT_SUM,1,filter_index)), & + trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index))) end if end do end do diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 776c7905aa..082b337a48 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -119,16 +119,16 @@ contains ! Score implicit absorption estimate of keff !$omp atomic - global_tallies(K_ABSORPTION) % value = & - global_tallies(K_ABSORPTION) % value + p % absorb_wgt * & + global_tallies(RESULT_VALUE, K_ABSORPTION) = & + global_tallies(RESULT_VALUE, K_ABSORPTION) + p % absorb_wgt * & material_xs % nu_fission / material_xs % absorption else ! See if disappearance reaction happens if (material_xs % absorption > prn() * material_xs % total) then ! Score absorption estimate of keff !$omp atomic - global_tallies(K_ABSORPTION) % value = & - global_tallies(K_ABSORPTION) % value + p % wgt * & + global_tallies(RESULT_VALUE, K_ABSORPTION) = & + global_tallies(RESULT_VALUE, K_ABSORPTION) + p % wgt * & material_xs % nu_fission / material_xs % absorption p % alive = .false. diff --git a/src/simulation.F90 b/src/simulation.F90 index f59ce371cc..ef055a5493 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -20,8 +20,7 @@ module simulation use source, only: initialize_source, sample_external_source use state_point, only: write_state_point, write_source_point use string, only: to_str - use tally, only: synchronize_tallies, setup_active_usertallies, & - reset_result + use tally, only: synchronize_tallies, setup_active_usertallies use trigger, only: check_triggers use tracking, only: transport use volume_calc, only: run_volume_calculations @@ -220,7 +219,7 @@ contains if (ufs) call count_source_for_ufs() ! Store current value of tracklength k - keff_generation = global_tallies(K_TRACKLENGTH) % value + keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) end if end subroutine initialize_generation @@ -237,24 +236,24 @@ contains !$omp parallel !$omp critical if (run_mode == MODE_EIGENVALUE) then - global_tallies(K_COLLISION) % value = & - global_tallies(K_COLLISION) % value + global_tally_collision - global_tallies(K_ABSORPTION) % value = & - global_tallies(K_ABSORPTION) % value + global_tally_absorption - global_tallies(K_TRACKLENGTH) % value = & - global_tallies(K_TRACKLENGTH) % value + global_tally_tracklength + global_tallies(RESULT_VALUE, K_COLLISION) = & + global_tallies(RESULT_VALUE, K_COLLISION) + global_tally_collision + global_tallies(RESULT_VALUE, K_ABSORPTION) = & + global_tallies(RESULT_VALUE, K_ABSORPTION) + global_tally_absorption + global_tallies(RESULT_VALUE, K_TRACKLENGTH) = & + global_tallies(RESULT_VALUE, K_TRACKLENGTH) + global_tally_tracklength end if - global_tallies(LEAKAGE) % value = & - global_tallies(LEAKAGE) % value + global_tally_leakage + global_tallies(RESULT_VALUE, LEAKAGE) = & + global_tallies(RESULT_VALUE, LEAKAGE) + global_tally_leakage !$omp end critical ! reset private tallies if (run_mode == MODE_EIGENVALUE) then - global_tally_collision = 0 - global_tally_absorption = 0 - global_tally_tracklength = 0 + global_tally_collision = ZERO + global_tally_absorption = ZERO + global_tally_tracklength = ZERO end if - global_tally_leakage = 0 + global_tally_leakage = ZERO !$omp end parallel if (run_mode == MODE_EIGENVALUE) then @@ -302,7 +301,7 @@ contains ! Reset global tally results if (.not. active_batches) then - call reset_result(global_tallies) + global_tallies(:,:) = ZERO n_realizations = 0 end if diff --git a/src/state_point.F90 b/src/state_point.F90 index 6b7a467d4f..f9616f2285 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -353,7 +353,7 @@ contains ! Write sum and sum_sq for each bin tally_group = open_group(tallies_group, "tally " & // to_str(tally % id)) - call write_dataset(tally_group, "results", tally % results) + call tally % write_results_hdf5(tally_group) call close_group(tally_group) end do TALLY_RESULTS @@ -481,7 +481,7 @@ contains integer :: n_bins ! total number of bins integer(HID_T) :: tallies_group, tally_group real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results - real(8), target :: global_temp(2,N_GLOBAL_TALLIES) + real(8), target :: global_temp(3,N_GLOBAL_TALLIES) #ifdef MPI real(8) :: dummy ! temporary receive buffer for non-root reduces #endif @@ -489,7 +489,7 @@ contains type(ElemKeyValueII), pointer :: current type(ElemKeyValueII), pointer :: next type(TallyObject), pointer :: tally - type(TallyResult), allocatable :: tallyresult_temp(:,:) + type(TallyObject) :: dummy_tally ! ========================================================================== ! COLLECT AND WRITE GLOBAL TALLIES @@ -505,9 +505,8 @@ contains end if ! Copy global tallies into temporary array for reducing - n_bins = 2 * N_GLOBAL_TALLIES - global_temp(1,:) = global_tallies(:)%sum - global_temp(2,:) = global_tallies(:)%sum_sq + n_bins = 3 * N_GLOBAL_TALLIES + global_temp(:,:) = global_tallies(:,:) if (master) then ! The MPI_IN_PLACE specifier allows the master to copy values into a @@ -519,20 +518,11 @@ contains ! Transfer values to value on master if (current_batch == n_max_batches .or. satisfy_triggers) then - global_tallies(:)%sum = global_temp(1,:) - global_tallies(:)%sum_sq = global_temp(2,:) + global_tallies(:,:) = global_temp(:,:) end if - ! Put reduced value in temporary tally result - allocate(tallyresult_temp(N_GLOBAL_TALLIES, 1)) - tallyresult_temp(:,1)%sum = global_temp(1,:) - tallyresult_temp(:,1)%sum_sq = global_temp(2,:) - ! Write out global tallies sum and sum_sq - call write_dataset(file_id, "global_tallies", tallyresult_temp) - - ! Deallocate temporary tally result - deallocate(tallyresult_temp) + call write_dataset(file_id, "global_tallies", global_temp) else ! Receive buffer not significant at other processors #ifdef MPI @@ -568,15 +558,15 @@ contains tally => tallies(i) ! Determine size of tally results array - m = size(tally%results, 1) - n = size(tally%results, 2) + m = size(tally%results, 2) + n = size(tally%results, 3) n_bins = m*n*2 ! Allocate array for storing sums and sums of squares, but ! contiguously in memory for each allocate(tally_temp(2,m,n)) - tally_temp(1,:,:) = tally%results(:,:)%sum - tally_temp(2,:,:) = tally%results(:,:)%sum_sq + tally_temp(1,:,:) = tally%results(RESULT_SUM,:,:) + tally_temp(2,:,:) = tally%results(RESULT_SUM_SQ,:,:) if (master) then tally_group = open_group(tallies_group, "tally " // & @@ -592,20 +582,20 @@ contains ! At the end of the simulation, store the results back in the ! regular TallyResults array if (current_batch == n_max_batches .or. satisfy_triggers) then - tally%results(:,:)%sum = tally_temp(1,:,:) - tally%results(:,:)%sum_sq = tally_temp(2,:,:) + tally%results(RESULT_SUM,:,:) = tally_temp(1,:,:) + tally%results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:) end if ! Put in temporary tally result - allocate(tallyresult_temp(m,n)) - tallyresult_temp(:,:)%sum = tally_temp(1,:,:) - tallyresult_temp(:,:)%sum_sq = tally_temp(2,:,:) + allocate(dummy_tally % results(3,m,n)) + dummy_tally % results(RESULT_SUM,:,:) = tally_temp(1,:,:) + dummy_tally % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:) ! Write reduced tally results to file - call write_dataset(tally_group, "results", tally%results) + call dummy_tally % write_results_hdf5(tally_group) ! Deallocate temporary tally result - deallocate(tallyresult_temp) + deallocate(dummy_tally % results) else ! Receive buffer not significant at other processors #ifdef MPI @@ -771,7 +761,7 @@ contains call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.) ! Read global tally data - call read_dataset(file_id, "global_tallies", global_tallies) + call read_dataset(global_tallies, file_id, "global_tallies") ! Check if tally results are present tallies_group = open_group(file_id, "tallies") @@ -787,7 +777,7 @@ contains ! Read sum, sum_sq, and N for each bin tally_group = open_group(tallies_group, "tally " // & trim(to_str(tally % id))) - call read_dataset(tally_group, "results", tally % results) + call tally % read_results_hdf5(tally_group) call read_dataset(tally % n_realizations, tally_group, & "n_realizations") call close_group(tally_group) diff --git a/src/tally.F90 b/src/tally.F90 index 69e4e8ecdc..d4968d9d68 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1,5 +1,7 @@ module tally + use, intrinsic :: ISO_C_BINDING + #ifdef MPI use message_passing #endif @@ -18,7 +20,6 @@ module tally use output, only: header use particle_header, only: LocalCoord, Particle use string, only: to_str - use tally_header, only: TallyResult use tally_filter implicit none @@ -1965,8 +1966,8 @@ contains score = score * calc_pn(t % moment_order(i), p % mu) endif !$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score + t % results(RESULT_VALUE, score_index, filter_index) = & + t % results(RESULT_VALUE, score_index, filter_index) + score case(SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN) @@ -1982,10 +1983,9 @@ contains ! multiply score by the angular flux moments and store !$omp critical (score_general_scatt_yn) - t % results(score_index: score_index + num_nm - 1, filter_index) & - % value = t & - % results(score_index: score_index + num_nm - 1, filter_index)& - % value & + t % results(RESULT_VALUE, score_index: score_index + num_nm - 1, & + filter_index) = t % results(RESULT_VALUE, & + score_index: score_index + num_nm - 1, filter_index) & + score * calc_pn(n, p % mu) * calc_rn(n, p % last_uvw) !$omp end critical (score_general_scatt_yn) end do @@ -2011,10 +2011,9 @@ contains ! multiply score by the angular flux moments and store !$omp critical (score_general_flux_tot_yn) - t % results(score_index: score_index + num_nm - 1, filter_index) & - % value = t & - % results(score_index: score_index + num_nm - 1, filter_index)& - % value & + t % results(RESULT_VALUE, score_index: score_index + num_nm - 1, & + filter_index) = t % results(RESULT_VALUE, & + score_index: score_index + num_nm - 1, filter_index) & + score * calc_rn(n, uvw) !$omp end critical (score_general_flux_tot_yn) end do @@ -2031,8 +2030,8 @@ contains ! get the score and tally it !$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value & + t % results(RESULT_VALUE, score_index, filter_index) = & + t % results(RESULT_VALUE, score_index, filter_index) & + score * calc_pn(n, p % mu) end do i = i + t % moment_order(i) @@ -2040,8 +2039,8 @@ contains case default !$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score + t % results(RESULT_VALUE, score_index, filter_index) = & + t % results(RESULT_VALUE, score_index, filter_index) + score end select @@ -2458,8 +2457,8 @@ contains ! Add score to tally !$omp atomic - t % results(i_score, i_filter) % value = & - t % results(i_score, i_filter) % value + score * filter_weight + t % results(RESULT_VALUE, i_score, i_filter) = & + t % results(RESULT_VALUE, i_score, i_filter) + score * filter_weight ! Case for tallying delayed emissions else if (score_bin == SCORE_DELAYED_NU_FISSION .and. g /= 0) then @@ -2498,8 +2497,8 @@ contains ! Add score to tally !$omp atomic - t % results(i_score, i_filter) % value = & - t % results(i_score, i_filter) % value + score * filter_weight + t % results(RESULT_VALUE, i_score, i_filter) = & + t % results(RESULT_VALUE, i_score, i_filter) + score * filter_weight end if end if end do @@ -2536,8 +2535,8 @@ contains filter_weight = product(filter_weights(:size(t % filters))) !$omp atomic - t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) % value + score * filter_weight + t % results(RESULT_VALUE, score_index, filter_index) = & + t % results(RESULT_VALUE, score_index, filter_index) + score * filter_weight ! reset original delayed group bin matching_bins(t % find_filter(FILTER_DELAYEDGROUP)) = bin_original @@ -2997,8 +2996,8 @@ contains filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt + t % results(RESULT_VALUE, 1, filter_index) = & + t % results(RESULT_VALUE, 1, filter_index) + p % wgt end if ! Inward current on d1 min surface @@ -3033,8 +3032,8 @@ contains filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt + t % results(RESULT_VALUE, 1, filter_index) = & + t % results(RESULT_VALUE, 1, filter_index) + p % wgt ijk0(d1) = ijk0(d1) - 1 end if @@ -3053,8 +3052,8 @@ contains filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt + t % results(RESULT_VALUE, 1, filter_index) = & + t % results(RESULT_VALUE, 1, filter_index) + p % wgt end if ! Inward current on d1 max surface @@ -3089,8 +3088,8 @@ contains filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt + t % results(RESULT_VALUE, 1, filter_index) = & + t % results(RESULT_VALUE, 1, filter_index) + p % wgt ijk0(d1) = ijk0(d1) + 1 end if @@ -3116,9 +3115,10 @@ contains subroutine synchronize_tallies() integer :: i - real(8) :: k_col ! Copy of batch collision estimate of keff - real(8) :: k_abs ! Copy of batch absorption estimate of keff - real(8) :: k_tra ! Copy of batch tracklength estimate of keff + real(C_DOUBLE) :: k_col ! Copy of batch collision estimate of keff + real(C_DOUBLE) :: k_abs ! Copy of batch absorption estimate of keff + real(C_DOUBLE) :: k_tra ! Copy of batch tracklength estimate of keff + real(C_DOUBLE) :: val #ifdef MPI ! Combine tally results onto master process @@ -3142,9 +3142,9 @@ contains if (run_mode == MODE_EIGENVALUE) then if (active_batches) then ! Accumulate products of different estimators of k - k_col = global_tallies(K_COLLISION) % value / total_weight - k_abs = global_tallies(K_ABSORPTION) % value / total_weight - k_tra = global_tallies(K_TRACKLENGTH) % value / total_weight + k_col = global_tallies(RESULT_VALUE, K_COLLISION) / total_weight + k_abs = global_tallies(RESULT_VALUE, K_ABSORPTION) / total_weight + k_tra = global_tallies(RESULT_VALUE, K_TRACKLENGTH) / total_weight k_col_abs = k_col_abs + k_col * k_abs k_col_tra = k_col_tra + k_col * k_tra k_abs_tra = k_abs_tra + k_abs * k_tra @@ -3152,7 +3152,14 @@ contains end if ! Accumulate results for global tallies - call accumulate_result(global_tallies) + do i = 1, size(global_tallies, 2) + val = global_tallies(RESULT_VALUE, i)/total_weight + global_tallies(RESULT_VALUE, i) = ZERO + + global_tallies(RESULT_SUM, i) = global_tallies(RESULT_SUM, i) + val + global_tallies(RESULT_SUM_SQ, i) = & + global_tallies(RESULT_SUM_SQ, i) + val*val + end do end if end subroutine synchronize_tallies @@ -3182,7 +3189,7 @@ contains allocate(tally_temp(m,n)) - tally_temp = t % results(:,:) % value + tally_temp = t % results(RESULT_VALUE,:,:) if (master) then ! The MPI_IN_PLACE specifier allows the master to copy values into @@ -3191,35 +3198,35 @@ contains MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Transfer values to value on master - t % results(:,:) % value = tally_temp + t % results(RESULT_VALUE,:,:) = tally_temp else ! Receive buffer not significant at other processors call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, & MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Reset value on other processors - t % results(:,:) % value = 0 + t % results(RESULT_VALUE,:,:) = ZERO end if deallocate(tally_temp) end do ! Copy global tallies into array to be reduced - global_temp = global_tallies(:) % value + global_temp = global_tallies(RESULT_VALUE, :) if (master) then call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, & MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Transfer values back to global_tallies on master - global_tallies(:) % value = global_temp + global_tallies(RESULT_VALUE, :) = global_temp else ! Receive buffer not significant at other processors call MPI_REDUCE(global_temp, dummy, N_GLOBAL_TALLIES, & MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Reset value on other processors - global_tallies(:) % value = ZERO + global_tallies(RESULT_VALUE, :) = ZERO end if ! We also need to determine the total starting weight of particles from the @@ -3244,6 +3251,9 @@ contains type(TallyObject), intent(inout) :: t + integer :: i, j + real(C_DOUBLE) :: val + ! Increment number of realizations if (reduce_tallies) then t % n_realizations = t % n_realizations + 1 @@ -3251,92 +3261,59 @@ contains t % n_realizations = t % n_realizations + n_procs end if - ! Accumulate each TallyResult - call accumulate_result(t % results) + ! Accumulate each result + do j = 1, size(t % results, 3) + do i = 1, size(t % results, 2) + val = t % results(RESULT_VALUE, i, j)/total_weight + t % results(RESULT_VALUE, i, j) = ZERO + + t % results(RESULT_SUM, i, j) = & + t % results(RESULT_SUM, i, j) + val + t % results(RESULT_SUM_SQ, i, j) = & + t % results(RESULT_SUM_SQ, i, j) + val*val + end do + end do end subroutine accumulate_tally !=============================================================================== ! TALLY_STATISTICS computes the mean and standard deviation of the mean of each -! tally and stores them in the val and val_sq attributes of the TallyResults -! respectively +! tally and stores them in the RESULT_SUM and RESULT_SUM_SQ positions !=============================================================================== subroutine tally_statistics() - integer :: i ! index in tallies array - type(TallyObject), pointer :: t - - ! Calculate statistics for user-defined tallies - do i = 1, n_tallies - t => tallies(i) - - call statistics_result(t % results, t % n_realizations) - end do - - ! Calculate statistics for global tallies - call statistics_result(global_tallies, n_realizations) - - end subroutine tally_statistics - -!=============================================================================== -! ACCUMULATE_RESULT accumulates results from many histories (or many generations) -! into a single realization of a random variable. -!=============================================================================== - - elemental subroutine accumulate_result(this) - - type(TallyResult), intent(inout) :: this - - real(8) :: val - - ! Add the sum and square of the sum of contributions from a tally result to - ! the variables sum and sum_sq. This will later allow us to calculate a - ! variance on the tallies. - - val = this % value/total_weight - this % sum = this % sum + val - this % sum_sq = this % sum_sq + val*val - - ! Reset the single batch estimate - this % value = ZERO - - end subroutine accumulate_result - -!=============================================================================== -! STATISTICS_RESULT determines the sample mean and the standard deviation of the -! mean for a TallyResult. -!=============================================================================== - - elemental subroutine statistics_result(this, n) - - type(TallyResult), intent(inout) :: this - integer, intent(in) :: n + integer :: j, k ! score/filter indices + integer :: n ! number of realizations ! Calculate sample mean and standard deviation of the mean -- note that we ! have used Bessel's correction so that the estimator of the variance of the ! sample mean is unbiased. - this % sum = this % sum/n - this % sum_sq = sqrt((this % sum_sq/n - this % sum * & - this % sum) / (n - 1)) + do i = 1, n_tallies + n = tallies(i) % n_realizations - end subroutine statistics_result + associate (r => tallies(i) % results) + do k = 1, size(r, 3) + do j = 1, size(r, 2) + r(RESULT_SUM, j, k) = r(RESULT_SUM, j, k) / n + r(RESULT_SUM_SQ, j, k) = sqrt((r(RESULT_SUM_SQ, j, k)/n - & + r(RESULT_SUM, j, k) * r(RESULT_SUM, j, k))/(n - 1)) + end do + end do + end associate + end do -!=============================================================================== -! RESET_RESULT zeroes out the value and accumulated sum and sum-squared for a -! single TallyResult. -!=============================================================================== - - elemental subroutine reset_result(this) - - type(TallyResult), intent(inout) :: this - - this % value = ZERO - this % sum = ZERO - this % sum_sq = ZERO - - end subroutine reset_result + ! Calculate statistics for global tallies + n = n_realizations + associate (r => global_tallies) + do j = 1, size(r, 2) + r(RESULT_SUM, j) = r(RESULT_SUM, j) / n + r(RESULT_SUM_SQ, j) = sqrt((r(RESULT_SUM_SQ, j)/n - & + r(RESULT_SUM, j) * r(RESULT_SUM, j))/(n - 1)) + end do + end associate + end subroutine tally_statistics !=============================================================================== ! SETUP_ACTIVE_USERTALLIES diff --git a/src/tally_header.F90 b/src/tally_header.F90 index fe385458b2..eeda03e762 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -1,23 +1,15 @@ module tally_header + use, intrinsic :: ISO_C_BINDING + + use hdf5 + use constants, only: NONE, N_FILTER_TYPES use tally_filter_header, only: TallyFilterContainer use trigger_header, only: TriggerObject - use, intrinsic :: ISO_C_BINDING - implicit none -!=============================================================================== -! TALLYRESULT provides accumulation of results in a particular tally bin -!=============================================================================== - - type, bind(C) :: TallyResult - real(C_DOUBLE) :: value = 0. - real(C_DOUBLE) :: sum = 0. - real(C_DOUBLE) :: sum_sq = 0. - end type TallyResult - !=============================================================================== ! TALLYOBJECT describes a user-specified tally. The region of phase space to ! tally in is given by the TallyFilters and the results are stored in a @@ -68,7 +60,7 @@ module tally_header integer :: total_filter_bins integer :: total_score_bins - type(TallyResult), allocatable :: results(:,:) + real(C_DOUBLE), allocatable :: results(:,:,:) ! reset property - allows a tally to be reset after every batch logical :: reset = .false. @@ -79,6 +71,79 @@ module tally_header ! Tally precision triggers integer :: n_triggers = 0 ! # of triggers type(TriggerObject), allocatable :: triggers(:) ! Array of triggers + + contains + procedure :: write_results_hdf5 + procedure :: read_results_hdf5 end type TallyObject +contains + + subroutine write_results_hdf5(this, group_id) + class(TallyObject), intent(in) :: this + integer(HID_T), intent(in) :: group_id + + integer :: hdf5_err + integer(HID_T) :: dset, dspace + integer(HID_T) :: memspace + integer(HSIZE_T) :: dims(3) + integer(HSIZE_T) :: dims_slab(3) + integer(HSIZE_T) :: offset(3) = [1,0,0] + + ! Create file dataspace + dims_slab(:) = shape(this % results) + dims_slab(1) = 2 + call h5screate_simple_f(3, dims_slab, dspace, hdf5_err) + + ! Create memory dataspace that contains only SUM and SUM_SQ values + dims(:) = shape(this % results) + call h5screate_simple_f(3, dims, memspace, hdf5_err) + call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offset, dims_slab, & + hdf5_err) + + ! Create and write to dataset + call h5dcreate_f(group_id, "results", H5T_NATIVE_DOUBLE, dspace, dset, & + hdf5_err) + call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, this % results, dims_slab, & + hdf5_err, mem_space_id=memspace) + + ! Close identifiers + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(memspace, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + end subroutine write_results_hdf5 + + subroutine read_results_hdf5(this, group_id) + class(TallyObject), intent(inout) :: this + integer(HID_T), intent(in) :: group_id + + integer :: hdf5_err + integer(HID_T) :: dset, dspace + integer(HID_T) :: memspace + integer(HSIZE_T) :: dims(3) + integer(HSIZE_T) :: dims_slab(3) + integer(HSIZE_T) :: offset(3) = [1,0,0] + + ! Create file dataspace + dims_slab(:) = shape(this % results) + dims_slab(1) = 2 + call h5screate_simple_f(3, dims_slab, dspace, hdf5_err) + + ! Create memory dataspace that contains only SUM and SUM_SQ values + dims(:) = shape(this % results) + call h5screate_simple_f(3, dims, memspace, hdf5_err) + call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offset, dims_slab, & + hdf5_err) + + ! Create and write to dataset + call h5dopen_f(group_id, "results", dset, hdf5_err) + call h5dread_f(dset, H5T_NATIVE_DOUBLE, this % results, dims_slab, & + hdf5_err, mem_space_id=memspace) + + ! Close identifiers + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(memspace, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + end subroutine read_results_hdf5 + end module tally_header diff --git a/src/tally_initialize.F90 b/src/tally_initialize.F90 index 1ae403bd65..0e7a50a9c5 100644 --- a/src/tally_initialize.F90 +++ b/src/tally_initialize.F90 @@ -20,7 +20,8 @@ contains subroutine configure_tallies() ! Allocate global tallies - allocate(global_tallies(N_GLOBAL_TALLIES)) + allocate(global_tallies(3, N_GLOBAL_TALLIES)) + global_tallies(:,:) = ZERO call setup_tally_arrays() @@ -62,7 +63,8 @@ contains t % total_score_bins = t % n_score_bins * t % n_nuclide_bins ! Allocate results array - allocate(t % results(t % total_score_bins, t % total_filter_bins)) + allocate(t % results(3, t % total_score_bins, t % total_filter_bins)) + t % results(:,:,:) = ZERO end do TALLY_LOOP diff --git a/src/trigger.F90 b/src/trigger.F90 index 39786809b4..f41cc46cf3 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -432,17 +432,19 @@ contains real(8), intent(inout) :: rel_err ! tally relative error integer, intent(in) :: score_index ! tally results score index integer, intent(in) :: filter_index ! tally results filter index - integer :: n ! number of realizations - real(8) :: mean ! tally mean - type(TallyResult) :: tally_result ! pointer to TallyResult - type(TallyObject), pointer :: t ! tally pointer + type(TallyObject), intent(in) :: t ! tally + + integer :: n ! number of realizations + real(8) :: mean ! tally mean + real(8) :: tally_sum, tally_sum_sq ! results for a single tally bin n = t % n_realizations - tally_result = t % results(score_index, filter_index) + tally_sum = t % results(RESULT_SUM, score_index, filter_index) + tally_sum_sq = t % results(RESULT_SUM_SQ, score_index, filter_index) ! Compute the tally mean and standard deviation - mean = tally_result % sum / n - std_dev = sqrt((tally_result % sum_sq / n - mean * mean) / (n - 1)) + mean = tally_sum / n + std_dev = sqrt((tally_sum_sq / n - mean * mean) / (n - 1)) ! Compute the relative error if the mean is non-zero if (mean == ZERO) then From 0bd1ddd569f944b300549911ea9235f52abcae95 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 16 Nov 2016 09:16:03 -0600 Subject: [PATCH 19/19] Add back 'end subroutine' that was accidentally deleted --- src/initialize.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/initialize.F90 b/src/initialize.F90 index d51f9ef8eb..692319acb8 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -211,6 +211,8 @@ contains call MPI_TYPE_CREATE_STRUCT(5, bank_blocks, bank_disp, & bank_types, MPI_BANK, mpi_err) call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) + + end subroutine initialize_mpi #endif !===============================================================================