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/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 8452046b3b..5f34d4700c 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -240,12 +240,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/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..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 @@ -149,6 +150,7 @@ def atomic_mass(isotope): """ if not _ATOMIC_MASS: + # Load data from AME2012 file mass_file = os.path.join(os.path.dirname(__file__), 'mass.mas12') with open(mass_file, 'r') as ame: @@ -159,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/data/function.py b/openmc/data/function.py index e9ecf82c58..0515a57c0c 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -397,6 +397,67 @@ 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 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 + combinations of Combination objects. The operations parameter must have + a length one less than the number of functions. + + + Attributes + ---------- + functions : Iterable of Callable + 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 + combinations of Combination objects. The operations parameter must have + a length one less than the number of functions. + + """ + + def __init__(self, functions, operations): + self.functions = functions + self.operations = operations + + def __call__(self, x): + ans = self.functions[0](x) + for i, operation in enumerate(self.operations): + ans = operation(ans, self.functions[i + 1](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) - 1 + cv.check_length('operations', operations, length, length_max=length) + self._operations = operations + + class Sum(EqualityMixin): """Sum of multiple functions. @@ -432,6 +493,63 @@ class Sum(EqualityMixin): self._functions = functions +class Regions1D(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 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 + ---------- + 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. diff --git a/openmc/data/library.py b/openmc/data/library.py index 0485af0641..9f6d931467 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -22,6 +22,21 @@ class DataLibrary(EqualityMixin): def __init__(self): self.libraries = [] + def get_by_material(self, value): + """Return the library dictionary containing a given material. + + Returns + ------- + library : dict or None + Dictionary summarizing cross section data from a single file; + the dictionary has keys 'path', 'type', and 'materials'. + + """ + for library in self.libraries: + if value in library['materials']: + return library + return None + def register_file(self, filename): """Register a file with the data library. @@ -78,3 +93,38 @@ 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. + + Returns + ------- + data : openmc.data.DataLibrary + Data library object initialized from the provided XML + + """ + + 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 diff --git a/openmc/material.py b/openmc/material.py index c265c90694..70a070ddd6 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -64,6 +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'). + 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. """ @@ -194,6 +198,27 @@ class Material(object): def distrib_otf_file(self): return self._distrib_otf_file + @property + def average_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): 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 5a677a8e2f..144146ba3a 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 @@ -107,12 +108,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 @@ -213,13 +208,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 1f02129c0d..d1789b78e9 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 @@ -280,20 +281,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) @@ -1763,8 +1758,6 @@ class Tally(object): # HDF5 binary file if format == 'hdf5': - import h5py - filename = directory + '/' + filename + '.h5' if append: diff --git a/openmc/universe.py b/openmc/universe.py index 545b05d5a1..b2e7dbb920 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, + **kwargs): """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. + **kwargs + All keyword arguments are passed to + :func:`matplotlib.pyplot.imshow`. """ 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), **kwargs) # Show or save the plot if filename is None: 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 93c01c0f7b..044ea1c6bc 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 9f66cc1a1f..02f6870189 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, TallyDerivative + use tally_header, only: TallyObject, TallyDerivative 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 @@ -276,10 +278,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 c8892c53c4..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 @@ -2589,7 +2459,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/initialize.F90 b/src/initialize.F90 index 4fe6f8769c..692319acb8 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. @@ -222,34 +212,6 @@ contains 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 +221,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 +229,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 c88f8a6805..1bbf68f48b 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 123cc9caa7..eac4ea94cc 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/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(:, :, :) diff --git a/src/output.F90 b/src/output.F90 index 883165556c..f5766187ee 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 @@ -898,8 +898,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) @@ -908,9 +908,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, & @@ -924,9 +923,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 @@ -938,8 +937,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 @@ -1042,16 +1041,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 @@ -1059,16 +1058,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 @@ -1078,16 +1077,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 @@ -1095,16 +1094,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 @@ -1114,16 +1113,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 @@ -1131,16 +1130,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 93950b1f9d..e6b20ddf8f 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 cf334787a6..23543775cb 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -393,7 +393,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 @@ -521,7 +521,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 @@ -529,7 +529,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 @@ -545,9 +545,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 @@ -559,20 +558,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 @@ -608,15 +598,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 " // & @@ -632,20 +622,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 @@ -811,7 +801,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") @@ -827,7 +817,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 df8efa0c02..79f22610e0 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 @@ -19,7 +21,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 @@ -1974,8 +1975,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) @@ -1991,10 +1992,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 @@ -2020,10 +2020,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 @@ -2040,8 +2039,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) @@ -2049,8 +2048,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 @@ -2474,8 +2473,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 @@ -2514,8 +2513,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 @@ -2552,8 +2551,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 @@ -3013,8 +3012,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 @@ -3049,8 +3048,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 @@ -3069,8 +3068,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 @@ -3105,8 +3104,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 @@ -3904,9 +3903,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 @@ -3930,9 +3930,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 @@ -3940,7 +3940,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 @@ -3970,7 +3977,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 @@ -3979,35 +3986,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 @@ -4032,6 +4039,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 @@ -4039,92 +4049,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 e292aba60b..1d34779aec 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 - !=============================================================================== ! TALLYDERIVATIVE describes a first-order derivative that can be applied to ! tallies. @@ -81,7 +73,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. @@ -95,6 +87,79 @@ module tally_header ! Index for the TallyDerivative for differential tallies. integer :: deriv = NONE + + 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