diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index e8b8f8b90..5f34d4700 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -162,6 +162,18 @@ if run_mode == 'k-eigenvalue': Width of each mesh cell in each dimension. +**/tallies/derivatives/derivative /independent variable** (*char[]*) + + Independent variable of tally derivative + +**/tallies/derivatives/derivative /material** (*int*) + + ID of the perturbed material + +**/tallies/derivatives/derivative /nuclide** (*char[]*) + + Alias of the perturbed nuclide + **/tallies/n_tallies** (*int*) Number of user-defined tallies. @@ -204,6 +216,10 @@ if run_mode == 'k-eigenvalue': Array of nuclides to tally. Note that if no nuclide is specified in the user input, a single 'total' nuclide appears here. +**/tallies/tally /derivative** (*int*) + + ID of the derivative applied to the tally. + **/tallies/tally /n_score_bins** (*int*) Number of scoring bins for a single nuclide. In general, this can be greater @@ -224,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/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 7e310857c..e364452b1 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -293,6 +293,7 @@ Multi-group Cross Sections openmc.mgxs.NuScatterXS openmc.mgxs.NuScatterMatrixXS openmc.mgxs.PromptNuFissionXS + openmc.mgxs.PromptNuFissionMatrixXS openmc.mgxs.ScatterXS openmc.mgxs.ScatterMatrixXS openmc.mgxs.TotalXS @@ -309,6 +310,7 @@ Multi-delayed-group Cross Sections openmc.mgxs.MDGXS openmc.mgxs.ChiDelayed openmc.mgxs.DelayedNuFissionXS + openmc.mgxs.DelayedNuFissionMatrixXS openmc.mgxs.Beta openmc.mgxs.DecayRate diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 92f53e6bf..4a81615a1 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1916,6 +1916,13 @@ The ```` element accepts the following sub-elements: *Default*: "all" + :derivative: + The id of a ``derivative`` element. This derivative will be applied to all + scores in the tally. Differential tallies are currently only implemented + for collision and analog estimators. + + *Default*: None + ```` Element ------------------ @@ -1944,6 +1951,39 @@ attributes/sub-elements: One of ```` or ```` must be specified, but not both (even if they are consistent with one another). +```` Element +------------------------ + +OpenMC can take the first-order derivative of many tallies with respect to +material perturbations. It works by propagating a derivative through the +transport equation. Essentially, OpenMC keeps track of how each particle's +weight would change as materials are perturbed, and then accounts for that +weight change in the tallies. Note that this assumes material perturbations are +small enough not to change the distribution of fission sites. This element has +the following attributes/sub-elements: + + :id: + A unique integer that can be used to identify the derivative. + + :variable: + The independent variable of the derivative. Accepted options are "density", + "nuclide_density", and "temperature". A "density" derivative will give the + derivative with respect to the density of the material in [g / cm^3]. A + "nuclide_density" derivative will give the derivative with respect to the + density of a particular nuclide in units of [atom / b / cm]. A + "temperature" derivative is with respect to a material temperature in units + of [K]. The temperature derivative requires windowed multipole to be + turned on. Note also that the temperature derivative only accounts for + resolved resonance Doppler broadening. It does not account for thermal + expansion, S(a, b) scattering, resonance scattering, or unresolved Doppler + broadening. + + :material: + The perturbed material. (Necessary for all derivative types) + + :nuclide: + The perturbed nuclide. (Necessary only for "nuclide_density") + ```` Element ----------------------------- diff --git a/openmc/__init__.py b/openmc/__init__.py index 0685a80b1..40a4dcae4 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -16,6 +16,7 @@ from openmc.universe import * from openmc.mesh import * from openmc.filter import * from openmc.trigger import * +from openmc.tally_derivative import * from openmc.tallies import * from openmc.mgxs_library import * from openmc.cmfd import * diff --git a/openmc/data/data.py b/openmc/data/data.py index a3e715812..92be7cade 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -186,3 +186,9 @@ K_BOLTZMANN = 8.6173324e-5 # Used for converting units in ACE data EV_PER_MEV = 1.0e6 + +# Avogadro's constant from CODATA 2010 +AVOGADRO = 6.02214129E23 + +# Neutron mass from CODATA 2010 in units of amu +NEUTRON_MASS = 1.008664916 diff --git a/openmc/data/library.py b/openmc/data/library.py index 3beb25543..c179f78f8 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -1,5 +1,6 @@ import os import xml.etree.ElementTree as ET +from six import string_types import h5py @@ -114,7 +115,7 @@ class DataLibrary(EqualityMixin): data = cls() - # If cross_sections is None, get the cross sections from the + # If path is None, get the cross sections from the # OPENMC_CROSS_SECTIONS environment variable if path is None: path = os.environ.get('OPENMC_CROSS_SECTIONS') @@ -124,7 +125,7 @@ class DataLibrary(EqualityMixin): raise ValueError("Either path or OPENMC_CROSS_SECTIONS " "environmental variable must be set") - check_type('path', path, str) + check_type('path', path, string_types) tree = ET.parse(path) root = tree.getroot() diff --git a/openmc/material.py b/openmc/material.py index 16d1fd147..08b001e47 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -606,15 +606,10 @@ class Material(object): return nuclides - def get_nuclide_atom_densities(self, cross_sections=None): + def get_nuclide_atom_densities(self): """Returns all nuclides in the material and their atomic densities in units of atom/b-cm - Parameters - ---------- - cross_sections : str, optional - Location of cross_sections.xml file. Default is None. - Returns ------- nuclides : dict @@ -623,11 +618,6 @@ class Material(object): """ - import scipy.constants as sc - - # Load the library - library = openmc.data.DataLibrary.from_xml(cross_sections) - # Expand elements in to nuclides nuclides = self.get_nuclide_densities() @@ -664,17 +654,12 @@ class Material(object): sum_percent = 0. awrs = [] - n = -1 - for nuclide in nuclides.items(): - n += 1 - lib = library.get_by_material(nuclide[0]) - if lib is not None: - nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) - awrs.append(nuc.atomic_weight_ratio) - if not percent_in_atom: - nuc_densities[n] = -nuc_densities[n] / awrs[-1] + for n, nuclide in enumerate(nuclides.items()): + awr = openmc.data.atomic_mass(nuclide[0]) + if awr is not None: + awrs.append(awr / openmc.data.NEUTRON_MASS) else: - raise ValueError(nuclide[0] + " not in library") + raise ValueError(nuclide[0] + " is invalid") # Now that we have the awr, lets finish calculating densities sum_percent = np.sum(nuc_densities) @@ -686,7 +671,7 @@ class Material(object): sum_percent += x * awrs[n] sum_percent = 1. / sum_percent density = -density * sum_percent * \ - sc.Avogadro / sc.value('neutron mass in u') * 1.E-24 + openmc.data.AVOGADRO / openmc.data.NEUTRON_MASS * 1.E-24 nuc_densities = density * nuc_densities nuclides = OrderedDict() diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 513884d96..2e2b7f908 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -582,7 +582,7 @@ class Library(object): ---------- domain : Material or Cell or Universe or Integral The material, cell, or universe object of interest (or its ID) - mgxs_type : {'total', 'transport', 'nu-transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'multiplicity matrix', 'nu-fission matrix', chi', 'chi-prompt', 'inverse-velocity', 'prompt-nu-fission', 'delayed-nu-fission', 'chi-delayed', 'beta'} + mgxs_type : {'total', 'transport', 'nu-transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'multiplicity matrix', 'nu-fission matrix', chi', 'chi-prompt', 'inverse-velocity', 'prompt-nu-fission', 'prompt-nu-fission matrix', 'delayed-nu-fission', 'delayed-nu-fission matrix', 'chi-delayed', 'beta'} The type of multi-group cross section object to return Returns @@ -1026,12 +1026,24 @@ class Library(object): nuclide=[nuclide], subdomain=subdomain) + if 'prompt-nu-fission matrix' in self.mgxs_types: + mymgxs = self.get_mgxs(domain, 'prompt-nu-fission matrix') + xsdata.set_prompt_nu_fission_mgxs(mymgxs, xs_type=xs_type, + nuclide=[nuclide], + subdomain=subdomain) + if 'delayed-nu-fission' in self.mgxs_types: mymgxs = self.get_mgxs(domain, 'delayed-nu-fission') xsdata.set_delayed_nu_fission_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], subdomain=subdomain) + if 'delayed-nu-fission matrix' in self.mgxs_types: + mymgxs = self.get_mgxs(domain, 'delayed-nu-fission matrix') + xsdata.set_delayed_nu_fission_mgxs(mymgxs, xs_type=xs_type, + nuclide=[nuclide], + subdomain=subdomain) + if 'beta' in self.mgxs_types: mymgxs = self.get_mgxs(domain, 'nu-fission') xsdata.set_beta_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], diff --git a/openmc/mgxs/mdgxs.py b/openmc/mgxs/mdgxs.py index a99a99bee..5b2f451c2 100644 --- a/openmc/mgxs/mdgxs.py +++ b/openmc/mgxs/mdgxs.py @@ -21,7 +21,8 @@ import openmc.checkvalue as cv MDGXS_TYPES = ['delayed-nu-fission', 'chi-delayed', 'beta', - 'decay-rate'] + 'decay-rate', + 'delayed-nu-fission matrix'] # Maximum number of delayed groups, from src/constants.F90 MAX_DELAYED_GROUPS = 8 @@ -211,7 +212,7 @@ class MDGXS(MGXS): Parameters ---------- - mdgxs_type : {'delayed-nu-fission', 'chi-delayed', 'beta', 'decay-rate'} + mdgxs_type : {'delayed-nu-fission', 'chi-delayed', 'beta', 'decay-rate', 'delayed-nu-fission matrix'} The type of multi-delayed-group cross section object to return domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh @@ -249,6 +250,9 @@ class MDGXS(MGXS): mdgxs = Beta(domain, domain_type, energy_groups, delayed_groups) elif mdgxs_type == 'decay-rate': mdgxs = DecayRate(domain, domain_type, energy_groups, delayed_groups) + elif mdgxs_type == 'delayed-nu-fission matrix': + mdgxs = DelayedNuFissionMatrixXS(domain, domain_type, energy_groups, + delayed_groups) mdgxs.by_nuclide = by_nuclide mdgxs.name = name @@ -1733,3 +1737,609 @@ class DecayRate(MDGXS): super(DecayRate, self)._compute_xs() return self._xs_tally + + +@add_metaclass(ABCMeta) +class MatrixMDGXS(MDGXS): + """An abstract multi-delayed-group cross section for some energy group and + delayed group structure within some spatial domain. This class is + specifically intended for cross sections which depend on both the incoming + and outgoing energy groups and are therefore represented by matrices. + An example of this is the delayed-nu-fission matrix. + + This class can be used for both OpenMC input generation and tally data + post-processing to compute spatially-homogenized and energy-integrated + multi-group and multi-delayed-group cross sections for downstream neutronics + calculations. + + NOTE: Users should instantiate the subclasses of this abstract class. + + Parameters + ---------- + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + The domain for spatial homogenization + domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} + The domain type for spatial homogenization + energy_groups : openmc.mgxs.EnergyGroups + The energy group structure for energy condensation + by_nuclide : bool + If true, computes cross sections for each nuclide in domain + name : str, optional + Name of the multi-group cross section. Used as a label to identify + tallies in OpenMC 'tallies.xml' file. + delayed_groups : list of int + Delayed groups to filter out the xs + + Attributes + ---------- + name : str, optional + Name of the multi-group cross section + rxn_type : str + Reaction type (e.g., 'total', 'nu-fission', etc.) + by_nuclide : bool + If true, computes cross sections for each nuclide in domain + domain : Material or Cell or Universe or Mesh + Domain for spatial homogenization + domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} + Domain type for spatial homogenization + energy_groups : openmc.mgxs.EnergyGroups + Energy group structure for energy condensation + delayed_groups : list of int + Delayed groups to filter out the xs + tally_trigger : openmc.Trigger + An (optional) tally precision trigger given to each tally used to + compute the cross section + scores : list of str + The scores in each tally used to compute the multi-group cross section + filters : list of openmc.Filter + The filters in each tally used to compute the multi-group cross section + tally_keys : list of str + The keys into the tallies dictionary for each tally used to compute + the multi-group cross section + estimator : {'tracklength', 'collision', 'analog'} + The tally estimator used to compute the multi-group cross section + tallies : collections.OrderedDict + OpenMC tallies needed to compute the multi-group cross section + rxn_rate_tally : openmc.Tally + Derived tally for the reaction rate tally used in the numerator to + compute the multi-group cross section. This attribute is None + unless the multi-group cross section has been computed. + xs_tally : openmc.Tally + Derived tally for the multi-group cross section. This attribute + is None unless the multi-group cross section has been computed. + num_subdomains : int + The number of subdomains is unity for 'material', 'cell' and 'universe' + domain types. This is equal to the number of cell instances + for 'distribcell' domain types (it is equal to unity prior to loading + tally data from a statepoint file) and the number of mesh cells for + 'mesh' domain types. + num_nuclides : int + The number of nuclides for which the multi-group cross section is + being tracked. This is unity if the by_nuclide attribute is False. + nuclides : Iterable of str or 'sum' + The optional user-specified nuclides for which to compute cross + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides + are not specified by the user, all nuclides in the spatial domain + are included. This attribute is 'sum' if by_nuclide is false. + sparse : bool + Whether or not the MGXS' tallies use SciPy's LIL sparse matrix format + for compressed data storage + loaded_sp : bool + Whether or not a statepoint file has been loaded with tally data + derived : bool + Whether or not the MGXS is merged from one or more other MGXS + hdf5_key : str + The key used to index multi-group cross sections in an HDF5 data store + + """ + + @property + def filters(self): + # Create the non-domain specific Filters for the Tallies + group_edges = self.energy_groups.group_edges + energy = openmc.EnergyFilter(group_edges) + energyout = openmc.EnergyoutFilter(group_edges) + + if self.delayed_groups is not None: + delayed = openmc.DelayedGroupFilter(self.delayed_groups) + return [[energy], [delayed, energy, energyout]] + else: + return [[energy], [energy, energyout]] + + def get_xs(self, in_groups='all', out_groups='all', + subdomains='all', nuclides='all', + xs_type='macro', order_groups='increasing', + row_column='inout', value='mean', delayed_groups='all', + squeeze=True, **kwargs): + """Returns an array of multi-group cross sections. + + This method constructs a 4D NumPy array for the requested + multi-group cross section data for one or more subdomains + (1st dimension), delayed groups (2nd dimension), energy groups in + (3rd dimension), energy groups out (4th dimension), and nuclides + (5th dimension). + + Parameters + ---------- + in_groups : Iterable of Integral or 'all' + Incoming energy groups of interest. Defaults to 'all'. + out_groups : Iterable of Integral or 'all' + Outgoing energy groups of interest. Defaults to 'all'. + subdomains : Iterable of Integral or 'all' + Subdomain IDs of interest. Defaults to 'all'. + nuclides : Iterable of str or 'all' or 'sum' + A list of nuclide name strings (e.g., ['U235', 'U238']). The + special string 'all' will return the cross sections for all + nuclides in the spatial domain. The special string 'sum' will + return the cross section summed over all nuclides. Defaults to + 'all'. + xs_type: {'macro', 'micro'} + Return the macro or micro cross section in units of cm^-1 or barns. + Defaults to 'macro'. + order_groups: {'increasing', 'decreasing'} + Return the cross section indexed according to increasing or + decreasing energy groups (decreasing or increasing energies). + Defaults to 'increasing'. + row_column: {'inout', 'outin'} + Return the cross section indexed first by incoming group and + second by outgoing group ('inout'), or vice versa ('outin'). + Defaults to 'inout'. + value : {'mean', 'std_dev', 'rel_err'} + A string for the type of value to return. Defaults to 'mean'. + delayed_groups : list of int or 'all' + Delayed groups of interest. Defaults to 'all'. + squeeze : bool + A boolean representing whether to eliminate the extra dimensions + of the multi-dimensional array to be returned. Defaults to True. + + Returns + ------- + numpy.ndarray + A NumPy array of the multi-group cross section indexed in the order + each group and subdomain is listed in the parameters. + + Raises + ------ + ValueError + When this method is called before the multi-group cross section is + computed from tally data. + + """ + + cv.check_value('value', value, ['mean', 'std_dev', 'rel_err']) + cv.check_value('xs_type', xs_type, ['macro', 'micro']) + + # FIXME: Unable to get microscopic xs for mesh domain because the mesh + # cells do not know the nuclide densities in each mesh cell. + if self.domain_type == 'mesh' and xs_type == 'micro': + msg = 'Unable to get micro xs for mesh domain since the mesh ' \ + 'cells do not know the nuclide densities in each mesh cell.' + raise ValueError(msg) + + filters = [] + filter_bins = [] + + # Construct a collection of the domain filter bins + if not isinstance(subdomains, string_types): + cv.check_iterable_type('subdomains', subdomains, Integral, + max_depth=3) + for subdomain in subdomains: + filters.append(_DOMAIN_TO_FILTER[self.domain_type]) + filter_bins.append((subdomain,)) + + # Construct list of energy group bounds tuples for all requested groups + if not isinstance(in_groups, string_types): + cv.check_iterable_type('groups', in_groups, Integral) + for group in in_groups: + filters.append(openmc.EnergyFilter) + filter_bins.append(( + self.energy_groups.get_group_bounds(group),)) + + # Construct list of energy group bounds tuples for all requested groups + if not isinstance(out_groups, string_types): + cv.check_iterable_type('groups', out_groups, Integral) + for group in out_groups: + filters.append(openmc.EnergyoutFilter) + filter_bins.append(( + self.energy_groups.get_group_bounds(group),)) + + # Construct list of delayed group tuples for all requested groups + if not isinstance(delayed_groups, string_types): + cv.check_type('delayed groups', delayed_groups, list, int) + for delayed_group in delayed_groups: + filters.append(openmc.DelayedGroupFilter) + filter_bins.append((delayed_group,)) + + # Construct a collection of the nuclides to retrieve from the xs tally + if self.by_nuclide: + if nuclides == 'all' or nuclides == 'sum' or nuclides == ['sum']: + query_nuclides = self.get_nuclides() + else: + query_nuclides = nuclides + else: + query_nuclides = ['total'] + + # Use tally summation if user requested the sum for all nuclides + if nuclides == 'sum' or nuclides == ['sum']: + xs_tally = self.xs_tally.summation(nuclides=query_nuclides) + xs = xs_tally.get_values(filters=filters, filter_bins=filter_bins, + value=value) + else: + xs = self.xs_tally.get_values(filters=filters, + filter_bins=filter_bins, + nuclides=query_nuclides, value=value) + + # Divide by atom number densities for microscopic cross sections + if xs_type == 'micro': + if self.by_nuclide: + densities = self.get_nuclide_densities(nuclides) + else: + densities = self.get_nuclide_densities('sum') + if value == 'mean' or value == 'std_dev': + xs /= densities[np.newaxis, :, np.newaxis] + + # Eliminate the trivial score dimension + xs = np.squeeze(xs, axis=len(xs.shape) - 1) + xs = np.nan_to_num(xs) + + if in_groups == 'all': + num_in_groups = self.num_groups + else: + num_in_groups = len(in_groups) + + if out_groups == 'all': + num_out_groups = self.num_groups + else: + num_out_groups = len(out_groups) + + if delayed_groups == 'all': + num_delayed_groups = self.num_delayed_groups + else: + num_delayed_groups = len(delayed_groups) + + # Reshape tally data array with separate axes for domain and energy + num_subdomains = int(xs.shape[0] / (num_in_groups * num_out_groups * + num_delayed_groups)) + new_shape = (num_subdomains, num_delayed_groups, num_in_groups, + num_out_groups) + new_shape += xs.shape[1:] + xs = np.reshape(xs, new_shape) + + # Transpose the matrix if requested by user + if row_column == 'outin': + xs = np.swapaxes(xs, 2, 3) + + # Reverse data if user requested increasing energy groups since + # tally data is stored in order of increasing energies + if order_groups == 'increasing': + xs = xs[:, :, ::-1, ::-1, :] + + if squeeze: + xs = np.squeeze(xs) + xs = np.atleast_2d(xs) + + return xs + + def get_slice(self, nuclides=[], in_groups=[], out_groups=[], + delayed_groups=[]): + """Build a sliced MatrixMDGXS object for the specified nuclides and + energy groups. + + This method constructs a new MdGXS to encapsulate a subset of the data + represented by this MdGXS. The subset of data to include in the tally + slice is determined by the nuclides, energy groups, and delayed groups + specified in the input parameters. + + Parameters + ---------- + nuclides : list of str + A list of nuclide name strings + (e.g., ['U235', 'U238']; default is []) + in_groups : list of int + A list of incoming energy group indices starting at 1 for the high + energies (e.g., [1, 2, 3]; default is []) + out_groups : list of int + A list of outgoing energy group indices starting at 1 for the high + energies (e.g., [1, 2, 3]; default is []) + delayed_groups : list of int + A list of delayed group indices + (e.g., [1, 2, 3]; default is []) + + Returns + ------- + openmc.mgxs.MatrixMDGXS + A new MatrixMDGXS object which encapsulates the subset of data + requested for the nuclide(s) and/or energy group(s) requested in + the parameters. + + """ + + # Call super class method and null out derived tallies + slice_xs = super(MatrixMDGXS, self).get_slice(nuclides, in_groups, + delayed_groups) + slice_xs._rxn_rate_tally = None + slice_xs._xs_tally = None + + # Slice outgoing energy groups if needed + if len(out_groups) != 0: + filter_bins = [] + for group in out_groups: + group_bounds = self.energy_groups.get_group_bounds(group) + filter_bins.append(group_bounds) + filter_bins = [tuple(filter_bins)] + + # Slice each of the tallies across energyout groups + for tally_type, tally in slice_xs.tallies.items(): + if tally.contains_filter(openmc.EnergyoutFilter): + tally_slice = tally.get_slice( + filters=[openmc.EnergyoutFilter], + filter_bins=filter_bins) + slice_xs.tallies[tally_type] = tally_slice + + slice_xs.sparse = self.sparse + return slice_xs + + def print_xs(self, subdomains='all', nuclides='all', xs_type='macro'): + """Prints a string representation for the multi-group cross section. + + Parameters + ---------- + subdomains : Iterable of Integral or 'all' + The subdomain IDs of the cross sections to include in the report. + Defaults to 'all'. + nuclides : Iterable of str or 'all' or 'sum' + The nuclides of the cross-sections to include in the report. This + may be a list of nuclide name strings (e.g., ['U235', 'U238']). + The special string 'all' will report the cross sections for all + nuclides in the spatial domain. The special string 'sum' will + report the cross sections summed over all nuclides. Defaults to + 'all'. + xs_type: {'macro', 'micro'} + Return the macro or micro cross section in units of cm^-1 or barns. + Defaults to 'macro'. + + """ + + # Construct a collection of the subdomains to report + if not isinstance(subdomains, string_types): + cv.check_iterable_type('subdomains', subdomains, Integral) + elif self.domain_type == 'distribcell': + subdomains = np.arange(self.num_subdomains, dtype=np.int) + elif self.domain_type == 'mesh': + xyz = [range(1, x+1) for x in self.domain.dimension] + subdomains = list(itertools.product(*xyz)) + else: + subdomains = [self.domain.id] + + # Construct a collection of the nuclides to report + if self.by_nuclide: + if nuclides == 'all': + nuclides = self.get_nuclides() + if nuclides == 'sum': + nuclides = ['sum'] + else: + cv.check_iterable_type('nuclides', nuclides, string_types) + else: + nuclides = ['sum'] + + cv.check_value('xs_type', xs_type, ['macro', 'micro']) + + # Build header for string with type and domain info + string = 'Multi-Delayed-Group XS\n' + string += '{0: <16}=\t{1}\n'.format('\tReaction Type', self.rxn_type) + string += '{0: <16}=\t{1}\n'.format('\tDomain Type', self.domain_type) + string += '{0: <16}=\t{1}\n'.format('\tDomain ID', self.domain.id) + + # Generate the header for an individual XS + xs_header = '\tCross Sections [{0}]:'.format(self.get_units(xs_type)) + + # If cross section data has not been computed, only print string header + if self.tallies is None: + print(string) + return + + string += '{0: <16}\n'.format('\tEnergy Groups:') + template = '{0: <12}Group {1} [{2: <10} - {3: <10}MeV]\n' + + # Loop over energy groups ranges + for group in range(1, self.num_groups + 1): + bounds = self.energy_groups.get_group_bounds(group) + string += template.format('', group, bounds[0], bounds[1]) + + # Loop over all subdomains + for subdomain in subdomains: + + if self.domain_type == 'distribcell': + string += '{: <16}=\t{}\n'.format('\tSubdomain', subdomain) + + # Loop over all Nuclides + for nuclide in nuclides: + + # Build header for nuclide type + if xs_type != 'sum': + string += '{: <16}=\t{}\n'.format('\tNuclide', nuclide) + + # Build header for cross section type + string += '{: <16}\n'.format(xs_header) + + if self.delayed_groups is not None: + + for delayed_group in self.delayed_groups: + + template = '{0: <12}Delayed Group {1}:\t' + string += template.format('', delayed_group) + string += '\n' + + template = '{0: <12}Group {1} -> Group {2}:\t\t' + + # Loop over incoming/outgoing energy groups ranges + for in_group in range(1, self.num_groups + 1): + for out_group in range(1, self.num_groups + 1): + string += template.format('', in_group, out_group) + average = self.get_xs([in_group], [out_group], + [subdomain], [nuclide], + xs_type=xs_type, + value='mean', + delayed_groups=[delayed_group]) + rel_err = self.get_xs([in_group], [out_group], + [subdomain], [nuclide], + xs_type=xs_type, + value='rel_err', + delayed_groups=[delayed_group]) + average = average.flatten()[0] + rel_err = rel_err.flatten()[0] * 100. + string += '{:.2e} +/- {:.2e}%'.format(average, + rel_err) + string += '\n' + string += '\n' + string += '\n' + else: + + template = '{0: <12}Group {1} -> Group {2}:\t\t' + + # Loop over incoming/outgoing energy groups ranges + for in_group in range(1, self.num_groups + 1): + for out_group in range(1, self.num_groups + 1): + string += template.format('', in_group, out_group) + average = self.get_xs([in_group], [out_group], + [subdomain], [nuclide], + xs_type=xs_type, value='mean') + rel_err = self.get_xs([in_group], [out_group], + [subdomain], [nuclide], + xs_type=xs_type, value='rel_err') + average = average.flatten()[0] + rel_err = rel_err.flatten()[0] * 100. + string += '{:.2e} +/- {:.2e}%'.format(average, + rel_err) + string += '\n' + string += '\n' + string += '\n' + string += '\n' + + print(string) + + +class DelayedNuFissionMatrixXS(MatrixMDGXS): + r"""A fission delayed neutron production matrix multi-group cross section. + + This class can be used for both OpenMC input generation and tally data + post-processing to compute spatially-homogenized and energy-integrated + multi-group fission neutron production cross sections for multi-group + neutronics calculations. At a minimum, one needs to set the + :attr:`DelayedNuFissionMatrixXS.energy_groups` and + :attr:`DelayedNuFissionMatrixXS.domain` properties. Tallies for the flux and + appropriate reaction rates over the specified domain are generated + automatically via the :attr:`DelayedNuFissionMatrixXS.tallies` property, + which can then be appended to a :class:`openmc.Tallies` instance. + + For post-processing, the :meth:`MGXS.load_from_statepoint` will pull in the + necessary data to compute multi-group cross sections from a + :class:`openmc.StatePoint` instance. The derived multi-group cross section + can then be obtained from the :attr:`DelayedNuFissionMatrixXS.xs_tally` + property. + + For a spatial domain :math:`V`, energy group :math:`[E_g,E_{g-1}]`, and + delayed group :math:`d`, the fission delayed neutron production cross + section is calculated as: + + .. math:: + + \langle \nu\sigma_{f,g'\rightarrow g} \phi \rangle &= \int_{r \in V} dr + \int_{4\pi} d\Omega' \int_{E_{g'}}^{E_{g'-1}} dE' \int_{E_g}^{E_{g-1}} dE + \; \chi(E) \nu\sigma_f^d (r, E') \psi(r, E', \Omega')\\ + \langle \phi \rangle &= \int_{r \in V} dr \int_{4\pi} d\Omega + \int_{E_g}^{E_{g-1}} dE \; \psi (r, E, \Omega) \\ + \nu\sigma_{f,g'\rightarrow g} &= \frac{\langle \nu\sigma_{f,g'\rightarrow + g}^d \phi \rangle}{\langle \phi \rangle} + + + Parameters + ---------- + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + The domain for spatial homogenization + domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} + The domain type for spatial homogenization + groups : openmc.mgxs.EnergyGroups + The energy group structure for energy condensation + by_nuclide : bool + If true, computes cross sections for each nuclide in domain + name : str, optional + Name of the multi-group cross section. Used as a label to identify + tallies in OpenMC 'tallies.xml' file. + delayed_groups : list of int + Delayed groups to filter out the xs + + Attributes + ---------- + name : str, optional + Name of the multi-group cross section + rxn_type : str + Reaction type (e.g., 'total', 'nu-fission', etc.) + by_nuclide : bool + If true, computes cross sections for each nuclide in domain + domain : Material or Cell or Universe or Mesh + Domain for spatial homogenization + domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} + Domain type for spatial homogenization + energy_groups : openmc.mgxs.EnergyGroups + Energy group structure for energy condensation + delayed_groups : list of int + Delayed groups to filter out the xs + tally_trigger : openmc.Trigger + An (optional) tally precision trigger given to each tally used to + compute the cross section + scores : list of str + The scores in each tally used to compute the multi-group cross section + filters : list of openmc.Filter + The filters in each tally used to compute the multi-group cross section + tally_keys : list of str + The keys into the tallies dictionary for each tally used to compute + the multi-group cross section + estimator : {'tracklength', 'analog'} + The tally estimator used to compute the multi-group cross section + tallies : collections.OrderedDict + OpenMC tallies needed to compute the multi-group cross section. The keys + are strings listed in the :attr:`DelayedNuFissionXS.tally_keys` property + and values are instances of :class:`openmc.Tally`. + rxn_rate_tally : openmc.Tally + Derived tally for the reaction rate tally used in the numerator to + compute the multi-group cross section. This attribute is None + unless the multi-group cross section has been computed. + xs_tally : openmc.Tally + Derived tally for the multi-group cross section. This attribute + is None unless the multi-group cross section has been computed. + num_subdomains : int + The number of subdomains is unity for 'material', 'cell' and 'universe' + domain types. When the This is equal to the number of cell instances + for 'distribcell' domain types (it is equal to unity prior to loading + tally data from a statepoint file). + num_nuclides : int + The number of nuclides for which the multi-group cross section is + being tracked. This is unity if the by_nuclide attribute is False. + nuclides : Iterable of str or 'sum' + The optional user-specified nuclides for which to compute cross + sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + are not specified by the user, all nuclides in the spatial domain + are included. This attribute is 'sum' if by_nuclide is false. + sparse : bool + Whether or not the MGXS' tallies use SciPy's LIL sparse matrix format + for compressed data storage + loaded_sp : bool + Whether or not a statepoint file has been loaded with tally data + derived : bool + Whether or not the MGXS is merged from one or more other MGXS + hdf5_key : str + The key used to index multi-group cross sections in an HDF5 data store + + """ + + def __init__(self, domain=None, domain_type=None, energy_groups=None, + delayed_groups=None, by_nuclide=False, name=''): + super(DelayedNuFissionMatrixXS, self).__init__(domain, domain_type, + energy_groups, + delayed_groups, + by_nuclide, name) + self._rxn_type = 'delayed-nu-fission' + self._hdf5_key = 'delayed-nu-fission matrix' + self._estimator = 'analog' + self._valid_estimators = ['analog'] diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 06527d5c2..db31e2f32 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -36,7 +36,8 @@ MGXS_TYPES = ['total', 'chi', 'chi-prompt', 'inverse-velocity', - 'prompt-nu-fission'] + 'prompt-nu-fission', + 'prompt-nu-fission matrix'] # Supported domain types DOMAIN_TYPES = ['cell', @@ -451,7 +452,7 @@ class MGXS(object): Parameters ---------- - mgxs_type : {'total', 'transport', 'nu-transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'multiplicity matrix', 'nu-fission matrix', 'chi', 'chi-prompt', 'inverse-velocity', 'prompt-nu-fission'} + mgxs_type : {'total', 'transport', 'nu-transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'multiplicity matrix', 'nu-fission matrix', 'chi', 'chi-prompt', 'inverse-velocity', 'prompt-nu-fission', 'prompt-nu-fission matrix'} The type of multi-group cross section object to return domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh The domain for spatial homogenization @@ -512,6 +513,8 @@ class MGXS(object): mgxs = InverseVelocity(domain, domain_type, energy_groups) elif mgxs_type == 'prompt-nu-fission': mgxs = PromptNuFissionXS(domain, domain_type, energy_groups) + elif mgxs_type == 'prompt-nu-fission matrix': + mgxs = PromptNuFissionMatrixXS(domain, domain_type, energy_groups) mgxs.by_nuclide = by_nuclide mgxs.name = name @@ -1755,7 +1758,7 @@ class MatrixMGXS(MGXS): Returns ------- - ndarray + numpy.ndarray A NumPy array of the multi-group cross section indexed in the order each group and subdomain is listed in the parameters. @@ -3654,7 +3657,7 @@ class ScatterMatrixXS(MatrixMGXS): Returns ------- - ndarray + numpy.ndarray A NumPy array of the multi-group cross section indexed in the order each group and subdomain is listed in the parameters. @@ -5232,3 +5235,120 @@ class PromptNuFissionXS(MGXS): super(PromptNuFissionXS, self).__init__(domain, domain_type, groups, by_nuclide, name) self._rxn_type = 'prompt-nu-fission' + + +class PromptNuFissionMatrixXS(MatrixMGXS): + r"""A prompt fission neutron production matrix multi-group cross section. + + This class can be used for both OpenMC input generation and tally data + post-processing to compute spatially-homogenized and energy-integrated + multi-group cross sections for multi-group neutronics calculations. At a + minimum, one needs to set the :attr:`PromptNuFissionMatrixXS.energy_groups` + and :attr:`PromptNuFissionMatrixXS.domain` properties. Tallies for the flux + and appropriate reaction rates over the specified domain are generated + automatically via the :attr:`PromptNuFissionMatrixXS.tallies` property, + which can then be appended to a :class:`openmc.Tallies` instance. + + For post-processing, the :meth:`MGXS.load_from_statepoint` will pull in the + necessary data to compute multi-group cross sections from a + :class:`openmc.StatePoint` instance. The derived multi-group cross section + can then be obtained from the :attr:`PromptNuFissionMatrixXS.xs_tally` + property. + + For a spatial domain :math:`V` and energy group :math:`[E_g,E_{g-1}]`, the + fission spectrum is calculated as: + + .. math:: + + \langle \nu\sigma_{f,g'\rightarrow g} \phi \rangle &= \int_{r \in V} dr + \int_{4\pi} d\Omega' \int_{E_{g'}}^{E_{g'-1}} dE' \int_{E_g}^{E_{g-1}} dE + \; \chi(E) \nu\sigma_f^p (r, E') \psi(r, E', \Omega')\\ + \langle \phi \rangle &= \int_{r \in V} dr \int_{4\pi} d\Omega + \int_{E_g}^{E_{g-1}} dE \; \psi (r, E, \Omega) \\ + \nu\sigma_{f,g'\rightarrow g} &= \frac{\langle \nu\sigma_{f,g'\rightarrow + g}^p \phi \rangle}{\langle \phi \rangle} + + Parameters + ---------- + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + The domain for spatial homogenization + domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} + The domain type for spatial homogenization + groups : openmc.mgxs.EnergyGroups + The energy group structure for energy condensation + by_nuclide : bool + If true, computes cross sections for each nuclide in domain + name : str, optional + Name of the multi-group cross section. Used as a label to identify + tallies in OpenMC 'tallies.xml' file. + + Attributes + ---------- + name : str, optional + Name of the multi-group cross section + rxn_type : str + Reaction type (e.g., 'total', 'nu-fission', etc.) + by_nuclide : bool + If true, computes cross sections for each nuclide in domain + domain : Material or Cell or Universe or Mesh + Domain for spatial homogenization + domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} + Domain type for spatial homogenization + energy_groups : openmc.mgxs.EnergyGroups + Energy group structure for energy condensation + tally_trigger : openmc.Trigger + An (optional) tally precision trigger given to each tally used to + compute the cross section + scores : list of str + The scores in each tally used to compute the multi-group cross section + filters : list of openmc.Filter + The filters in each tally used to compute the multi-group cross section + tally_keys : list of str + The keys into the tallies dictionary for each tally used to compute + the multi-group cross section + estimator : {'tracklength', 'collision', 'analog'} + The tally estimator used to compute the multi-group cross section + tallies : collections.OrderedDict + OpenMC tallies needed to compute the multi-group cross section. The keys + are strings listed in the :attr:`PromptNuFissionXS.tally_keys` property + and values are instances of :class:`openmc.Tally`. + rxn_rate_tally : openmc.Tally + Derived tally for the reaction rate tally used in the numerator to + compute the multi-group cross section. This attribute is None + unless the multi-group cross section has been computed. + xs_tally : openmc.Tally + Derived tally for the multi-group cross section. This attribute + is None unless the multi-group cross section has been computed. + num_subdomains : int + The number of subdomains is unity for 'material', 'cell' and 'universe' + domain types. This is equal to the number of cell instances + for 'distribcell' domain types (it is equal to unity prior to loading + tally data from a statepoint file). + num_nuclides : int + The number of nuclides for which the multi-group cross section is + being tracked. This is unity if the by_nuclide attribute is False. + nuclides : Iterable of str or 'sum' + The optional user-specified nuclides for which to compute cross + sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + are not specified by the user, all nuclides in the spatial domain + are included. This attribute is 'sum' if by_nuclide is false. + sparse : bool + Whether or not the MGXS' tallies use SciPy's LIL sparse matrix format + for compressed data storage + loaded_sp : bool + Whether or not a statepoint file has been loaded with tally data + derived : bool + Whether or not the MGXS is merged from one or more other MGXS + hdf5_key : str + The key used to index multi-group cross sections in an HDF5 data store + + """ + + def __init__(self, domain=None, domain_type=None, + groups=None, by_nuclide=False, name=''): + super(PromptNuFissionMatrixXS, self).__init__(domain, domain_type, + groups, by_nuclide, name) + self._rxn_type = 'prompt-nu-fission' + self._hdf5_key = 'prompt-nu-fission matrix' + self._estimator = 'analog' + self._valid_estimators = ['analog'] diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 6e99030e0..ef1dd236c 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1,3 +1,4 @@ +from collections import Iterable from numbers import Real, Integral import os @@ -15,7 +16,7 @@ from openmc.checkvalue import check_type, check_value, check_greater_than, \ _REPRESENTATIONS = ['isotropic', 'angle'] _SCATTER_TYPES = ['tabular', 'legendre', 'histogram'] _XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]", - "[G'][DG]"] + "[G'][DG]", "[G][G'][DG]"] class XSdata(object): @@ -136,11 +137,11 @@ class XSdata(object): [G][G'][Order]: scatter_matrix [G]: total, absorption, fission, kappa_fission, nu_fission, - prompt_nu_fission, inverse_velocity + prompt_nu_fission, delayed_nu_fission, inverse_velocity [G']: chi, chi_prompt, chi_delayed - [G][G']: multiplicity_matrix, nu_fission + [G][G']: multiplicity_matrix, nu_fission, prompt_nu_fission [DG]: beta, decay_rate @@ -148,6 +149,8 @@ class XSdata(object): [G'][DG]: chi_delayed + [G][G'][DG]: delayed_nu_fission + """ def __init__(self, name, energy_groups, temperatures=[294.], @@ -297,6 +300,10 @@ class XSdata(object): self.num_delayed_groups) self._xs_shapes["[G'][DG]"] = (self.energy_groups.num_groups, self.num_delayed_groups) + self._xs_shapes["[G][G'][DG]"] = (self.energy_groups.num_groups, + self.energy_groups.num_groups, + self.num_delayed_groups) + self._xs_shapes["[G][G'][Order]"] \ = (self.energy_groups.num_groups, self.energy_groups.num_groups, self.num_orders) @@ -822,7 +829,7 @@ class XSdata(object): """ # Get the accepted shapes for this xs - shapes = [self.xs_shapes["[G]"]] + shapes = [self.xs_shapes["[G]"], self.xs_shapes["[G][G']"]] # Convert to a numpy array so we can easily get the shape for checking prompt_nu_fission = np.asarray(prompt_nu_fission) @@ -856,7 +863,7 @@ class XSdata(object): """ # Get the accepted shapes for this xs - shapes = [self.xs_shapes["[G][DG]"]] + shapes = [self.xs_shapes["[G][DG]"], self.xs_shapes["[G][G'][DG]"]] # Convert to a numpy array so we can easily get the shape for checking delayed_nu_fission = np.asarray(delayed_nu_fission) @@ -1092,12 +1099,15 @@ class XSdata(object): def set_prompt_nu_fission_mgxs(self, prompt_nu_fission, temperature=294., nuclide='total', xs_type='macro', subdomain=None): - """This method allows for an openmc.mgxs.PromptNuFissionXS to be used to - set the prompt-nu-fission cross section for this XSdata object. + """Sets the prompt-nu-fission cross section. + + This method allows for an openmc.mgxs.PromptNuFissionXS or + openmc.mgxs.PromptNuFissionMatrixXS to be used to set the + prompt-nu-fission cross section for this XSdata object. Parameters ---------- - prompt_nu_fission: openmc.mgxs.PromptNuFissionXS + prompt_nu_fission: openmc.mgxs.PromptNuFissionXS or openmc.mgxs.PromptNuFissionMatrixXS MGXS Object containing the prompt-nu-fission cross section for the domain of interest. temperature : float @@ -1121,7 +1131,8 @@ class XSdata(object): """ check_type('prompt_nu_fission', prompt_nu_fission, - (openmc.mgxs.PromptNuFissionXS,)) + (openmc.mgxs.PromptNuFissionXS, + openmc.mgxs.PromptNuFissionMatrixXS)) check_value('energy_groups', prompt_nu_fission.energy_groups, [self.energy_groups]) check_value('domain_type', prompt_nu_fission.domain_type, @@ -1145,12 +1156,13 @@ class XSdata(object): def set_delayed_nu_fission_mgxs(self, delayed_nu_fission, temperature=294., nuclide='total', xs_type='macro', subdomain=None): - """This method allows for an openmc.mgxs.DelayedNuFissionXS to be used - to set the delayed-nu-fission cross section for this XSdata object. + """This method allows for an openmc.mgxs.DelayedNuFissionXS or + openmc.mgxs.DelayedNuFissionMatrixXS to be used to set the + delayed-nu-fission cross section for this XSdata object. Parameters ---------- - delayed_nu_fission: openmc.mgxs.DelayedNuFissionXS + delayed_nu_fission: openmc.mgxs.DelayedNuFissionXS or openmc.mgxs.DelayedNuFissionMatrixXS MGXS Object containing the delayed-nu-fission cross section for the domain of interest. temperature : float @@ -1174,11 +1186,11 @@ class XSdata(object): """ check_type('delayed_nu_fission', delayed_nu_fission, - (openmc.mgxs.DelayedNuFissionXS,)) + (openmc.mgxs.DelayedNuFissionXS, + openmc.mgxs.DelayedNuFissionMatrixXS)) check_value('energy_groups', delayed_nu_fission.energy_groups, [self.energy_groups]) - check_value('num_delayed_groups', - delayed_nu_fission.num_delayed_groups, + check_value('num_delayed_groups', delayed_nu_fission.num_delayed_groups, [self.num_delayed_groups]) check_value('domain_type', delayed_nu_fission.domain_type, openmc.mgxs.DOMAIN_TYPES) @@ -2120,46 +2132,3 @@ class MGXSLibrary(object): xsdata.to_hdf5(file) file.close() - - @classmethod - def from_hdf5(cls, filename=None): - """Generate an MGXS Library from an HDF5 group or file - - Parameters - ---------- - filename : str, optional - Name of HDF5 file containing MGXS data. Default is None. - If not provided, the value of the OPENMC_MG_CROSS_SECTIONS - environmental variable will be used - - Returns - ------- - openmc.MGXSLibrary - Multi-group cross section data object. - - """ - - # If filename is None, get the cross sections from the - # OPENMC_CROSS_SECTIONS environment variable - if filename is None: - filename = os.environ.get('OPENMC_MG_CROSS_SECTIONS') - - # Check to make sure there was an environmental variable. - if filename is None: - raise ValueError("Either path or OPENMC_MG_CROSS_SECTIONS " - "environmental variable must be set") - - check_type('filename', filename, str) - file = h5py.File(filename, 'r') - - group_structure = file.attrs['group structure'] - num_delayed_groups = file.attrs['delayed_groups'] - energy_groups = openmc.mgxs.EnergyGroups(group_structure) - data = cls(energy_groups, num_delayed_groups) - - for group_name, group in file.items(): - data.add_xsdata(openmc.XSdata.from_hdf5(group, group_name, - energy_groups, - num_delayed_groups)) - - return data diff --git a/openmc/particle_restart.py b/openmc/particle_restart.py index a0410298b..cb1326586 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/plotter.py b/openmc/plotter.py index c19d544b2..aa48ef944 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -1,4 +1,5 @@ from numbers import Integral, Real +from six import string_types import numpy as np @@ -19,42 +20,18 @@ UNITY_MT = -1 XI_MT = -2 # MTs to combine to generate associated plot_types -PLOT_TYPES_MT = {'total': (2, 3,), - 'scatter': (2, 4, 11, 16, 17, 22, 23, 24, 25, 28, 29, 30, - 32, 33, 34, 35, 36, 37, 41, 42, 44, 45, 152, - 153, 154, 156, 157, 158, 159, 160, 161, 162, - 163, 164, 165, 166, 167, 168, 169, 170, 171, - 172, 173, 174, 175, 176, 177, 178, 179, 180, - 181, 183, 184, 190, 194, 196, 198, 199, 200, - 875, 891), - 'elastic': (2,), - 'inelastic': (4, 11, 16, 17, 22, 23, 24, 25, 28, 29, 30, - 32, 33, 34, 35, 36, 37, 41, 42, 44, 45, 152, - 153, 154, 156, 157, 158, 159, 160, 161, 162, - 163, 164, 165, 166, 167, 168, 169, 170, 171, - 172, 173, 174, 175, 176, 177, 178, 179, 180, - 181, 183, 184, 190, 194, 196, 198, 199, 200, - 875, 891), - 'fission': (18,), - 'absorption': (27,), 'capture': (101,), - 'nu-fission': (18,), - 'nu-scatter': (2, 4, 11, 16, 17, 22, 23, 24, 25, 28, 29, 30, - 32, 33, 34, 35, 36, 37, 41, 42, 44, 45, 152, - 153, 154, 156, 157, 158, 159, 160, 161, 162, - 163, 164, 165, 166, 167, 168, 169, 170, 171, - 172, 173, 174, 175, 176, 177, 178, 179, 180, - 181, 183, 184, 190, 194, 196, 198, 199, 200, - 875, 891), - 'unity': (UNITY_MT,), - 'slowing-down power': (2, 4, 11, 16, 17, 22, 23, 24, 25, 28, - 29, 30, 32, 33, 34, 35, 36, 37, 41, 42, - 44, 45, 152, 153, 154, 156, 157, 158, - 159, 160, 161, 162, 163, 164, 165, 166, - 167, 168, 169, 170, 171, 172, 173, 174, - 175, 176, 177, 178, 179, 180, 181, 183, - 184, 190, 194, 196, 198, 199, 200, 875, - 891, XI_MT), - 'damage': (444,)} +_INELASTIC = [mt for mt in openmc.data.SUM_RULES[3] if mt != 27] +PLOT_TYPES_MT = {'total': openmc.data.SUM_RULES[1], + 'scatter': [2] + _INELASTIC, + 'elastic': [2], + 'inelastic': _INELASTIC, + 'fission': [18], + 'absorption': [27], 'capture': [101], + 'nu-fission': [18], + 'nu-scatter': [2] + _INELASTIC, + 'unity': [UNITY_MT], + 'slowing-down power': [2] + _INELASTIC + [XI_MT], + 'damage': [444]} # Operations to use when combining MTs the first np.add is used in reference # to zero PLOT_TYPES_OP = {'total': (np.add,), @@ -69,27 +46,13 @@ PLOT_TYPES_OP = {'total': (np.add,), (np.add,) * (len(PLOT_TYPES_MT['slowing-down power']) - 2) + (np.multiply,), 'damage': ()} -# Whether or not to multiply the reaction by the yield as well -PLOT_TYPES_YIELD = {'total': (False, False), - 'scatter': (False,) * len(PLOT_TYPES_MT['scatter']), - 'elastic': (False,), - 'inelastic': (False,) * len(PLOT_TYPES_MT['inelastic']), - 'fission': (False,), 'absorption': (False,), - 'capture': (False,), 'nu-fission': (True,), - 'nu-scatter': (True,) * len(PLOT_TYPES_MT['nu-scatter']), - 'unity': (False,), - 'slowing-down power': - (True,) * len(PLOT_TYPES_MT['slowing-down power']), - 'damage': (False,)} - # Types of plots to plot linearly in y PLOT_TYPES_LINEAR = {'nu-fission / fission', 'nu-scatter / scatter', 'nu-fission / absorption', 'fission / absorption'} def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, - energy_range=(1.E-5, 20.E6), sab_name=None, cross_sections=None, - enrichment=None, **kwargs): + sab_name=None, cross_sections=None, enrichment=None, **kwargs): """Creates a figure of continuous-energy cross sections for this item Parameters @@ -110,8 +73,6 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, axis : matplotlib.axes, optional A previously generated axis to use for plotting. If not specified, a new axis and figure will be generated. - energy_range : tuple of floats - Energy range (in eV) to plot the cross section within sab_name : str, optional Name of S(a,b) library to apply to MT=2 data when applicable; only used for items which are instances of openmc.Element or openmc.Nuclide @@ -196,6 +157,7 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, to_plot = data[i](E) else: to_plot = data[i, :] + to_plot = np.nan_to_num(to_plot) if np.sum(to_plot) > 0.: plot_func(E, to_plot, label=types[i]) @@ -216,10 +178,10 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, ylabel = 'Macroscopic Cross Section [1/cm]' ax.set_ylabel(ylabel) ax.legend(loc='best') - ax.set_xlim(energy_range) + # Set to the most likely expected range + ax.set_xlim((1.E-5, 20.E6)) if this.name is not None: - title = 'Cross Section for ' + this.name - ax.set_title(title) + ax.set_title('Cross Section for ' + this.name) return fig @@ -384,7 +346,7 @@ def calculate_xs(this, types, temperature=294., sab_name=None, Returns ------- - energy_grid : numpy.array + energy_grid : numpy.ndarray Energies at which cross sections are calculated, in units of eV data : numpy.ndarray Cross sections calculated at the energy grid described by energy_grid @@ -394,7 +356,7 @@ def calculate_xs(this, types, temperature=294., sab_name=None, # Check types cv.check_type('temperature', temperature, Real) if sab_name: - cv.check_type('sab_name', sab_name, str) + cv.check_type('sab_name', sab_name, string_types) if enrichment: cv.check_type('enrichment', enrichment, Real) @@ -440,7 +402,7 @@ def _calculate_xs_element(this, types, temperature=294., sab_name=None, Returns ------- - energy_grid : numpy.array + energy_grid : numpy.ndarray Energies at which cross sections are calculated, in units of eV data : numpy.ndarray Macroscopic cross sections calculated at the energy grid described @@ -452,7 +414,7 @@ def _calculate_xs_element(this, types, temperature=294., sab_name=None, library = openmc.data.DataLibrary.from_xml(cross_sections) # Expand elements in to nuclides with atomic densities - nuclides = this.expand(100., 'ao', enrichment=enrichment, + nuclides = this.expand(1., 'ao', enrichment=enrichment, cross_sections=cross_sections) # For ease of processing split out nuc and nuc_density @@ -519,7 +481,7 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, Returns ------- - energy_grid : numpy.array + energy_grid : numpy.ndarray Energies at which cross sections are calculated, in units of eV data : numpy.ndarray Cross sections calculated at the energy grid described by @@ -534,15 +496,18 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, for line in types: if line in PLOT_TYPES: mts.append(PLOT_TYPES_MT[line]) - yields.append(PLOT_TYPES_YIELD[line]) + if line.startswith('nu'): + yields.append(True) + else: + yields.append(False) ops.append(PLOT_TYPES_OP[line]) else: # Not a built-in type, we have to parse it ourselves cv.check_type('MT in types', line, Integral) cv.check_greater_than('MT in types', line, 0) mts.append((line,)) - yields.append((False,)) ops.append(()) + yields.append(False) # Load the library library = openmc.data.DataLibrary.from_xml(cross_sections) @@ -565,7 +530,7 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, for t in range(len(data_Ts)): # Take off the "K" and convert to a float data_Ts[t] = float(data_Ts[t][:-1]) - min_delta = np.finfo(np.float64).max + min_delta = float('inf') closest_t = -1 for t in data_Ts: if abs(data_Ts[t] - T) < min_delta: @@ -620,7 +585,7 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, # Get the reaction xs data from the nuclide funcs = [] op = ops[i] - for mt, yield_check in zip(mt_set, yields[i]): + for mt in mt_set: if mt == 2: if sab_name: # Then we need to do a piece-wise function of @@ -632,29 +597,36 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, funcs.append(pw_funcs) else: funcs.append(nuc[mt].xs[nucT]) + elif mt == 5 and mt in nuc: + # Only consider the (n,misc) products if neutrons are + # included in the outgoing channel since (n,misc) is only + # explicitly needed for scatter cross sections. + for prod in nuc[mt].products: + if prod.particle == 'neutron' and \ + prod.emission_mode in ('total', 'prompt'): + if yields[i]: + func = openmc.data.Combination( + [nuc[mt].xs[nucT], prod.yield_], + [np.multiply]) + else: + func = nuc[mt].xs[nucT] + + funcs.append(func) + break + else: + funcs.append(lambda x: 0.) + elif mt in nuc: - if yield_check: - found_it = False + if yields[i]: for prod in nuc[mt].products: if prod.particle == 'neutron' and \ - prod.emission_mode == 'total': + prod.emission_mode in ('total', 'prompt'): func = openmc.data.Combination( [nuc[mt].xs[nucT], prod.yield_], [np.multiply]) funcs.append(func) - found_it = True break - if not found_it: - for prod in nuc[mt].products: - if prod.particle == 'neutron' and \ - prod.emission_mode == 'prompt': - func = openmc.data.Combination( - [nuc[mt].xs[nucT], - prod.yield_], [np.multiply]) - funcs.append(func) - found_it = True - break - if not found_it: + else: # Assume the yield is 1 funcs.append(nuc[mt].xs[nucT]) else: @@ -695,7 +667,7 @@ def _calculate_xs_material(this, types, temperature=294., cross_sections=None): Returns ------- - energy_grid : numpy.array + energy_grid : numpy.ndarray Energies at which cross sections are calculated, in units of eV data : numpy.ndarray Macroscopic cross sections calculated at the energy grid described @@ -712,7 +684,7 @@ def _calculate_xs_material(this, types, temperature=294., cross_sections=None): library = openmc.data.DataLibrary.from_xml(cross_sections) # Expand elements in to nuclides with atomic densities - nuclides = this.get_nuclide_atom_densities(cross_sections) + nuclides = this.get_nuclide_atom_densities() # For ease of processing split out nuc and nuc_density nuc_densities = [nuclide[1][1] for nuclide in nuclides.items()] @@ -740,8 +712,8 @@ def _calculate_xs_material(this, types, temperature=294., cross_sections=None): # Condense the data for every nuclide # First create a union energy grid energy_grid = E[0] - for n in range(1, len(E)): - energy_grid = np.union1d(energy_grid, E[n]) + for grid in E[1:]: + energy_grid = np.union1d(energy_grid, grid) # Now we can combine all the nuclidic data data = np.zeros((len(types), len(energy_grid))) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index aaa5c0b84..144146ba3 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 @@ -96,6 +97,9 @@ class StatePoint(object): Dictionary whose keys are tally IDs and whose values are Tally objects tallies_present : bool Indicate whether user-defined tallies are present + tally_derivatives : dict + Dictionary whose keys are tally derivative IDs and whose values are + TallyDerivative objects version: tuple of Integral Version of OpenMC summary : None or openmc.Summary @@ -104,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 @@ -134,6 +132,7 @@ class StatePoint(object): self._summary = None self._global_tallies = None self._sparse = False + self._derivs_read = False # Automatically link in a summary file if one exists if autolink: @@ -209,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 @@ -396,6 +395,12 @@ class StatePoint(object): base, tally_key)].value.decode() tally.num_realizations = n_realizations + # Read derivative information. + if 'derivative' in self._f['{0}{1}'.format(base, tally_key)]: + deriv_id = self._f['{0}{1}/derivative'.format( + base, tally_key)].value + tally.derivative = self.tally_derivatives[deriv_id] + # Read the number of Filters n_filters = \ self._f['{0}{1}/n_filters'.format(base, tally_key)].value @@ -457,6 +462,43 @@ class StatePoint(object): def tallies_present(self): return self._f['tallies/tallies_present'].value + @property + def tally_derivatives(self): + if not self._derivs_read: + # Initialize dictionaries for the Meshes + # Keys - Derivative IDs + # Values - TallyDerivative objects + self._derivs = {} + + # Populate the dictionary if any derivatives are present. + if 'derivatives' in self._f['tallies']: + # Read the derivative ids. + base = 'tallies/derivatives' + deriv_ids = [int(k.split(' ')[1]) for k in self._f[base]] + + # Create each derivative object and add it to the dictionary. + for d_id in deriv_ids: + base = 'tallies/derivatives/derivative {:d}'.format(d_id) + deriv = openmc.TallyDerivative(derivative_id=d_id) + deriv.variable = \ + self._f[base + '/independent variable'].value.decode() + if deriv.variable == 'density': + deriv.material = self._f[base + '/material'].value + elif deriv.variable == 'nuclide_density': + deriv.material = self._f[base + '/material'].value + deriv.nuclide = \ + self._f[base + '/nuclide'].value.decode() + elif deriv.variable == 'temperature': + deriv.material = self._f[base + '/material'].value + else: + raise RuntimeError('Unrecognized tally differential ' + 'variable') + self._derivs[d_id] = deriv + + self._derivs_read = True + + return self._derivs + @property def version(self): return (self._f['version_major'].value, diff --git a/openmc/summary.py b/openmc/summary.py index fd1689eb2..37509ec37 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 ac959080e..979929a26 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 @@ -103,6 +104,8 @@ class Tally(object): sparse : bool Whether or not the tally uses SciPy's LIL sparse matrix format for compressed data storage + derivative : openmc.TallyDerivative + A material perturbation derivative to apply to all scores in the tally. """ @@ -115,6 +118,7 @@ class Tally(object): self._scores = cv.CheckedList(_SCORE_CLASSES, 'tally scores') self._estimator = None self._triggers = cv.CheckedList(openmc.Trigger, 'tally triggers') + self._derivative = None self._num_realizations = 0 self._with_summary = False @@ -150,6 +154,10 @@ class Tally(object): if nuclide not in other.nuclides: return False + # Check derivatives + if self.derivative != other.derivative: + return False + # Check all scores if len(self.scores) != len(other.scores): return False @@ -171,27 +179,31 @@ class Tally(object): def __repr__(self): string = 'Tally\n' - string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self.id) - string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self.name) + string += '{: <16}=\t{}\n'.format('\tID', self.id) + string += '{: <16}=\t{}\n'.format('\tName', self.name) - string += '{0: <16}{1}\n'.format('\tFilters', '=\t') + if self.derivative is not None: + string += '{: <16}=\t{}\n'.format('\tDerivative ID', + str(self.derivative.id)) + + string += '{: <16}=\n'.format('\tFilters') for self_filter in self.filters: - string += '{0: <16}\t\t{1}\t{2}\n'.format('', + string += '{: <16}\t\t{}\t{}\n'.format('', type(self_filter).__name__, self_filter.bins) - string += '{0: <16}{1}'.format('\tNuclides', '=\t') + string += '{: <16}=\t'.format('\tNuclides') for nuclide in self.nuclides: if isinstance(nuclide, openmc.Nuclide): - string += '{0} '.format(nuclide.name) + string += nuclide.name + ' ' else: - string += '{0} '.format(nuclide) + string += nuclide + ' ' string += '\n' - string += '{0: <16}{1}{2}\n'.format('\tScores', '=\t', self.scores) - string += '{0: <16}{1}{2}\n'.format('\tEstimator', '=\t', self.estimator) + string += '{: <16}=\t{}\n'.format('\tScores', self.scores) + string += '{: <16}=\t{}\n'.format('\tEstimator', self.estimator) return string @@ -269,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) @@ -380,6 +386,10 @@ class Tally(object): def derived(self): return self._derived + @property + def derivative(self): + return self._derivative + @property def sparse(self): return self._sparse @@ -434,6 +444,12 @@ class Tally(object): else: self._name = '' + @derivative.setter + def derivative(self, deriv): + if deriv is not None: + cv.check_type('tally derivative', deriv, openmc.TallyDerivative) + self._derivative = deriv + @filters.setter def filters(self, filters): cv.check_type('tally filters', filters, MutableSequence) @@ -1060,6 +1076,11 @@ class Tally(object): for trigger in self.triggers: trigger.get_trigger_xml(element) + # Optional derivatives + if self.derivative is not None: + subelement = ET.SubElement(element, "derivative") + subelement.text = str(self.derivative.id) + return element def contains_filter(self, filter_type): @@ -1495,7 +1516,8 @@ class Tally(object): return data def get_pandas_dataframe(self, filters=True, nuclides=True, scores=True, - distribcell_paths=True, float_format='{:.2e}'): + derivative=True, distribcell_paths=True, + float_format='{:.2e}'): """Build a Pandas DataFrame for the Tally data. This method constructs a Pandas DataFrame object for the Tally data @@ -1513,6 +1535,8 @@ class Tally(object): Include columns with nuclide bin information (default is True). scores : bool Include columns with score bin information (default is True). + derivative : bool + Include columns with differential tally info (default is True). distribcell_paths : bool, optional Construct columns for distribcell tally filters (default is True). The geometric information in the Summary object is embedded into a @@ -1593,6 +1617,14 @@ class Tally(object): tile_factor = data_size / len(self.scores) df[column_name] = np.tile(scores, int(tile_factor)) + # Include columns for derivatives if user requested it + if derivative and (self.derivative is not None): + df['d_variable'] = self.derivative.variable + if self.derivative.material is not None: + df['d_material'] = self.derivative.material + if self.derivative.nuclide is not None: + df['d_nuclide'] = self.derivative.nuclide + # Append columns with mean, std. dev. for each tally bin df['mean'] = self.mean.ravel() df['std. dev.'] = self.std_dev.ravel() @@ -1727,8 +1759,6 @@ class Tally(object): # HDF5 binary file if format == 'hdf5': - import h5py - filename = directory + '/' + filename + '.h5' if append: @@ -3562,6 +3592,18 @@ class Tallies(cv.CheckedList): self._tallies_file.append(xml_element) already_written.add(f.mesh) + def _create_derivative_subelements(self): + # Get a list of all derivatives referenced in a tally. + derivs = [] + for tally in self: + deriv = tally.derivative + if deriv is not None and deriv not in derivs: + derivs.append(deriv) + + # Add the derivatives to the XML tree. + for d in derivs: + self._tallies_file.append(d.to_xml_element()) + def export_to_xml(self): """Create a tallies.xml file that can be used for a simulation. @@ -3572,6 +3614,7 @@ class Tallies(cv.CheckedList): self._create_mesh_subelements() self._create_tally_subelements() + self._create_derivative_subelements() # Clean the indentation in the file to be user-readable clean_xml_indentation(self._tallies_file) diff --git a/openmc/tally_derivative.py b/openmc/tally_derivative.py new file mode 100644 index 000000000..1c8a316cf --- /dev/null +++ b/openmc/tally_derivative.py @@ -0,0 +1,141 @@ +from __future__ import division + +import sys +from numbers import Integral +from xml.etree import ElementTree as ET + +from six import string_types + +import openmc.checkvalue as cv +from openmc.mixin import EqualityMixin + + +# "Static" variable for auto-generated TallyDerivative IDs +AUTO_TALLY_DERIV_ID = 10000 + +def reset_auto_tally_deriv_id(): + global AUTO_TALLY_ID + AUTO_TALLY_DERIV_ID = 10000 + + +class TallyDerivative(EqualityMixin): + """A material perturbation derivative to apply to a tally. + + Parameters + ---------- + derivative_id : Integral, optional + Unique identifier for the tally derivative. If none is specified, an + identifier will automatically be assigned + variable : str, optional + Accepted values are 'density', 'nuclide_density', and 'temperature' + material : Integral, optional + The perturubed material ID + nuclide : str, optional + The perturbed nuclide. Only needed for 'nuclide_density' derivatives. + Ex: 'Xe135' + + Attributes + ---------- + id : Integral + Unique identifier for the tally derivative + variable : str + Accepted values are 'density', 'nuclide_density', and 'temperature' + material : Integral + The perturubed material ID + nuclide : str + The perturbed nuclide. Only needed for 'nuclide_density' derivatives. + Ex: 'Xe135' + + """ + + def __init__(self, derivative_id=None, variable=None, material=None, + nuclide=None): + # Initialize Tally class attributes + self.id = derivative_id + self.variable = variable + self.material = material + self.nuclide = nuclide + + def __hash__(self): + return hash(repr(self)) + + def __repr__(self): + string = 'Tally Derivative\n' + string += '{: <16}=\t{}\n'.format('\tID', self.id) + string += '{: <16}=\t{}\n'.format('\tVariable', self.variable) + + if self.variable == 'density': + string += '{: <16}=\t{}\n'.format('\tMaterial', self.material) + elif self.variable == 'nuclide_density': + string += '{: <16}=\t{}\n'.format('\tMaterial', self.material) + string += '{: <16}=\t{}\n'.format('\tNuclide', self.nuclide) + elif self.variable == 'temperature': + string += '{: <16}=\t{}\n'.format('\tMaterial', self.material) + + return string + + @property + def id(self): + return self._id + + @property + def variable(self): + return self._variable + + @property + def material(self): + return self._material + + @property + def nuclide(self): + return self._nuclide + + @id.setter + def id(self, deriv_id): + if deriv_id is None: + global AUTO_TALLY_DERIV_ID + self._id = AUTO_TALLY_DERIV_ID + AUTO_TALLY_DERIV_ID += 1 + else: + cv.check_type('tally derivative ID', deriv_id, Integral) + cv.check_greater_than('tally derivative ID', deriv_id, 0, + equality=True) + self._id = deriv_id + + @variable.setter + def variable(self, var): + if var is not None: + cv.check_type('derivative variable', var, string_types) + cv.check_value('derivative variable', var, + ('density', 'nuclide_density', 'temperature')) + self._variable = var + + @material.setter + def material(self, mat): + if mat is not None: + cv.check_type('derivative material', mat, Integral) + self._material = mat + + @nuclide.setter + def nuclide(self, nuc): + if nuc is not None: + cv.check_type('derivative nuclide', nuc, string_types) + self._nuclide = nuc + + def to_xml_element(self): + """Return XML representation of the tally derivative + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing derivative data + + """ + + element = ET.Element("derivative") + element.set("id", str(self.id)) + element.set("variable", self.variable) + element.set("material", str(self.material)) + if self.variable == 'nuclide_density': + element.set("nuclide", self.nuclide) + return element diff --git a/openmc/volume.py b/openmc/volume.py index 2a2a92f61..b1391c687 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 e1bf1f9c3..e5ed5df70 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 d2631254b..a5d92002b 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 ac157d7ad..044ea1c6b 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, & @@ -385,6 +391,12 @@ module constants K_TRACKLENGTH = 3, & LEAKAGE = 4 + ! Differential tally independent variables + integer, parameter :: & + DIFF_DENSITY = 1, & + DIFF_NUCLIDE_DENSITY = 2, & + DIFF_TEMPERATURE = 3 + ! ============================================================================ ! RANDOM NUMBER STREAM CONSTANTS diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 79ac84a54..39552a716 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -7,7 +7,7 @@ module cross_section use global use list_header, only: ListElemInt use material_header, only: Material - use math, only: faddeeva, broaden_wmp_polynomials + use math, only: faddeeva, w_derivative, broaden_wmp_polynomials use multipole_header, only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, RM_RF, & MLBW_RT, MLBW_RX, MLBW_RA, MLBW_RF, FIT_T, FIT_A,& FIT_F, MultipoleArray @@ -706,6 +706,94 @@ contains end if end subroutine multipole_eval +!=============================================================================== +! MULTIPOLE_DERIV_EVAL evaluates the windowed multipole equations for the +! derivative of cross sections in the resolved resonance regions with respect to +! temperature. +!=============================================================================== + + subroutine multipole_deriv_eval(multipole, E, sqrtkT, sigT, sigA, sigF) + type(MultipoleArray), intent(in) :: multipole ! The windowed multipole + ! object to process. + real(8), intent(in) :: E ! The energy at which to + ! evaluate the cross section + real(8), intent(in) :: sqrtkT ! The temperature in the form + ! sqrt(kT), at which to + ! evaluate the XS. + real(8), intent(out) :: sigT ! Total cross section + real(8), intent(out) :: sigA ! Absorption cross section + real(8), intent(out) :: sigF ! Fission cross section + complex(8) :: w_val ! The faddeeva function evaluated at Z + complex(8) :: Z ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole) + complex(8) :: sigT_factor(multipole % num_l) + real(8) :: sqrtE ! sqrt(E), eV + real(8) :: invE ! 1/E, eV + real(8) :: dopp ! sqrt(atomic weight ratio / kT) + integer :: i_pole ! index of pole + integer :: i_window ! index of window + integer :: startw ! window start pointer (for poles) + integer :: endw ! window end pointer + real(8) :: T + + ! ========================================================================== + ! Bookkeeping + + ! Define some frequently used variables. + sqrtE = sqrt(E) + invE = ONE / E + dopp = multipole % sqrtAWR / sqrtkT + T = sqrtkT**2 / K_BOLTZMANN + + if (sqrtkT == ZERO) call fatal_error("Windowed multipole temperature & + &derivatives are not implemented for 0 Kelvin cross sections.") + + ! Locate us + i_window = floor((sqrtE - sqrt(multipole % start_E)) / multipole % spacing & + + ONE) + startw = multipole % w_start(i_window) + endw = multipole % w_end(i_window) + + ! Fill in factors. + if (startw <= endw) then + call compute_sigT_factor(multipole, sqrtE, sigT_factor) + end if + + ! Initialize the ouptut cross sections. + sigT = ZERO + sigA = ZERO + sigF = ZERO + + ! TODO Polynomials: Some of the curvefit polynomials Doppler broaden so + ! rigorously we should be computing the derivative of those. But in + ! practice, those derivatives are only large at very low energy and they + ! have no effect on reactor calculations. + + ! ========================================================================== + ! Add the contribution from the poles in this window. + + if (endw >= startw) then + do i_pole = startw, endw + Z = (sqrtE - multipole % data(MP_EA, i_pole)) * dopp + w_val = -invE * SQRT_PI * HALF * w_derivative(Z, 2) + if (multipole % formalism == FORM_MLBW) then + sigT = sigT + real((multipole % data(MLBW_RT, i_pole) * & + sigT_factor(multipole%l_value(i_pole)) + & + multipole % data(MLBW_RX, i_pole)) * w_val) + sigA = sigA + real(multipole % data(MLBW_RA, i_pole) * w_val) + sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * w_val) + else if (multipole % formalism == FORM_RM) then + sigT = sigT + real(multipole % data(RM_RT, i_pole) * w_val * & + sigT_factor(multipole % l_value(i_pole))) + sigA = sigA + real(multipole % data(RM_RA, i_pole) * w_val) + sigF = sigF + real(multipole % data(RM_RF, i_pole) * w_val) + end if + end do + sigT = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigT + sigA = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigA + sigF = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigF + end if + end subroutine multipole_deriv_eval + !=============================================================================== ! COMPUTE_SIGT_FACTOR calculates the sigT_factor, a factor inside of the sigT ! equation not present in the sigA and sigF equations. diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 713bbc351..05c8190a6 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 86100d195..83a5ac5fa 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 8111fd514..02f687018 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, 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 @@ -183,6 +185,10 @@ module global integer :: n_tallies = 0 ! # of tallies integer :: n_user_tallies = 0 ! # of user tallies + ! Tally derivatives + type(TallyDerivative), allocatable :: tally_derivs(:) +!$omp threadprivate(tally_derivs) + ! Normalization for statistics integer :: n_realizations = 0 ! # of independent realizations real(8) :: total_weight ! total starting particle weight in realization @@ -272,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 2a0711954..5e14605ae 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 4fe6f8769..692319acb 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 ff530c618..1bbf68f48 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 @@ -2629,11 +2628,13 @@ contains type(Node), pointer :: node_mesh => null() type(Node), pointer :: node_tal => null() type(Node), pointer :: node_filt => null() - type(Node), pointer :: node_trigger=>null() + type(Node), pointer :: node_trigger => null() + type(Node), pointer :: node_deriv => null() type(NodeList), pointer :: node_mesh_list => null() type(NodeList), pointer :: node_tal_list => null() type(NodeList), pointer :: node_filt_list => null() type(NodeList), pointer :: node_trigger_list => null() + type(NodeList), pointer :: node_deriv_list => null() type(ElemKeyValueCI), pointer :: scores type(ElemKeyValueCI), pointer :: next @@ -2641,6 +2642,12 @@ contains filename = trim(path_input) // "tallies.xml" inquire(FILE=filename, EXIST=file_exists) if (.not. file_exists) then + ! We need to allocate tally_derivs to avoid segfaults. Also needs to be + ! done in parallel because tally derivs are threadprivate. +!$omp parallel + allocate(tally_derivs(0)) +!$omp end parallel + ! Since a tallies.xml file is optional, no error is issued here return end if @@ -2819,6 +2826,94 @@ contains ! We only need the mesh info for plotting if (run_mode == MODE_PLOTTING) return + ! ========================================================================== + ! READ DATA FOR DERIVATIVES + + ! Get pointer list to XML nodes and allocate global array. + ! The array is threadprivate so it must be allocated in parallel. + call get_node_list(doc, "derivative", node_deriv_list) +!$omp parallel + allocate(tally_derivs(get_list_size(node_deriv_list))) +!$omp end parallel + + ! Make sure this is not an MG run. + if (.not. run_CE .and. get_list_size(node_deriv_list) > 0) then + call fatal_error("Differential tallies not supported in multi-group mode") + end if + + ! Read derivative attributes. + do i = 1, get_list_size(node_deriv_list) + associate(deriv => tally_derivs(i)) + ! Get pointer to derivative node. + call get_list_item(node_deriv_list, i, node_deriv) + + ! Copy the derivative id. + if (check_for_node(node_deriv, "id")) then + call get_node_value(node_deriv, "id", deriv % id) + else + call fatal_error("Must specify an ID for elements in the& + & tally XML file") + end if + + ! Make sure the id is > 0. + if (deriv % id <= 0) then + call fatal_error(" IDs must be an integer greater than & + &zero") + end if + + ! Make sure this id has not already been used. + do j = 1, i-1 + if (tally_derivs(j) % id == deriv % id) then + call fatal_error("Two or more 's use the same unique & + &ID: " // trim(to_str(deriv % id))) + end if + end do + + ! Read the independent variable name. + temp_str = "" + call get_node_value(node_deriv, "variable", temp_str) + temp_str = to_lower(temp_str) + + select case(temp_str) + + case("density") + deriv % variable = DIFF_DENSITY + call get_node_value(node_deriv, "material", deriv % diff_material) + + case("nuclide_density") + deriv % variable = DIFF_NUCLIDE_DENSITY + call get_node_value(node_deriv, "material", deriv % diff_material) + + call get_node_value(node_deriv, "nuclide", word) + word = trim(to_lower(word)) + pair_list => nuclide_dict % keys() + do while (associated(pair_list)) + if (starts_with(pair_list % key, word)) then + word = pair_list % key(1:150) + exit + end if + + ! Advance to next + pair_list => pair_list % next + end do + + ! Check if no nuclide was found + if (.not. associated(pair_list)) then + call fatal_error("Could not find the nuclide " & + // trim(word) // " specified in derivative " & + // trim(to_str(deriv % id)) // " in any material.") + end if + deallocate(pair_list) + + deriv % diff_nuclide = nuclide_dict % get_key(word) + + case("temperature") + deriv % variable = DIFF_TEMPERATURE + call get_node_value(node_deriv, "material", deriv % diff_material) + end select + end associate + end do + ! ========================================================================== ! READ TALLY DATA @@ -2840,7 +2935,7 @@ contains t % estimator = ESTIMATOR_TRACKLENGTH - ! Copy material id + ! Copy tally id if (check_for_node(node_tal, "id")) then call get_node_value(node_tal, "id", t % id) else @@ -3430,6 +3525,7 @@ contains call fatal_error("Cannot tally flux with an outgoing energy & &filter.") end if + case ('flux-yn') ! Prohibit user from tallying flux for an individual nuclide if (.not. (t % n_nuclide_bins == 1 .and. & @@ -3767,6 +3863,48 @@ contains // trim(to_str(t % id)) // ".") end if + ! Check for a tally derivative. + if (check_for_node(node_tal, "derivative")) then + ! Temporarily store the derivative id. + call get_node_value(node_tal, "derivative", t % deriv) + + ! Find the derivative with the given id, and store it's index. + do j = 1, size(tally_derivs) + if (tally_derivs(j) % id == t % deriv) then + t % deriv = j + ! Only analog or collision estimators are supported for differential + ! tallies. + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + t % estimator = ESTIMATOR_COLLISION + end if + ! We found the derivative we were looking for; exit the do loop. + exit + end if + if (j == size(tally_derivs)) then + call fatal_error("Could not find derivative " & + // trim(to_str(t % deriv)) // " specified on tally " & + // trim(to_str(t % id))) + end if + end do + + if (tally_derivs(t % deriv) % variable == DIFF_NUCLIDE_DENSITY & + .or. tally_derivs(t % deriv) % variable == DIFF_TEMPERATURE) then + if (any(t % nuclide_bins == -1)) then + if (t % find_filter(FILTER_ENERGYOUT) > 0) then + call fatal_error("Error on tally " // trim(to_str(t % id)) & + // ": Cannot use a 'nuclide_density' or 'temperature' & + &derivative on a tally with an outgoing energy filter and & + &'total' nuclide rate. Instead, tally each nuclide in the & + &material individually.") + ! Note that diff tallies with these characteristics would work + ! correctly if no tally events occur in the perturbed material + ! (e.g. pertrubing moderator but only tallying fuel), but this + ! case would be hard to check for by only reading inputs. + end if + end if + end if + end if + ! If settings.xml trigger is turned on, create tally triggers if (trigger_on) then @@ -4788,8 +4926,8 @@ contains sum_percent = sum(mat % atom_density) mat % atom_density = mat % atom_density / sum_percent - ! Change density in g/cm^3 to atom/b-cm. Since all values are now in atom - ! percent, the sum needs to be re-evaluated as 1/sum(x*awr) + ! Change density in g/cm^3 to atom/b-cm. Since all values are now in + ! atom percent, the sum needs to be re-evaluated as 1/sum(x*awr) if (.not. density_in_atom) then sum_percent = ZERO do j = 1, mat % n_nuclides @@ -4808,6 +4946,18 @@ contains ! Calculate nuclide atom densities mat % atom_density = mat % density * mat % atom_density + + ! Calculate density in g/cm^3. + mat % density_gpcc = ZERO + do j = 1, mat % n_nuclides + if (run_CE) then + awr = nuclides(mat % nuclide(j)) % awr + else + awr = ONE + end if + mat % density_gpcc = mat % density_gpcc & + + mat % atom_density(j) * awr * MASS_NEUTRON / N_AVOGADRO + end do end associate end do diff --git a/src/material_header.F90 b/src/material_header.F90 index 772e3a415..f7c1b5db0 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -13,6 +13,7 @@ module material_header integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm + real(8) :: density_gpcc ! total density in g/cm^3 ! Energy grid information integer :: n_grid ! # of union material grid points diff --git a/src/math.F90 b/src/math.F90 index 3ab8195aa..eac4ea94c 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 @@ -752,6 +753,22 @@ contains end function faddeeva + recursive function w_derivative(z, order) result(wv) + complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at + integer, intent(in) :: order + complex(8) :: wv ! The resulting w(z) value + + select case(order) + case (0) + wv = faddeeva(z) + case (1) + wv = -TWO * z * faddeeva(z) + TWO * ONEI / SQRT_PI + case default + wv = -TWO * z * w_derivative(z, order-1) & + - TWO * (order-1) * w_derivative(z, order-2) + end select + end function w_derivative + !=============================================================================== ! BROADEN_WMP_POLYNOMIALS Doppler broadens the windowed multipole curvefit. The ! curvefit is a polynomial of the form diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index 4fb1f7272..0fe48f7c9 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -434,7 +434,7 @@ module mgxs_header integer :: ndims integer(HSIZE_T) :: dims(2) real(8), allocatable :: temp_arr(:), temp_2d(:, :) - real(8), allocatable :: temp_beta(:, :) + real(8), allocatable :: temp_beta(:, :), temp_3d(:, :, :) real(8) :: dmu, mu, norm, chi_sum integer :: order, order_dim, gin, gout, l, imu, length type(VectorInt) :: temps_to_read @@ -767,9 +767,57 @@ module mgxs_header ! If prompt-nu-fission present, set prompt-nu-fission if (object_exists(xsdata_grp, "prompt-nu-fission")) then - ! Set prompt-nu-fission - call read_dataset(xs % prompt_nu_fission, xsdata_grp, & - "prompt-nu-fission") + ! Get the dimensions of the prompt-nu-fission dataset + xsdata = open_dataset(xsdata_grp, "prompt-nu-fission") + call get_ndims(xsdata, ndims) + + ! If prompt-nu-fission is a vector + if (ndims == 1) then + + ! Set prompt_nu_fission + call read_dataset(xs % prompt_nu_fission, xsdata_grp, & + "prompt-nu-fission") + + ! If prompt-nu-fission is a matrix, set prompt_nu_fission and + ! chi_prompt. + else if (ndims == 2) then + + ! chi_prompt is embedded in prompt_nu_fission -> extract + ! chi_prompt + allocate(temp_arr(energy_groups * energy_groups)) + call read_dataset(temp_arr, xsdata_grp, "prompt-nu-fission") + allocate(temp_2d(energy_groups, energy_groups)) + temp_2d = reshape(temp_arr, (/energy_groups, energy_groups/)) + + ! Deallocate temporary 1D array for prompt_nu_fission matrix + deallocate(temp_arr) + + ! Set the vector prompt-nu-fission from the matrix + ! prompt-nu-fission + do gin = 1, energy_groups + xs % prompt_nu_fission(gin) = sum(temp_2d(:, gin)) + end do + + ! Now pull out information needed for chi + xs % chi_prompt(:, :) = temp_2d + + ! Deallocate temporary 2D array for nu_fission matrix + deallocate(temp_2d) + + ! Normalize chi so its CDF goes to 1 + do gin = 1, energy_groups + chi_sum = sum(xs % chi_prompt(:, gin)) + if (chi_sum == ZERO) then + call fatal_error("Encountered chi prompt for a group & + &that sums to zero") + else + xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / chi_sum + end if + end do + else + call fatal_error("prompt-nu-fission must be provided as a 1D & + &or 2D array") + end if end if ! If delayed-nu-fission provided, set delayed-nu-fission. If @@ -848,9 +896,52 @@ module mgxs_header ! Deallocate temporary array for delayed-nu-fission matrix deallocate(temp_arr) + ! If delayed nu-fission is a 3D matrix, set delayed_nu_fission + ! and chi_delayed. + else if (ndims == 3) then + + ! chi_delayed is embedded in delayed_nu_fission -> extract + ! chi_delayed + allocate(temp_arr(delayed_groups * energy_groups * & + energy_groups)) + call read_dataset(temp_arr, xsdata_grp, "delayed-nu-fission") + allocate(temp_3d(delayed_groups, energy_groups, energy_groups)) + temp_3d = reshape(temp_arr, (/delayed_groups, energy_groups, & + energy_groups/)) + + ! Deallocate temporary 1D array for delayed_nu_fission matrix + deallocate(temp_arr) + + ! Set the 2D delayed-nu-fission matrix and 3D chi_dealyed matrix + ! from the 3D delayed-nu-fission matrix + do dg = 1, delayed_groups + do gin = 1, energy_groups + xs % delayed_nu_fission(dg, gin) = sum(temp_3d(dg, :, gin)) + do gout = 1, energy_groups + xs % chi_delayed(dg, gout, gin) = temp_3d(dg, gout, gin) + end do + end do + end do + + ! Normalize chi_delayed so its CDF goes to 1 + do dg = 1, delayed_groups + do gin = 1, energy_groups + chi_sum = sum(xs % chi_delayed(dg, :, gin)) + if (chi_sum == ZERO) then + call fatal_error("Encountered chi delayed for a group & + &that sums to zero") + else + xs % chi_delayed(dg, :, gin) = & + xs % chi_delayed(dg, :, gin) / chi_sum + end if + end do + end do + + ! Deallocate temporary 3D matrix for delayed_nu_fission + deallocate(temp_3d) else call fatal_error("delayed-nu-fission must be provided as a & - &1D or 2D array") + &1D, 2D, or 3D array") end if end if @@ -1119,7 +1210,8 @@ module mgxs_header integer(HSIZE_T) :: dims(4) integer, allocatable :: int_arr(:) real(8), allocatable :: temp_1d(:), temp_3d(:, :, :) - real(8), allocatable :: temp_4d(:, :, :, :), temp_beta(:, :, :, :) + real(8), allocatable :: temp_4d(:, :, :, :), temp_5d(:, :, :, :, :) + real(8), allocatable :: temp_beta(:, :, :, :) real(8) :: dmu, mu, norm, chi_sum integer :: order, order_dim, gin, gout, l, imu, dg type(VectorInt) :: temps_to_read @@ -1540,16 +1632,70 @@ module mgxs_header ! If prompt-nu-fission present, set prompt-nu-fission if (object_exists(xsdata_grp, "prompt-nu-fission")) then - ! Allocate temporary array for prompt-nu-fission - allocate(temp_1d(energy_groups * this % n_azi * this % n_pol)) + ! Get the dimensions of the prompt-nu-fission dataset + xsdata = open_dataset(xsdata_grp, "prompt-nu-fission") + call get_ndims(xsdata, ndims) - ! Read prompt-nu-fission - call read_dataset(temp_1d, xsdata_grp, "prompt-nu-fission") - xs % prompt_nu_fission = reshape(temp_1d, (/energy_groups, & - this % n_azi, this % n_pol/)) + ! If prompt-nu-fission is a vector for each azi and pol + if (ndims == 3) then - ! Deallocate temporary array for prompt-nu-fission - deallocate(temp_1d) + ! Set prompt_nu_fission + call read_dataset(xs % prompt_nu_fission, xsdata_grp, & + "prompt-nu-fission") + + ! If prompt-nu-fission is a matrix for each azi and pol, + ! set prompt_nu_fission and chi_prompt. + else if (ndims == 4) then + + ! chi_prompt is embedded in prompt_nu_fission -> extract + ! chi_prompt + allocate(temp_1d(energy_groups * energy_groups & + * this % n_azi * this % n_pol)) + allocate(temp_4d(energy_groups, energy_groups, this % n_azi, & + this % n_pol)) + call read_dataset(temp_1d, xsdata_grp, "prompt-nu-fission") + temp_4d = reshape(temp_1d, (/energy_groups, energy_groups, & + this % n_azi, this % n_pol/)) + + ! Deallocate temporary 1D array for prompt_nu_fission matrix + deallocate(temp_1d) + + ! Set the vector prompt-nu-fission from the matrix + ! prompt-nu-fission + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gin = 1, energy_groups + xs % prompt_nu_fission(gin, iazi, ipol) = & + sum(temp_4d(:, gin, iazi, ipol)) + end do + end do + end do + + ! Now pull out information needed for chi + xs % chi_prompt(:, :, :, :) = temp_4d + + ! Deallocate temporary 4D array for nu_fission matrix + deallocate(temp_4d) + + ! Normalize chi so its CDF goes to 1 + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gin = 1, energy_groups + chi_sum = sum(xs % chi_prompt(:, gin, iazi, ipol)) + if (chi_sum == ZERO) then + call fatal_error("Encountered chi prompt for a group & + &that sums to zero") + else + xs % chi_prompt(:, gin, iazi, ipol) = & + xs % chi_prompt(:, gin, iazi, ipol) / chi_sum + end if + end do + end do + end do + else + call fatal_error("prompt-nu-fission must be provided as a 3D & + &or 4D array") + end if end if ! If delayed-nu-fission provided, set delayed-nu-fission. If @@ -1639,9 +1785,64 @@ module mgxs_header ! Deallocate temporary array for delayed-nu-fission matrix deallocate(temp_1d) + ! If delayed nu-fission is a 5D matrix, set delayed_nu_fission + ! and chi_delayed. + else if (ndims == 5) then + + ! chi_delayed is embedded in delayed_nu_fission -> extract + ! chi_delayed + allocate(temp_1d(delayed_groups * energy_groups * & + energy_groups * this % n_azi * this % n_pol)) + allocate(temp_5d(delayed_groups, energy_groups, energy_groups, & + this % n_azi, this % n_pol)) + call read_dataset(temp_1d, xsdata_grp, "delayed-nu-fission") + temp_5d = reshape(temp_1d, (/delayed_groups, energy_groups, & + energy_groups, this % n_azi, this % n_pol/)) + + ! Deallocate temporary 1D array for delayed_nu_fission matrix + deallocate(temp_1d) + + ! Set the 4D delayed-nu-fission matrix and 5D chi_delayed matrix + ! from the 5D delayed-nu-fission matrix + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do dg = 1, delayed_groups + do gin = 1, energy_groups + xs % delayed_nu_fission(dg, gin, iazi, ipol) = & + sum(temp_5d(dg, :, gin, iazi, ipol)) + do gout = 1, energy_groups + xs % chi_delayed(dg, gout, gin, iazi, ipol) = & + temp_5d(dg, gout, gin, iazi, ipol) + end do + end do + end do + end do + end do + + ! Normalize chi_delayed so its CDF goes to 1 + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do dg = 1, delayed_groups + do gin = 1, energy_groups + chi_sum = sum(xs % chi_delayed(dg, :, gin, iazi, ipol)) + if (chi_sum == ZERO) then + call fatal_error("Encountered chi delayed for a group& + & that sums to zero") + else + xs % chi_delayed(dg, :, gin, iazi, ipol) = & + xs % chi_delayed(dg, :, gin, iazi, ipol) / & + chi_sum + end if + end do + end do + end do + end do + + ! Deallocate temporary 5D matrix for delayed_nu_fission + deallocate(temp_5d) else call fatal_error("delayed-nu-fission must be provided as a & - &1D or 2D array") + &3D, 4D, or 5D array") end if end if diff --git a/src/output.F90 b/src/output.F90 index 8127725a6..f5766187e 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 @@ -777,6 +777,28 @@ contains // trim(t % name), unit=unit_tally, level=3) endif + ! Write derivative information. + if (t % deriv /= NONE) then + associate(deriv => tally_derivs(t % deriv)) + select case (deriv % variable) + case (DIFF_DENSITY) + write(unit=unit_tally, fmt="(' Density derivative Material ',A)") & + to_str(deriv % diff_material) + case (DIFF_NUCLIDE_DENSITY) + write(unit=unit_tally, fmt="(' Nuclide density derivative & + &Material ',A,' Nuclide ',A)") & + trim(to_str(deriv % diff_material)), & + trim(nuclides(deriv % diff_nuclide) % name) + case (DIFF_TEMPERATURE) + write(unit=unit_tally, fmt="(' Temperature derivative Material ',& + &A)") to_str(deriv % diff_material) + case default + call fatal_error("Differential tally dependent variable for tally "& + // trim(to_str(t % id)) // " not defined in output.F90.") + end select + end associate + end if + ! Handle surface current tallies separately if (t % type == TALLY_SURFACE_CURRENT) then call write_surface_current(t, unit_tally) @@ -876,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) @@ -886,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, & @@ -902,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 @@ -916,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 @@ -1020,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 @@ -1037,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 @@ -1056,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 @@ -1073,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 @@ -1092,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 @@ -1109,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 776c7905a..082b337a4 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/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index cc327f804..2143c8e47 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -15,6 +15,23 @@ element tallies { ) }* & + element derivative { + (element id { xsd:int } | attribute id { xsd:int }) & + (element material { xsd:int } | attribute material { xsd:int }) & + ( (element variable { ( "density") } + | attribute variable { ( "density" ) } ) | + ( + (element variable { ( "nuclide_density" ) } + | attribute variable { ( "nuclide_density" ) } ) + & + (element nuclide { xsd:string { maxLength = "12" } } + | attribute nuclide { xsd:string { maxLength = "12" } } ) | + ) + (element variable { ( "temperature") } + | attribute variable { ( "temperature" ) } ) + ) + }* & + element tally { (element id { xsd:int } | attribute id { xsd:int }) & (element name { xsd:string { maxLength="52" } } | @@ -41,7 +58,8 @@ element tallies { (element type { xsd:string } | attribute type { xsd:string }) & (element threshold { xsd:double} | attribute threshold { xsd:double }) & (element scores { list { xsd:string { maxLength = "20" }+ } } | attribute scores { list { xsd:string { maxLength = "20"}+ } } )? - }* + }* & + (element derivative { xsd:int } | attribute derivative { xsd:int } )? }* & element assume_separate { xsd:boolean }? diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index 755e9e90d..fc60b8373 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -89,6 +89,68 @@ + + + + + + + + + + + + + + + + + + + + + + + density + + + density + + + + + + nuclide_density + + + nuclide_density + + + + + + 12 + + + + + 12 + + + + + + + temperature + + + temperature + + + + + + @@ -254,6 +316,16 @@ + + + + + + + + + + diff --git a/src/simulation.F90 b/src/simulation.F90 index f59ce371c..e6b20ddf8 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 @@ -84,7 +83,7 @@ contains ! ==================================================================== ! LOOP OVER PARTICLES -!$omp parallel do schedule(static) firstprivate(p) +!$omp parallel do schedule(static) firstprivate(p) copyin(tally_derivs) PARTICLE_LOOP: do i_work = 1, work current_work = i_work @@ -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 6b7a467d4..23543775c 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -50,7 +50,8 @@ contains integer, allocatable :: key_array(:) integer(HID_T) :: file_id integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, & - mesh_group, filter_group, runtime_group + mesh_group, filter_group, derivs_group, deriv_group, & + runtime_group character(MAX_WORD_LEN), allocatable :: str_array(:) character(MAX_FILE_LEN) :: filename type(RegularMesh), pointer :: meshp @@ -198,6 +199,38 @@ contains call close_group(meshes_group) + ! Write information for derivatives. + if (size(tally_derivs) > 0) then + derivs_group = create_group(tallies_group, "derivatives") + do i = 1, size(tally_derivs) + associate(deriv => tally_derivs(i)) + deriv_group = create_group(derivs_group, "derivative " & + // trim(to_str(deriv % id))) + select case (deriv % variable) + case (DIFF_DENSITY) + call write_dataset(deriv_group, "independent variable", "density") + call write_dataset(deriv_group, "material", deriv % diff_material) + case (DIFF_NUCLIDE_DENSITY) + call write_dataset(deriv_group, "independent variable", & + "nuclide_density") + call write_dataset(deriv_group, "material", deriv % diff_material) + call write_dataset(deriv_group, "nuclide", & + nuclides(deriv % diff_nuclide) % name) + case (DIFF_TEMPERATURE) + call write_dataset(deriv_group, "independent variable", & + "temperature") + call write_dataset(deriv_group, "material", deriv % diff_material) + case default + call fatal_error("Independent variable for derivative " & + // trim(to_str(deriv % id)) // " not defined in & + &state_point.F90.") + end select + call close_group(deriv_group) + end associate + end do + call close_group(derivs_group) + end if + ! Write number of tallies call write_dataset(tallies_group, "n_tallies", n_tallies) @@ -273,6 +306,13 @@ contains call write_dataset(tally_group, "nuclides", str_array) deallocate(str_array) + ! Write derivative information. + if (tally % deriv /= NONE) then + call write_dataset(tally_group, "derivative", & + tally_derivs(tally % deriv) % id) + end if + + ! Write scores. call write_dataset(tally_group, "n_score_bins", tally % n_score_bins) allocate(str_array(size(tally % score_bins))) do j = 1, size(tally % score_bins) @@ -353,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 @@ -481,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 @@ -489,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 @@ -505,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 @@ -519,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 @@ -568,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 " // & @@ -592,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 @@ -771,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") @@ -787,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 69e4e8ecd..79f22610e 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1,11 +1,14 @@ module tally + use, intrinsic :: ISO_C_BINDING + #ifdef MPI use message_passing #endif use algorithm, only: binary_search use constants + use cross_section, only: multipole_deriv_eval use error, only: fatal_error use geometry_header use global @@ -18,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 @@ -1086,6 +1088,14 @@ contains end select + !######################################################################### + ! Add derivative information on score for differential tallies. + + if (t % deriv /= NONE) then + call apply_derivative_to_score(p, t, i_nuclide, atom_density, & + score_bin, score) + end if + !######################################################################### ! Expand score if necessary and add to tally results. call expand_and_score(p, t, score_index, filter_index, score_bin, & @@ -1965,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) @@ -1982,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 @@ -2011,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 @@ -2031,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) @@ -2040,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 @@ -2421,6 +2429,13 @@ contains ! determine score based on bank site weight and keff score = keff * fission_bank(n_bank - p % n_bank + k) % wgt + ! Add derivative information for differential tallies. Note that the + ! i_nuclide and atom_density arguments do not matter since this is an + ! analog estimator. + if (t % deriv /= NONE) then + call apply_derivative_to_score(p, t, 0, ZERO, SCORE_NU_FISSION, score) + end if + if (.not. run_CE .and. eo_filt % matches_transport_groups) then ! determine outgoing energy from fission bank @@ -2458,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 @@ -2498,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 @@ -2536,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 @@ -2997,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 @@ -3033,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 @@ -3053,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 @@ -3089,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 @@ -3108,6 +3123,778 @@ contains end subroutine score_surface_current +!=============================================================================== +! APPLY_DERIVATIVE_TO_SCORE multiply the given score by its relative derivative +!=============================================================================== + + subroutine apply_derivative_to_score(p, t, i_nuclide, atom_density, & + score_bin, score) + type(Particle), intent(in) :: p + type(TallyObject), intent(in) :: t + integer, intent(in) :: i_nuclide + real(8), intent(in) :: atom_density ! atom/b-cm + integer, intent(in) :: score_bin + real(8), intent(inout) :: score + + integer :: l + logical :: scoring_diff_nuclide + real(8) :: flux_deriv + real(8) :: dsigT, dsigA, dsigF, cum_dsig + + if (score == ZERO) return + + ! If our score was previously c then the new score is + ! c * (1/f * d_f/d_p + 1/c * d_c/d_p) + ! where (1/f * d_f/d_p) is the (logarithmic) flux derivative and p is the + ! perturbated variable. + + associate(deriv => tally_derivs(t % deriv)) + flux_deriv = deriv % flux_deriv + + select case (tally_derivs(t % deriv) % variable) + + !========================================================================= + ! Density derivative: + ! c = Sigma_MT + ! c = sigma_MT * N + ! c = sigma_MT * rho * const + ! d_c / d_rho = sigma_MT * const + ! (1 / c) * (d_c / d_rho) = 1 / rho + + case (DIFF_DENSITY) + select case (t % estimator) + + case (ESTIMATOR_ANALOG) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & + SCORE_NU_FISSION) + if (materials(p % material) % id == deriv % diff_material) then + score = score * (flux_deriv + ONE & + / materials(p % material) % density_gpcc) + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case (ESTIMATOR_COLLISION) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & + SCORE_NU_FISSION) + if (materials(p % material) % id == deriv % diff_material) then + score = score * (flux_deriv + ONE & + / materials(p % material) % density_gpcc) + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case default + call fatal_error("Differential tallies are only implemented for & + &analog and collision estimators.") + end select + + !========================================================================= + ! Nuclide density derivative: + ! If we are scoring a reaction rate for a single nuclide then + ! c = Sigma_MT_i + ! c = sigma_MT_i * N_i + ! d_c / d_N_i = sigma_MT_i + ! (1 / c) * (d_c / d_N_i) = 1 / N_i + ! If the score is for the total material (i_nuclide = -1) + ! c = Sum_i(Sigma_MT_i) + ! d_c / d_N_i = sigma_MT_i + ! (1 / c) * (d_c / d_N) = sigma_MT_i / Sigma_MT + ! where i is the perturbed nuclide. + + case (DIFF_NUCLIDE_DENSITY) + select case (t % estimator) + + case (ESTIMATOR_ANALOG) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & + SCORE_NU_FISSION) + if (materials(p % material) % id == deriv % diff_material & + .and. p % event_nuclide == deriv % diff_nuclide) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == deriv % diff_nuclide) exit + end do + + score = score * (flux_deriv & + + ONE / mat % atom_density(l)) + end associate + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case (ESTIMATOR_COLLISION) + scoring_diff_nuclide = & + (materials(p % material) % id == deriv % diff_material) & + .and. (i_nuclide == deriv % diff_nuclide) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % total /= ZERO) then + score = score * (flux_deriv & + + micro_xs(deriv % diff_nuclide) % total & + / material_xs % total) + else if (scoring_diff_nuclide .and. & + micro_xs(deriv % diff_nuclide) % total /= ZERO) then + score = score * (flux_deriv + ONE / atom_density) + else + score = score * flux_deriv + end if + + case (SCORE_SCATTER) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % total - material_xs % absorption /= ZERO) then + score = score * (flux_deriv & + + (micro_xs(deriv % diff_nuclide) % total & + - micro_xs(deriv % diff_nuclide) % absorption) & + / (material_xs % total - material_xs % absorption)) + else if (scoring_diff_nuclide .and. & + (micro_xs(deriv % diff_nuclide) % total & + - micro_xs(deriv % diff_nuclide) % absorption) /= ZERO) then + score = score * (flux_deriv + ONE / atom_density) + else + score = score * flux_deriv + end if + + case (SCORE_ABSORPTION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % absorption /= ZERO) then + score = score * (flux_deriv & + + micro_xs(deriv % diff_nuclide) % absorption & + / material_xs % absorption ) + else if (scoring_diff_nuclide .and. & + micro_xs(deriv % diff_nuclide) % absorption /= ZERO) then + score = score * (flux_deriv + ONE / atom_density) + else + score = score * flux_deriv + end if + + case (SCORE_FISSION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % fission /= ZERO) then + score = score * (flux_deriv & + + micro_xs(deriv % diff_nuclide) % fission & + / material_xs % fission) + else if (scoring_diff_nuclide .and. & + micro_xs(deriv % diff_nuclide) % fission /= ZERO) then + score = score * (flux_deriv + ONE / atom_density) + else + score = score * flux_deriv + end if + + case (SCORE_NU_FISSION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % nu_fission /= ZERO) then + score = score * (flux_deriv & + + micro_xs(deriv % diff_nuclide) % nu_fission & + / material_xs % nu_fission) + else if (scoring_diff_nuclide .and. & + micro_xs(deriv % diff_nuclide) % nu_fission /= ZERO) then + score = score * (flux_deriv + ONE / atom_density) + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case default + call fatal_error("Differential tallies are only implemented for & + &analog and collision estimators.") + end select + + !========================================================================= + ! Temperature derivative: + ! If we are scoring a reaction rate for a single nuclide then + ! c = Sigma_MT_i + ! c = sigma_MT_i * N_i + ! d_c / d_T = (d_sigma_Mt_i / d_T) * N_i + ! (1 / c) * (d_c / d_T) = (d_sigma_MT_i / d_T) / sigma_MT_i + ! If the score is for the total material (i_nuclide = -1) + ! (1 / c) * (d_c / d_T) = Sum_i((d_sigma_MT_i / d_T) * N_i) / Sigma_MT_i + ! where i is the perturbed nuclide. The d_sigma_MT_i / d_T term is + ! computed by multipole_deriv_eval. It only works for the resolved + ! resonance range and requires multipole data. + + case (DIFF_TEMPERATURE) + select case (t % estimator) + + case (ESTIMATOR_ANALOG) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL) + if (materials(p % material) % id == deriv % diff_material .and. & + micro_xs(p % event_nuclide) % total > ZERO) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == p % event_nuclide) exit + end do + + dsigT = ZERO + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigT * mat % atom_density(l) / material_xs % total) + end associate + else + score = score * flux_deriv + end if + + case (SCORE_SCATTER) + if (materials(p % material) % id == deriv % diff_material .and. & + (micro_xs(p % event_nuclide) % total & + - micro_xs(p % event_nuclide) % absorption) > ZERO) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == p % event_nuclide) exit + end do + + dsigT = ZERO + dsigA = ZERO + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv + (dsigT - dsigA) & + * mat % atom_density(l) / & + (material_xs % total - material_xs % absorption)) + end associate + else + score = score * flux_deriv + end if + + case (SCORE_ABSORPTION) + if (materials(p % material) % id == deriv % diff_material .and. & + micro_xs(p % event_nuclide) % absorption > ZERO) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == p % event_nuclide) exit + end do + + dsigA = ZERO + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigA * mat % atom_density(l) / material_xs % absorption) + end associate + else + score = score * flux_deriv + end if + + case (SCORE_FISSION) + if (materials(p % material) % id == deriv % diff_material .and. & + micro_xs(p % event_nuclide) % fission > ZERO) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == p % event_nuclide) exit + end do + + dsigF = ZERO + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigF * mat % atom_density(l) / material_xs % fission) + end associate + else + score = score * flux_deriv + end if + + case (SCORE_NU_FISSION) + if (materials(p % material) % id == deriv % diff_material .and. & + micro_xs(p % event_nuclide) % nu_fission > ZERO) then + associate(mat => materials(p % material)) + ! Search for the index of the perturbed nuclide. + do l = 1, mat % n_nuclides + if (mat % nuclide(l) == p % event_nuclide) exit + end do + + dsigF = ZERO + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigF * mat % atom_density(l) / material_xs % nu_fission& + * micro_xs(p % event_nuclide) % nu_fission & + / micro_xs(p % event_nuclide) % fission) + end associate + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case (ESTIMATOR_COLLISION) + + select case (score_bin) + + case (SCORE_FLUX) + score = score * flux_deriv + + case (SCORE_TOTAL) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % total > ZERO) then + cum_dsig = ZERO + associate(mat => materials(p % material)) + do l = 1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E .and. & + micro_xs(mat % nuclide(l)) % total > ZERO) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + cum_dsig = cum_dsig + dsigT * mat % atom_density(l) + end if + end associate + end do + end associate + score = score * (flux_deriv & + + cum_dsig / material_xs % total) + else if (materials(p % material) % id == deriv % diff_material & + .and. material_xs % total > ZERO) then + dsigT = ZERO + associate (nuc => nuclides(i_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigT / micro_xs(i_nuclide) % total) + else + score = score * flux_deriv + end if + + case (SCORE_SCATTER) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + (material_xs % total - material_xs % absorption) > ZERO) then + cum_dsig = ZERO + associate(mat => materials(p % material)) + do l = 1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E .and. & + (micro_xs(mat % nuclide(l)) % total & + - micro_xs(mat % nuclide(l)) % absorption) > ZERO) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + cum_dsig = cum_dsig & + + (dsigT - dsigA) * mat % atom_density(l) + end if + end associate + end do + end associate + score = score * (flux_deriv + cum_dsig & + / (material_xs % total - material_xs % absorption)) + else if ( materials(p % material) % id == deriv % diff_material & + .and. (material_xs % total - material_xs % absorption) > ZERO)& + then + dsigT = ZERO + dsigA = ZERO + associate (nuc => nuclides(i_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv + (dsigT - dsigA) & + / (micro_xs(i_nuclide) % total & + - micro_xs(i_nuclide) % absorption)) + else + score = score * flux_deriv + end if + + case (SCORE_ABSORPTION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % absorption > ZERO) then + cum_dsig = ZERO + associate(mat => materials(p % material)) + do l = 1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E .and. & + micro_xs(mat % nuclide(l)) % absorption > ZERO) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + cum_dsig = cum_dsig + dsigA * mat % atom_density(l) + end if + end associate + end do + end associate + score = score * (flux_deriv & + + cum_dsig / material_xs % absorption) + else if (materials(p % material) % id == deriv % diff_material & + .and. material_xs % absorption > ZERO) then + dsigA = ZERO + associate (nuc => nuclides(i_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigA / micro_xs(i_nuclide) % absorption) + else + score = score * flux_deriv + end if + + case (SCORE_FISSION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % fission > ZERO) then + cum_dsig = ZERO + associate(mat => materials(p % material)) + do l = 1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E .and. & + micro_xs(mat % nuclide(l)) % fission > ZERO) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + cum_dsig = cum_dsig + dsigF * mat % atom_density(l) + end if + end associate + end do + end associate + score = score * (flux_deriv & + + cum_dsig / material_xs % fission) + else if (materials(p % material) % id == deriv % diff_material & + .and. material_xs % fission > ZERO) then + dsigF = ZERO + associate (nuc => nuclides(i_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigF / micro_xs(i_nuclide) % fission) + else + score = score * flux_deriv + end if + + case (SCORE_NU_FISSION) + if (i_nuclide == -1 .and. & + materials(p % material) % id == deriv % diff_material .and. & + material_xs % nu_fission > ZERO) then + cum_dsig = ZERO + associate(mat => materials(p % material)) + do l = 1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E .and. & + micro_xs(mat % nuclide(l)) % nu_fission > ZERO) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + cum_dsig = cum_dsig + dsigF * mat % atom_density(l) & + * micro_xs(mat % nuclide(l)) % nu_fission & + / micro_xs(mat % nuclide(l)) % fission + end if + end associate + end do + end associate + score = score * (flux_deriv & + + cum_dsig / material_xs % nu_fission) + else if (materials(p % material) % id == deriv % diff_material & + .and. material_xs % nu_fission > ZERO) then + dsigF = ZERO + associate (nuc => nuclides(i_nuclide)) + if (nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + end if + end associate + score = score * (flux_deriv & + + dsigF / micro_xs(i_nuclide) % fission) + else + score = score * flux_deriv + end if + + case default + call fatal_error('Tally derivative not defined for a score on & + &tally ' // trim(to_str(t % id))) + end select + + case default + call fatal_error("Differential tallies are only implemented for & + &analog and collision estimators.") + end select + end select + end associate + end subroutine apply_derivative_to_score + +!=============================================================================== +! SCORE_TRACK_DERIVATIVE Adjust flux derivatives on differential tallies to +! account for a neutron travelling through a perturbed material. +!=============================================================================== + + subroutine score_track_derivative(p, distance) + type(Particle), intent(in) :: p + real(8), intent(in) :: distance ! Neutron flight distance + + integer :: i, l + real(8) :: dsigT, dsigA, dsigF + + ! A void material cannot be perturbed so it will not affect flux derivatives + if (p % material == MATERIAL_VOID) return + + do i = 1, size(tally_derivs) + associate(deriv => tally_derivs(i)) + select case (deriv % variable) + + case (DIFF_DENSITY) + associate (mat => materials(p % material)) + if (mat % id == deriv % diff_material) then + ! phi is proportional to e^(-Sigma_tot * dist) + ! (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist + ! (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist + deriv % flux_deriv = deriv % flux_deriv & + - distance * material_xs % total / mat % density_gpcc + end if + end associate + + case (DIFF_NUCLIDE_DENSITY) + associate (mat => materials(p % material)) + if (mat % id == deriv % diff_material) then + ! phi is proportional to e^(-Sigma_tot * dist) + ! (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist + ! (1 / phi) * (d_phi / d_N) = - sigma_tot * dist + deriv % flux_deriv = deriv % flux_deriv & + - distance * micro_xs(deriv % diff_nuclide) % total + end if + end associate + + case (DIFF_TEMPERATURE) + associate (mat => materials(p % material)) + if (mat % id == deriv % diff_material) then + do l=1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (nuc % mp_present .and. & + p % E >= nuc % multipole % start_E .and. & + p % E <= nuc % multipole % end_E) then + ! phi is proportional to e^(-Sigma_tot * dist) + ! (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist + ! (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist + call multipole_deriv_eval(nuc % multipole, p % E, & + p % sqrtkT, dsigT, dsigA, dsigF) + deriv % flux_deriv = deriv % flux_deriv & + - distance * dsigT * mat % atom_density(l) + end if + end associate + end do + end if + end associate + end select + end associate + end do + end subroutine score_track_derivative + +!=============================================================================== +! SCORE_COLLISION_DERIVATIVE Adjust flux derivatives on differential tallies to +! account for a neutron scattering in the perturbed material. Note that this +! subroutine will be called after absorption events in addition to scattering +! events, but any flux derivatives scored after an absorption will never be +! tallied. This is because the order of operations is +! 1. Particle is moved. +! 2. score_track_derivative is called. +! 3. Collision physics are computed, and the particle is labeled absorbed. +! 4. Analog- and collision-estimated tallies are scored. +! 5. This subroutine is called. +! 6. Particle is killed and no more tallies are scored. +! Hence, it is safe to assume that only derivative of the scattering cross +! section need to be computed here. +!=============================================================================== + + subroutine score_collision_derivative(p) + type(Particle), intent(in) :: p + + integer :: i, j, l + real(8) :: dsigT, dsigA, dsigF + + ! A void material cannot be perturbed so it will not affect flux derivatives + if (p % material == MATERIAL_VOID) return + + do i = 1, size(tally_derivs) + associate(deriv => tally_derivs(i)) + select case (deriv % variable) + + case (DIFF_DENSITY) + associate (mat => materials(p % material)) + if (mat % id == deriv % diff_material) then + ! phi is proportional to Sigma_s + ! (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s + ! (1 / phi) * (d_phi / d_rho) = 1 / rho + deriv % flux_deriv = deriv % flux_deriv & + + ONE / mat % density_gpcc + end if + end associate + + case (DIFF_NUCLIDE_DENSITY) + associate (mat => materials(p % material)) + if (mat % id == deriv % diff_material & + .and. p % event_nuclide == deriv % diff_nuclide) then + ! Find the index in this material for the diff_nuclide. + do j = 1, mat % n_nuclides + if (mat % nuclide(j) == deriv % diff_nuclide) exit + end do + ! Make sure we found the nuclide. + if (mat % nuclide(j) /= deriv % diff_nuclide) then + call fatal_error("Couldn't find the right nuclide.") + end if + ! phi is proportional to Sigma_s + ! (1 / phi) * (d_phi / d_N) = (d_Sigma_s / d_N) / Sigma_s + ! (1 / phi) * (d_phi / d_N) = sigma_s / Sigma_s + ! (1 / phi) * (d_phi / d_N) = 1 / N + deriv % flux_deriv = deriv % flux_deriv & + + ONE / mat % atom_density(j) + end if + end associate + + case (DIFF_TEMPERATURE) + associate (mat => materials(p % material)) + if (mat % id == deriv % diff_material) then + do l=1, mat % n_nuclides + associate (nuc => nuclides(mat % nuclide(l))) + if (mat % nuclide(l) == p % event_nuclide .and. & + nuc % mp_present .and. & + p % last_E >= nuc % multipole % start_E .and. & + p % last_E <= nuc % multipole % end_E) then + ! phi is proportional to Sigma_s + ! (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s + ! (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s + call multipole_deriv_eval(nuc % multipole, p % last_E, & + p % sqrtkT, dsigT, dsigA, dsigF) + deriv % flux_deriv = deriv % flux_deriv + (dsigT - dsigA)& + / (micro_xs(mat % nuclide(l)) % total & + - micro_xs(mat % nuclide(l)) % absorption) + ! Note that this is an approximation! The real scattering + ! cross section is Sigma_s(E'->E, uvw'->uvw) = + ! Sigma_s(E') * P(E'->E, uvw'->uvw). We are assuming that + ! d_P(E'->E, uvw'->uvw) / d_T = 0 and only computing + ! d_S(E') / d_T. Using this approximation in the vicinity + ! of low-energy resonances causes errors (~2-5% for PWR + ! pincell eigenvalue derivatives). + end if + end associate + end do + end if + end associate + end select + end associate + end do + end subroutine score_collision_derivative + +!=============================================================================== +! ZERO_FLUX_DERIVS Set the flux derivatives on differential tallies to zero. +!=============================================================================== + + subroutine zero_flux_derivs() + integer :: i + do i = 1, size(tally_derivs) + tally_derivs(i) % flux_deriv = ZERO + end do + end subroutine zero_flux_derivs + !=============================================================================== ! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history ! within the batch to a new random variable @@ -3116,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 @@ -3142,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 @@ -3152,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 @@ -3182,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 @@ -3191,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 @@ -3244,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 @@ -3251,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 fe385458b..1d34779ae 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -1,22 +1,27 @@ 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 +! TALLYDERIVATIVE describes a first-order derivative that can be applied to +! tallies. !=============================================================================== - type, bind(C) :: TallyResult - real(C_DOUBLE) :: value = 0. - real(C_DOUBLE) :: sum = 0. - real(C_DOUBLE) :: sum_sq = 0. - end type TallyResult + type TallyDerivative + integer :: id + integer :: variable + integer :: diff_material + integer :: diff_nuclide + real(8) :: flux_deriv + end type TallyDerivative !=============================================================================== ! TALLYOBJECT describes a user-specified tally. The region of phase space to @@ -68,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. @@ -79,6 +84,82 @@ module tally_header ! Tally precision triggers integer :: n_triggers = 0 ! # of triggers type(TriggerObject), allocatable :: triggers(:) ! Array of triggers + + ! 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 1ae403bd6..0e7a50a9c 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/tracking.F90 b/src/tracking.F90 index b62cfc85a..a45f8e448 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -14,7 +14,9 @@ module tracking use random_lcg, only: prn use string, only: to_str use tally, only: score_analog_tally, score_tracklength_tally, & - score_collision_tally, score_surface_current + score_collision_tally, score_surface_current, & + score_track_derivative, & + score_collision_derivative, zero_flux_derivs use track_output, only: initialize_particle_track, write_particle_track, & add_particle_track, finalize_particle_track @@ -63,6 +65,9 @@ contains call initialize_particle_track() endif + ! Every particle starts with no accumulated flux derivative. + if (active_tallies % size() > 0) call zero_flux_derivs() + EVENT_LOOP: do ! If the cell hasn't been determined based on the particle's location, ! initiate a search for the current cell. This generally happens at the @@ -137,6 +142,9 @@ contains distance * material_xs % nu_fission end if + ! Score flux derivative accumulators for differential tallies. + if (active_tallies % size() > 0) call score_track_derivative(p, distance) + if (d_collision > d_boundary) then ! ==================================================================== ! PARTICLE CROSSES SURFACE @@ -213,6 +221,9 @@ contains p % coord(j + 1) % uvw = p % coord(j) % uvw end if end do + + ! Score flux derivative accumulators for differential tallies. + if (active_tallies % size() > 0) call score_collision_derivative(p) end if ! Save coordinates for tallying purposes diff --git a/src/trigger.F90 b/src/trigger.F90 index 39786809b..f41cc46cf 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 diff --git a/tests/test_diff_tally/inputs_true.dat b/tests/test_diff_tally/inputs_true.dat new file mode 100644 index 000000000..9a6747579 --- /dev/null +++ b/tests/test_diff_tally/inputs_true.dat @@ -0,0 +1 @@ +df5a0ee7d353fe25344496058cea42e3244a7ab32c13c34c03b4b2f6adf88579f08a71ec0d22d8d264e50e3663717068ff54308e2d16aed8e085c222a8c2fdbd \ No newline at end of file diff --git a/tests/test_diff_tally/results_true.dat b/tests/test_diff_tally/results_true.dat new file mode 100644 index 000000000..7c4f7f776 --- /dev/null +++ b/tests/test_diff_tally/results_true.dat @@ -0,0 +1,177 @@ +d_material,d_nuclide,d_variable,score,mean,std. dev. +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,3.9902949e-02,9.1258428e-03 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,9.8213745e-04,9.8213745e-04 +1,,density,nu-fission,2.7945389e-02,2.0999368e-02 +1,,density,scatter,4.0626155e-01,1.7017674e-02 +1,,density,nu-fission,2.4595279e-02,2.1224443e-02 +1,,density,scatter,2.6505962e-03,2.3698970e-03 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,-1.2721728e-01,1.8810986e-01 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,0.0000000e+00,0.0000000e+00 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,4.1950592e-02,7.3680059e-02 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,scatter,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,scatter,-1.9798576e-02,1.9798576e-02 +1,O16,nuclide_density,nu-fission,9.3632921e-01,2.0322286e+00 +1,O16,nuclide_density,scatter,-2.2042881e-01,1.6781767e-01 +1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,scatter,2.3761276e+00,2.3761276e+00 +1,U235,nuclide_density,nu-fission,7.7441690e+02,8.6460933e+01 +1,U235,nuclide_density,scatter,9.6917103e+01,1.0750582e+01 +1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00 +1,,temperature,scatter,-2.0607998e-06,2.0607998e-06 +1,,temperature,nu-fission,3.8845881e-06,3.5405083e-05 +1,,temperature,scatter,-6.6772431e-07,2.0750826e-07 +1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00 +1,,temperature,scatter,0.0000000e+00,0.0000000e+00 +1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00 +1,,temperature,scatter,0.0000000e+00,0.0000000e+00 +3,,density,flux,-7.6179774e+00,5.0326447e+00 +3,,density,flux,-1.6242019e+01,6.6231774e+00 +1,,density,flux,-2.2394746e-01,5.4534297e-02 +1,,density,flux,-5.7132854e-02,1.6898134e-01 +1,O16,nuclide_density,flux,-1.3403658e+01,1.3017127e+01 +1,O16,nuclide_density,flux,-1.0499715e+01,2.1684953e+01 +1,U235,nuclide_density,flux,-2.4741168e+03,5.7986812e+01 +1,U235,nuclide_density,flux,-1.9955533e+03,4.3254820e+02 +1,,temperature,flux,1.0601815e-04,3.3076550e-04 +1,,temperature,flux,1.6664450e-04,2.8761112e-04 +3,,density,total,-3.3161585e+00,2.5615932e+00 +3,,density,absorption,-4.3491248e-01,5.1559021e-01 +3,,density,scatter,-2.8812460e+00,2.0540305e+00 +3,,density,fission,-2.3060366e-01,3.1191109e-01 +3,,density,nu-fission,-5.6817321e-01,7.6143389e-01 +3,,density,total,-2.8908007e-01,3.9383184e-01 +3,,density,absorption,-2.5237485e-01,3.6565053e-01 +3,,density,scatter,-3.6705221e-02,2.9275371e-02 +3,,density,fission,-2.1271120e-01,3.0873698e-01 +3,,density,nu-fission,-5.1866339e-01,7.5235533e-01 +3,,density,total,2.7320569e+00,8.7709465e+00 +3,,density,absorption,8.4322466e-02,1.2807898e-01 +3,,density,scatter,2.6477344e+00,8.6432567e+00 +3,,density,fission,0.0000000e+00,0.0000000e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,total,0.0000000e+00,0.0000000e+00 +3,,density,absorption,0.0000000e+00,0.0000000e+00 +3,,density,scatter,0.0000000e+00,0.0000000e+00 +3,,density,fission,0.0000000e+00,0.0000000e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,total,4.7523439e-01,2.5652004e-02 +1,,density,absorption,3.1301234e-02,7.6793142e-03 +1,,density,scatter,4.4393315e-01,1.8023252e-02 +1,,density,fission,1.5559355e-02,3.2310294e-03 +1,,density,nu-fission,3.8658109e-02,7.8012233e-03 +1,,density,total,2.2062939e-02,4.3689061e-03 +1,,density,absorption,1.6845626e-02,4.0750601e-03 +1,,density,scatter,5.2173132e-03,2.9663452e-04 +1,,density,fission,1.3503649e-02,3.2642875e-03 +1,,density,nu-fission,3.2945474e-02,7.9507685e-03 +1,,density,total,-9.4735691e-02,2.5584665e-01 +1,,density,absorption,-4.0246397e-03,5.2111770e-03 +1,,density,scatter,-9.0711052e-02,2.5073762e-01 +1,,density,fission,0.0000000e+00,0.0000000e+00 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,density,total,0.0000000e+00,0.0000000e+00 +1,,density,absorption,0.0000000e+00,0.0000000e+00 +1,,density,scatter,0.0000000e+00,0.0000000e+00 +1,,density,fission,0.0000000e+00,0.0000000e+00 +1,,density,nu-fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,total,4.4488215e+01,5.3981603e+00 +1,O16,nuclide_density,absorption,-2.9372611e-01,9.3282813e-01 +1,O16,nuclide_density,scatter,4.4781941e+01,4.5257555e+00 +1,O16,nuclide_density,fission,9.9827872e-02,4.0698927e-01 +1,O16,nuclide_density,nu-fission,2.3342880e-01,9.8749212e-01 +1,O16,nuclide_density,total,4.0914051e-02,6.2301812e-01 +1,O16,nuclide_density,absorption,8.8289354e-02,5.5570922e-01 +1,O16,nuclide_density,scatter,-4.7375303e-02,6.7714071e-02 +1,O16,nuclide_density,fission,1.3678733e-01,4.4006413e-01 +1,O16,nuclide_density,nu-fission,3.3257520e-01,1.0719313e+00 +1,O16,nuclide_density,total,-1.6055918e+01,2.3439335e+01 +1,O16,nuclide_density,absorption,-3.9601231e-01,3.0131622e-01 +1,O16,nuclide_density,scatter,-1.5659906e+01,2.3140407e+01 +1,O16,nuclide_density,fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,total,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,absorption,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,fission,0.0000000e+00,0.0000000e+00 +1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,total,-6.1444902e+02,6.0111692e+01 +1,U235,nuclide_density,absorption,2.3977134e+02,4.5868405e+01 +1,U235,nuclide_density,scatter,-8.5422036e+02,2.4877920e+01 +1,U235,nuclide_density,fission,3.1319132e+02,3.1465126e+01 +1,U235,nuclide_density,nu-fission,7.6415465e+02,7.6501604e+01 +1,U235,nuclide_density,total,5.0916179e+02,3.5499319e+01 +1,U235,nuclide_density,absorption,3.9477360e+02,3.4326165e+01 +1,U235,nuclide_density,scatter,1.1438818e+02,2.7378217e+00 +1,U235,nuclide_density,fission,3.1360739e+02,3.2243605e+01 +1,U235,nuclide_density,nu-fission,7.6527228e+02,7.8683697e+01 +1,U235,nuclide_density,total,-4.1815711e+03,8.7796611e+02 +1,U235,nuclide_density,absorption,-1.0563226e+02,2.5308627e+01 +1,U235,nuclide_density,scatter,-4.0759388e+03,8.5332185e+02 +1,U235,nuclide_density,fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,total,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,absorption,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,fission,0.0000000e+00,0.0000000e+00 +1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00 +1,,temperature,total,2.1367130e-04,1.8632816e-04 +1,,temperature,absorption,6.9456249e-05,3.2199390e-05 +1,,temperature,scatter,1.4421505e-04,1.5754180e-04 +1,,temperature,fission,3.0674981e-06,1.8048377e-05 +1,,temperature,nu-fission,7.4768877e-06,4.3976766e-05 +1,,temperature,total,5.8223116e-06,2.3786474e-05 +1,,temperature,absorption,5.2142698e-06,2.1901911e-05 +1,,temperature,scatter,6.0804185e-07,1.9834890e-06 +1,,temperature,fission,3.0674200e-06,1.8048542e-05 +1,,temperature,nu-fission,7.4767028e-06,4.3977153e-05 +1,,temperature,total,2.1510156e-04,4.6388707e-04 +1,,temperature,absorption,2.2409527e-06,8.2901060e-06 +1,,temperature,scatter,2.1286061e-04,4.5560034e-04 +1,,temperature,fission,0.0000000e+00,0.0000000e+00 +1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00 +1,,temperature,total,0.0000000e+00,0.0000000e+00 +1,,temperature,absorption,0.0000000e+00,0.0000000e+00 +1,,temperature,scatter,0.0000000e+00,0.0000000e+00 +1,,temperature,fission,0.0000000e+00,0.0000000e+00 +1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,absorption,-1.6517201e-01,3.0984738e-01 +3,,density,absorption,8.0402344e-03,1.4689335e-01 +1,,density,absorption,2.9069882e-02,2.2139609e-03 +1,,density,absorption,-9.4690065e-03,8.6605392e-03 +1,O16,nuclide_density,absorption,7.6911962e-01,4.1945687e-01 +1,O16,nuclide_density,absorption,-6.3724795e-01,6.6341488e-01 +1,U235,nuclide_density,absorption,1.4109543e+02,1.8518890e+01 +1,U235,nuclide_density,absorption,-1.2052822e+02,3.2084002e+01 +1,,temperature,absorption,3.9995404e-05,2.1703384e-05 +1,,temperature,absorption,2.7274361e-06,8.9339257e-06 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,-8.7934551e-01,1.0210652e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,-3.5236516e-03,3.5236516e-03 +3,,density,nu-fission,7.2953012e-02,2.9725863e-01 +3,,density,scatter,-2.2716410e+00,1.2546347e+00 +3,,density,nu-fission,8.9171481e-02,3.0634856e-01 +3,,density,scatter,-1.0151747e-03,9.6296775e-03 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,2.3163923e+00,4.8766690e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,0.0000000e+00,0.0000000e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,4.0762434e-01,3.8504141e+00 +3,,density,nu-fission,0.0000000e+00,0.0000000e+00 +3,,density,scatter,0.0000000e+00,0.0000000e+00 diff --git a/tests/test_diff_tally/test_diff_tally.py b/tests/test_diff_tally/test_diff_tally.py new file mode 100644 index 000000000..36287a66a --- /dev/null +++ b/tests/test_diff_tally/test_diff_tally.py @@ -0,0 +1,140 @@ +#!/usr/bin/env python + +import glob +import os +import sys + +import pandas as pd + +sys.path.insert(0, os.pardir) +from testing_harness import PyAPITestHarness +import openmc + +class DiffTallyTestHarness(PyAPITestHarness): + def _build_inputs(self): + # Build default materials/geometry + self._input_set.build_default_materials_and_geometry() + + # Set settings explicitly + self._input_set.settings.batches = 3 + self._input_set.settings.inactive = 0 + self._input_set.settings.particles = 100 + self._input_set.settings.source = openmc.Source(space=openmc.stats.Box( + [-160, -160, -183], [160, 160, 183])) + self._input_set.settings.temperature['multipole'] = True + + self._input_set.tallies = openmc.Tallies() + + filt_mats = openmc.MaterialFilter((1, 3)) + filt_eout = openmc.EnergyoutFilter((0.0, 0.625, 20.0e6)) + + # We want density derivatives for both water and fuel to get coverage + # for both fissile and non-fissile materials. + d1 = openmc.TallyDerivative(derivative_id=1) + d1.variable = 'density' + d1.material = 3 + d2 = openmc.TallyDerivative(derivative_id=2) + d2.variable = 'density' + d2.material = 1 + + # O-16 is a good nuclide to test against because it is present in both + # water and fuel. Some routines need to recognize that they have the + # perturbed nuclide but not the perturbed material. + d3 = openmc.TallyDerivative(derivative_id=3) + d3.variable = 'nuclide_density' + d3.material = 1 + d3.nuclide = 'O16' + + # A fissile nuclide, just for good measure. + d4 = openmc.TallyDerivative(derivative_id=4) + d4.variable = 'nuclide_density' + d4.material = 1 + d4.nuclide = 'U235' + + # Temperature derivatives. + d5 = openmc.TallyDerivative(derivative_id=5) + d5.variable = 'temperature' + d5.material = 1 + + derivs = [d1, d2, d3, d4, d5] + + # Cover the flux score. + for i in range(5): + t = openmc.Tally() + t.add_score('flux') + t.add_filter(filt_mats) + t.derivative = derivs[i] + self._input_set.tallies.append(t) + + # Cover supported scores with a collision estimator. + for i in range(5): + t = openmc.Tally() + t.add_score('total') + t.add_score('absorption') + t.add_score('scatter') + t.add_score('fission') + t.add_score('nu-fission') + t.add_filter(filt_mats) + t.add_nuclide('total') + t.add_nuclide('U235') + t.derivative = derivs[i] + self._input_set.tallies.append(t) + + # Cover an analog estimator. + for i in range(5): + t = openmc.Tally() + t.add_score('absorption') + t.add_filter(filt_mats) + t.estimator = 'analog' + t.derivative = derivs[i] + self._input_set.tallies.append(t) + + # Energyout filter and total nuclide for the density derivatives. + for i in range(2): + t = openmc.Tally() + t.add_score('nu-fission') + t.add_score('scatter') + t.add_filter(filt_mats) + t.add_filter(filt_eout) + t.add_nuclide('total') + t.add_nuclide('U235') + t.derivative = derivs[i] + self._input_set.tallies.append(t) + + # Energyout filter without total nuclide for other derivatives. + for i in range(2, 5): + t = openmc.Tally() + t.add_score('nu-fission') + t.add_score('scatter') + t.add_filter(filt_mats) + t.add_filter(filt_eout) + t.add_nuclide('U235') + t.derivative = derivs[i] + self._input_set.tallies.append(t) + + self._input_set.export() + + def _get_results(self): + # Read the statepoint and summary files. + statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0] + sp = openmc.StatePoint(statepoint) + + # Extract the tally data as a Pandas DataFrame. + df = pd.DataFrame() + for t in sp.tallies.values(): + df = df.append(t.get_pandas_dataframe(), ignore_index=True) + + # Extract the relevant data as a CSV string. + cols = ('d_material', 'd_nuclide', 'd_variable', 'score', 'mean', + 'std. dev.') + return df.to_csv(None, columns=cols, index=False, float_format='%.7e') + + def _cleanup(self): + super(DiffTallyTestHarness, self)._cleanup() + f = os.path.join(os.getcwd(), 'tallies.xml') + if os.path.exists(f): os.remove(f) + + +if __name__ == '__main__': + harness = DiffTallyTestHarness('statepoint.3.h5', True) + harness.main() diff --git a/tests/test_mgxs_library_condense/inputs_true.dat b/tests/test_mgxs_library_condense/inputs_true.dat index e66d3d111..6b3f5ab3a 100644 --- a/tests/test_mgxs_library_condense/inputs_true.dat +++ b/tests/test_mgxs_library_condense/inputs_true.dat @@ -1 +1 @@ -e86f24e20f37096c7898f459fc5bf336f3e0670fb30f321443f2bb3a01a154ef3d234bb48d4cee8e0d3730fd182b6a2051ec1b4c09d771420886a56ba928fdd9 \ No newline at end of file +4291b40470e7d59383c9e51c4178ca923b698cb1aaea16c1982fe3789ca980df10a65b84fb5021dacd4d290ebc232c2579f99b6c990fb6b8f28f67eef2aabcb2 \ No newline at end of file diff --git a/tests/test_mgxs_library_condense/results_true.dat b/tests/test_mgxs_library_condense/results_true.dat index 525d76efe..6f29fbc4c 100644 --- a/tests/test_mgxs_library_condense/results_true.dat +++ b/tests/test_mgxs_library_condense/results_true.dat @@ -40,6 +40,8 @@ 0 10000 1 total 4.996730e-07 3.650635e-08 material group in nuclide mean std. dev. 0 10000 1 total 0.090004 0.006367 + material group in group out nuclide mean std. dev. +0 10000 1 1 total 0.084542 0.005716 material delayedgroup group in nuclide mean std. dev. 0 10000 1 1 total 0.000021 0.000001 1 10000 2 1 total 0.000110 0.000008 @@ -68,6 +70,13 @@ 3 10000 4 1 total 0.302780 0.109110 4 10000 5 1 total 0.000000 0.000000 5 10000 6 1 total 0.000000 0.000000 + material delayedgroup group in group out nuclide mean std. dev. +0 10000 1 1 1 total 0.000000 0.000000 +1 10000 2 1 1 total 0.000384 0.000236 +2 10000 3 1 1 total 0.000179 0.000180 +3 10000 4 1 1 total 0.000730 0.000188 +4 10000 5 1 1 total 0.000000 0.000000 +5 10000 6 1 1 total 0.000000 0.000000 material group in nuclide mean std. dev. 0 10001 1 total 0.311594 0.013793 material group in nuclide mean std. dev. @@ -110,6 +119,8 @@ 0 10001 1 total 5.454760e-07 4.949800e-08 material group in nuclide mean std. dev. 0 10001 1 total 0.0 0.0 + material group in group out nuclide mean std. dev. +0 10001 1 1 total 0.0 0.0 material delayedgroup group in nuclide mean std. dev. 0 10001 1 1 total 0.0 0.0 1 10001 2 1 total 0.0 0.0 @@ -138,6 +149,13 @@ 3 10001 4 1 total 0.0 0.0 4 10001 5 1 total 0.0 0.0 5 10001 6 1 total 0.0 0.0 + material delayedgroup group in group out nuclide mean std. dev. +0 10001 1 1 1 total 0.0 0.0 +1 10001 2 1 1 total 0.0 0.0 +2 10001 3 1 1 total 0.0 0.0 +3 10001 4 1 1 total 0.0 0.0 +4 10001 5 1 1 total 0.0 0.0 +5 10001 6 1 1 total 0.0 0.0 material group in nuclide mean std. dev. 0 10002 1 total 0.904999 0.043964 material group in nuclide mean std. dev. @@ -180,6 +198,8 @@ 0 10002 1 total 5.773006e-07 5.322132e-08 material group in nuclide mean std. dev. 0 10002 1 total 0.0 0.0 + material group in group out nuclide mean std. dev. +0 10002 1 1 total 0.0 0.0 material delayedgroup group in nuclide mean std. dev. 0 10002 1 1 total 0.0 0.0 1 10002 2 1 total 0.0 0.0 @@ -208,3 +228,10 @@ 3 10002 4 1 total 0.0 0.0 4 10002 5 1 total 0.0 0.0 5 10002 6 1 total 0.0 0.0 + material delayedgroup group in group out nuclide mean std. dev. +0 10002 1 1 1 total 0.0 0.0 +1 10002 2 1 1 total 0.0 0.0 +2 10002 3 1 1 total 0.0 0.0 +3 10002 4 1 1 total 0.0 0.0 +4 10002 5 1 1 total 0.0 0.0 +5 10002 6 1 1 total 0.0 0.0 diff --git a/tests/test_mgxs_library_distribcell/inputs_true.dat b/tests/test_mgxs_library_distribcell/inputs_true.dat index f03c41fe7..f00c2133d 100644 --- a/tests/test_mgxs_library_distribcell/inputs_true.dat +++ b/tests/test_mgxs_library_distribcell/inputs_true.dat @@ -1 +1 @@ -67ee414eed54f596464831e734ac365e43b88bf03297e2b08adfb97510e1d8e631ff45061d5b99aa0ad85ec08d5b51e341866e60f35cdbfae078030bcde6c794 \ No newline at end of file +df187239f7481867cc09138709da90bc28eeb647b4bcbb08b894e7fdf6ff45510341b77e34754c2358f0e0665f8e552bb1eafed8c42cd2c50595b60366320ce4 \ No newline at end of file diff --git a/tests/test_mgxs_library_distribcell/results_true.dat b/tests/test_mgxs_library_distribcell/results_true.dat index ef4d93eb9..2b43005aa 100644 --- a/tests/test_mgxs_library_distribcell/results_true.dat +++ b/tests/test_mgxs_library_distribcell/results_true.dat @@ -40,6 +40,8 @@ 0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 5.139437e-07 2.133314e-08 avg(distribcell) group in nuclide mean std. dev. 0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.091725 0.003604 + avg(distribcell) group in group out nuclide mean std. dev. +0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.093985 0.005872 avg(distribcell) delayedgroup group in nuclide mean std. dev. 0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.000021 8.253907e-07 1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 total 0.000112 4.284000e-06 @@ -68,3 +70,10 @@ 3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 total 0.000000 0.000000 4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 total 0.000000 0.000000 5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 total 2.853000 4.034751 + avg(distribcell) delayedgroup group in group out nuclide mean std. dev. +0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 1 total 0.000000 0.000000 +1 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 2 1 1 total 0.000175 0.000175 +2 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 3 1 1 total 0.000178 0.000178 +3 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 4 1 1 total 0.000000 0.000000 +4 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 5 1 1 total 0.000000 0.000000 +5 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 6 1 1 total 0.000178 0.000178 diff --git a/tests/test_mgxs_library_hdf5/inputs_true.dat b/tests/test_mgxs_library_hdf5/inputs_true.dat index e66d3d111..6b3f5ab3a 100644 --- a/tests/test_mgxs_library_hdf5/inputs_true.dat +++ b/tests/test_mgxs_library_hdf5/inputs_true.dat @@ -1 +1 @@ -e86f24e20f37096c7898f459fc5bf336f3e0670fb30f321443f2bb3a01a154ef3d234bb48d4cee8e0d3730fd182b6a2051ec1b4c09d771420886a56ba928fdd9 \ No newline at end of file +4291b40470e7d59383c9e51c4178ca923b698cb1aaea16c1982fe3789ca980df10a65b84fb5021dacd4d290ebc232c2579f99b6c990fb6b8f28f67eef2aabcb2 \ No newline at end of file diff --git a/tests/test_mgxs_library_hdf5/results_true.dat b/tests/test_mgxs_library_hdf5/results_true.dat index 8b02203c6..d491e72ba 100644 --- a/tests/test_mgxs_library_hdf5/results_true.dat +++ b/tests/test_mgxs_library_hdf5/results_true.dat @@ -72,6 +72,11 @@ domain=10000 type=inverse-velocity domain=10000 type=prompt-nu-fission [1.92392215e-02 4.66719027e-01] [1.30950595e-03 4.14108704e-02] +domain=10000 type=prompt-nu-fission matrix +[[2.01424282e-02 0.00000000e+00] + [4.45819177e-01 0.00000000e+00]] +[[3.14909168e-03 0.00000000e+00] + [2.86750787e-02 0.00000000e+00]] domain=10000 type=delayed-nu-fission [[2.29808234e-05 1.06974158e-04] [1.43606337e-04 5.52167907e-04] @@ -124,6 +129,41 @@ domain=10000 type=decay-rate [0.00000000e+00 1.09109511e-01] [0.00000000e+00 0.00000000e+00] [0.00000000e+00 0.00000000e+00]] +domain=10000 type=delayed-nu-fission matrix +[[[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [2.53814542e-03 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [1.18579166e-03 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [8.59787018e-04 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]]] +[[[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [1.56094584e-03 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [1.18610401e-03 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [2.22194634e-04 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]]] domain=10001 type=total [3.13737671e-01 3.00821402e-01] [1.55819024e-02 2.80524484e-02] @@ -198,6 +238,11 @@ domain=10001 type=inverse-velocity domain=10001 type=prompt-nu-fission [0.00000000e+00 0.00000000e+00] [0.00000000e+00 0.00000000e+00] +domain=10001 type=prompt-nu-fission matrix +[[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] +[[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] domain=10001 type=delayed-nu-fission [[0.00000000e+00 0.00000000e+00] [0.00000000e+00 0.00000000e+00] @@ -250,6 +295,41 @@ domain=10001 type=decay-rate [0.00000000e+00 0.00000000e+00] [0.00000000e+00 0.00000000e+00] [0.00000000e+00 0.00000000e+00]] +domain=10001 type=delayed-nu-fission matrix +[[[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]]] +[[[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]]] domain=10002 type=total [6.64572261e-01 2.05238401e+00] [3.12147519e-02 2.24342907e-01] @@ -324,6 +404,11 @@ domain=10002 type=inverse-velocity domain=10002 type=prompt-nu-fission [0.00000000e+00 0.00000000e+00] [0.00000000e+00 0.00000000e+00] +domain=10002 type=prompt-nu-fission matrix +[[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] +[[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] domain=10002 type=delayed-nu-fission [[0.00000000e+00 0.00000000e+00] [0.00000000e+00 0.00000000e+00] @@ -376,3 +461,38 @@ domain=10002 type=decay-rate [0.00000000e+00 0.00000000e+00] [0.00000000e+00 0.00000000e+00] [0.00000000e+00 0.00000000e+00]] +domain=10002 type=delayed-nu-fission matrix +[[[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]]] +[[[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]] + + [[0.00000000e+00 0.00000000e+00] + [0.00000000e+00 0.00000000e+00]]] diff --git a/tests/test_mgxs_library_mesh/inputs_true.dat b/tests/test_mgxs_library_mesh/inputs_true.dat index 3b297354d..70799f147 100644 --- a/tests/test_mgxs_library_mesh/inputs_true.dat +++ b/tests/test_mgxs_library_mesh/inputs_true.dat @@ -1 +1 @@ -a0d62fc011ae33756cd87202432f5c49f847793a40ec50efc6e4366755a70196f9ae970b24244766c011e2991961ea92317b53f6839b89363262b4e67ed4b289 \ No newline at end of file +03d894a7995ac40f7971b17349b460f4b563cb1c94abaa7e7241e6f7edad7f0cb80d3f80829db14ea6b7d70d18197cb15e308047cd10f699d834178eeb6396be \ No newline at end of file diff --git a/tests/test_mgxs_library_mesh/results_true.dat b/tests/test_mgxs_library_mesh/results_true.dat index 2fc403c89..b4681d25c 100644 --- a/tests/test_mgxs_library_mesh/results_true.dat +++ b/tests/test_mgxs_library_mesh/results_true.dat @@ -130,6 +130,12 @@ 1 1 2 1 1 total 0.020397 0.008086 2 2 1 1 1 total 0.025824 0.003192 3 2 2 1 1 total 0.020865 0.004879 + mesh 1 group in group out nuclide mean std. dev. + x y z +0 1 1 1 1 1 total 0.020874 0.002977 +1 1 2 1 1 1 total 0.017348 0.008786 +2 2 1 1 1 1 total 0.020409 0.003354 +3 2 2 1 1 1 total 0.011105 0.003806 mesh 1 delayedgroup group in nuclide mean std. dev. x y z 0 1 1 1 1 1 total 0.000005 1.004627e-06 @@ -234,3 +240,29 @@ 21 2 2 1 4 1 total 0.00000 0.00000 22 2 2 1 5 1 total 0.00000 0.00000 23 2 2 1 6 1 total 0.00000 0.00000 + mesh 1 delayedgroup group in group out nuclide mean std. dev. + x y z +0 1 1 1 1 1 1 total 0.000000 0.000000 +1 1 1 1 2 1 1 total 0.000000 0.000000 +2 1 1 1 3 1 1 total 0.000000 0.000000 +3 1 1 1 4 1 1 total 0.000000 0.000000 +4 1 1 1 5 1 1 total 0.000185 0.000186 +5 1 1 1 6 1 1 total 0.000000 0.000000 +6 1 2 1 1 1 1 total 0.000000 0.000000 +7 1 2 1 2 1 1 total 0.000000 0.000000 +8 1 2 1 3 1 1 total 0.000000 0.000000 +9 1 2 1 4 1 1 total 0.000000 0.000000 +10 1 2 1 5 1 1 total 0.000000 0.000000 +11 1 2 1 6 1 1 total 0.000000 0.000000 +12 2 1 1 1 1 1 total 0.000000 0.000000 +13 2 1 1 2 1 1 total 0.000000 0.000000 +14 2 1 1 3 1 1 total 0.000000 0.000000 +15 2 1 1 4 1 1 total 0.000000 0.000000 +16 2 1 1 5 1 1 total 0.000000 0.000000 +17 2 1 1 6 1 1 total 0.000000 0.000000 +18 2 2 1 1 1 1 total 0.000000 0.000000 +19 2 2 1 2 1 1 total 0.000000 0.000000 +20 2 2 1 3 1 1 total 0.000000 0.000000 +21 2 2 1 4 1 1 total 0.000000 0.000000 +22 2 2 1 5 1 1 total 0.000000 0.000000 +23 2 2 1 6 1 1 total 0.000000 0.000000 diff --git a/tests/test_mgxs_library_no_nuclides/inputs_true.dat b/tests/test_mgxs_library_no_nuclides/inputs_true.dat index e66d3d111..6b3f5ab3a 100644 --- a/tests/test_mgxs_library_no_nuclides/inputs_true.dat +++ b/tests/test_mgxs_library_no_nuclides/inputs_true.dat @@ -1 +1 @@ -e86f24e20f37096c7898f459fc5bf336f3e0670fb30f321443f2bb3a01a154ef3d234bb48d4cee8e0d3730fd182b6a2051ec1b4c09d771420886a56ba928fdd9 \ No newline at end of file +4291b40470e7d59383c9e51c4178ca923b698cb1aaea16c1982fe3789ca980df10a65b84fb5021dacd4d290ebc232c2579f99b6c990fb6b8f28f67eef2aabcb2 \ No newline at end of file diff --git a/tests/test_mgxs_library_no_nuclides/results_true.dat b/tests/test_mgxs_library_no_nuclides/results_true.dat index 07813519a..fac47af46 100644 --- a/tests/test_mgxs_library_no_nuclides/results_true.dat +++ b/tests/test_mgxs_library_no_nuclides/results_true.dat @@ -84,6 +84,11 @@ material group in nuclide mean std. dev. 1 10000 1 total 0.019239 0.001310 0 10000 2 total 0.466719 0.041411 + material group in group out nuclide mean std. dev. +3 10000 1 1 total 0.020142 0.003149 +2 10000 1 2 total 0.000000 0.000000 +1 10000 2 1 total 0.445819 0.028675 +0 10000 2 2 total 0.000000 0.000000 material delayedgroup group in nuclide mean std. dev. 1 10000 1 1 total 0.000023 0.000002 3 10000 2 1 total 0.000144 0.000011 @@ -136,6 +141,31 @@ 6 10000 4 2 total 0.302780 0.109110 8 10000 5 2 total 0.000000 0.000000 10 10000 6 2 total 0.000000 0.000000 + material delayedgroup group in group out nuclide mean std. dev. +3 10000 1 1 1 total 0.000000 0.000000 +7 10000 2 1 1 total 0.000000 0.000000 +11 10000 3 1 1 total 0.000000 0.000000 +15 10000 4 1 1 total 0.000000 0.000000 +19 10000 5 1 1 total 0.000000 0.000000 +23 10000 6 1 1 total 0.000000 0.000000 +2 10000 1 1 2 total 0.000000 0.000000 +6 10000 2 1 2 total 0.000000 0.000000 +10 10000 3 1 2 total 0.000000 0.000000 +14 10000 4 1 2 total 0.000000 0.000000 +18 10000 5 1 2 total 0.000000 0.000000 +22 10000 6 1 2 total 0.000000 0.000000 +1 10000 1 2 1 total 0.000000 0.000000 +5 10000 2 2 1 total 0.002538 0.001561 +9 10000 3 2 1 total 0.001186 0.001186 +13 10000 4 2 1 total 0.000860 0.000222 +17 10000 5 2 1 total 0.000000 0.000000 +21 10000 6 2 1 total 0.000000 0.000000 +0 10000 1 2 2 total 0.000000 0.000000 +4 10000 2 2 2 total 0.000000 0.000000 +8 10000 3 2 2 total 0.000000 0.000000 +12 10000 4 2 2 total 0.000000 0.000000 +16 10000 5 2 2 total 0.000000 0.000000 +20 10000 6 2 2 total 0.000000 0.000000 material group in nuclide mean std. dev. 1 10001 1 total 0.313738 0.015582 0 10001 2 total 0.300821 0.028052 @@ -222,6 +252,11 @@ material group in nuclide mean std. dev. 1 10001 1 total 0.0 0.0 0 10001 2 total 0.0 0.0 + material group in group out nuclide mean std. dev. +3 10001 1 1 total 0.0 0.0 +2 10001 1 2 total 0.0 0.0 +1 10001 2 1 total 0.0 0.0 +0 10001 2 2 total 0.0 0.0 material delayedgroup group in nuclide mean std. dev. 1 10001 1 1 total 0.0 0.0 3 10001 2 1 total 0.0 0.0 @@ -274,6 +309,31 @@ 6 10001 4 2 total 0.0 0.0 8 10001 5 2 total 0.0 0.0 10 10001 6 2 total 0.0 0.0 + material delayedgroup group in group out nuclide mean std. dev. +3 10001 1 1 1 total 0.0 0.0 +7 10001 2 1 1 total 0.0 0.0 +11 10001 3 1 1 total 0.0 0.0 +15 10001 4 1 1 total 0.0 0.0 +19 10001 5 1 1 total 0.0 0.0 +23 10001 6 1 1 total 0.0 0.0 +2 10001 1 1 2 total 0.0 0.0 +6 10001 2 1 2 total 0.0 0.0 +10 10001 3 1 2 total 0.0 0.0 +14 10001 4 1 2 total 0.0 0.0 +18 10001 5 1 2 total 0.0 0.0 +22 10001 6 1 2 total 0.0 0.0 +1 10001 1 2 1 total 0.0 0.0 +5 10001 2 2 1 total 0.0 0.0 +9 10001 3 2 1 total 0.0 0.0 +13 10001 4 2 1 total 0.0 0.0 +17 10001 5 2 1 total 0.0 0.0 +21 10001 6 2 1 total 0.0 0.0 +0 10001 1 2 2 total 0.0 0.0 +4 10001 2 2 2 total 0.0 0.0 +8 10001 3 2 2 total 0.0 0.0 +12 10001 4 2 2 total 0.0 0.0 +16 10001 5 2 2 total 0.0 0.0 +20 10001 6 2 2 total 0.0 0.0 material group in nuclide mean std. dev. 1 10002 1 total 0.664572 0.031215 0 10002 2 total 2.052384 0.224343 @@ -360,6 +420,11 @@ material group in nuclide mean std. dev. 1 10002 1 total 0.0 0.0 0 10002 2 total 0.0 0.0 + material group in group out nuclide mean std. dev. +3 10002 1 1 total 0.0 0.0 +2 10002 1 2 total 0.0 0.0 +1 10002 2 1 total 0.0 0.0 +0 10002 2 2 total 0.0 0.0 material delayedgroup group in nuclide mean std. dev. 1 10002 1 1 total 0.0 0.0 3 10002 2 1 total 0.0 0.0 @@ -412,3 +477,28 @@ 6 10002 4 2 total 0.0 0.0 8 10002 5 2 total 0.0 0.0 10 10002 6 2 total 0.0 0.0 + material delayedgroup group in group out nuclide mean std. dev. +3 10002 1 1 1 total 0.0 0.0 +7 10002 2 1 1 total 0.0 0.0 +11 10002 3 1 1 total 0.0 0.0 +15 10002 4 1 1 total 0.0 0.0 +19 10002 5 1 1 total 0.0 0.0 +23 10002 6 1 1 total 0.0 0.0 +2 10002 1 1 2 total 0.0 0.0 +6 10002 2 1 2 total 0.0 0.0 +10 10002 3 1 2 total 0.0 0.0 +14 10002 4 1 2 total 0.0 0.0 +18 10002 5 1 2 total 0.0 0.0 +22 10002 6 1 2 total 0.0 0.0 +1 10002 1 2 1 total 0.0 0.0 +5 10002 2 2 1 total 0.0 0.0 +9 10002 3 2 1 total 0.0 0.0 +13 10002 4 2 1 total 0.0 0.0 +17 10002 5 2 1 total 0.0 0.0 +21 10002 6 2 1 total 0.0 0.0 +0 10002 1 2 2 total 0.0 0.0 +4 10002 2 2 2 total 0.0 0.0 +8 10002 3 2 2 total 0.0 0.0 +12 10002 4 2 2 total 0.0 0.0 +16 10002 5 2 2 total 0.0 0.0 +20 10002 6 2 2 total 0.0 0.0 diff --git a/tests/test_mgxs_library_nuclides/inputs_true.dat b/tests/test_mgxs_library_nuclides/inputs_true.dat index d7c6eeab0..ed6f4d443 100644 --- a/tests/test_mgxs_library_nuclides/inputs_true.dat +++ b/tests/test_mgxs_library_nuclides/inputs_true.dat @@ -1 +1 @@ -1e8d0f408ba9d47fc9278e750dbf68d3104d248f7783ca1185c8b7666c9ee2ba3ec6fcf4a49bb614956bb8eb3a57dc1c104725b4160171dfbe6dc972a268df56 \ No newline at end of file +eed0190893105747f0146472dfea86ca58b9fc2e8d039961d9fef8f235f79da32215f9b59b840419e6b69c5a04544b4659f2634b63e3a76b6b0f8c75d05e416c \ No newline at end of file diff --git a/tests/test_mgxs_library_nuclides/results_true.dat b/tests/test_mgxs_library_nuclides/results_true.dat index b55ab34ab..bdc0d1768 100644 --- a/tests/test_mgxs_library_nuclides/results_true.dat +++ b/tests/test_mgxs_library_nuclides/results_true.dat @@ -1 +1 @@ -ce682f577fd65dd9b6e5da58763cf2ae4cf8e5041b38b5d1ed64d3ba7a3a4df7338aadfa8830a0764fc9ffb737f05fc989a494158c714047186dd18605a7d2b8 \ No newline at end of file +8bc6694ee99cc05ec59143dac065f3f51026713044c48a36ac486f46014289c1f18d4d8e14ef0bb468050aea0bac5628425fd76154a01c7438c7d85139aa2f06 \ No newline at end of file