diff --git a/docs/source/io_formats/mgxs_library.rst b/docs/source/io_formats/mgxs_library.rst index fd00ad8c34..57a3770402 100644 --- a/docs/source/io_formats/mgxs_library.rst +++ b/docs/source/io_formats/mgxs_library.rst @@ -15,8 +15,9 @@ to be used in the ``mgxs.h5`` file. The multi-group library is provided in the HDF5_ format. This library must provide some meta-data about the library itself (such as the number of -groups and the group structure, etc.) as well as the actual cross section -data itself for each of the necessary nuclides or materials. +energy groups, delayed groups, and the energy group structure, etc.) as +well as the actual cross section data itself for each of the necessary +nuclides or materials. .. _HDF5: http://www.hdfgroup.org/HDF5/ @@ -28,7 +29,8 @@ MGXS Library Specification **/** -:Attributes: - **groups** (*int*) -- Number of energy groups +:Attributes: - **energy_groups** (*int*) -- Number of energy groups + - **delayed_groups** (*int*) -- Number of delayed groups (optional) - **group structure** (*double[]*) -- Monotonically increasing list of group boundaries, in units of eV. The length of this array should be the number of groups plus 1. @@ -119,8 +121,8 @@ Temperature-dependent data, provided for temperature K. the `nu-fission` data must represent the fission neutron energy spectra as well and thus will have one additional dimension for the outgoing energy group. In this case, `nu-fission` has the - same dimensionality as `multiplicity_matrix`. - - **inverse_velocities** (*double[]*) -- Average inverse velocity + same dimensionality as `multiplicity matrix`. + - **inverse-velocity** (*double[]*) -- Average inverse velocity for each of the groups in the library. This dataset is optional. **//K/scatter_data/** diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 9df70ed973..92f53e6bf9 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1629,8 +1629,6 @@ The ```` element accepts the following sub-elements: - .. note:: This filter type is not used in the multi-group :ref:`energy_mode`. - :nuclides: If specified, the scores listed will be for particular nuclides, not the summation of reactions from all nuclides. The format for nuclides should be @@ -1811,12 +1809,10 @@ The ```` element accepts the following sub-elements: |Score | Description | +======================+===================================================+ |delayed-nu-fission |Total production of delayed neutrons due to | - | |fission. This score type is not used in the | - | |multi-group :ref:`energy_mode`. | + | |fission. | +----------------------+---------------------------------------------------+ |prompt-nu-fission |Total production of prompt neutrons due to | - | |fission. This score type is not used in the | - | |multi-group :ref:`energy_mode`. | + | |fission. | +----------------------+---------------------------------------------------+ |nu-fission |Total production of neutrons due to fission. | +----------------------+---------------------------------------------------+ @@ -1845,8 +1841,6 @@ The ```` element accepts the following sub-elements: +----------------------+---------------------------------------------------+ |inverse-velocity |The flux-weighted inverse velocity where the | | |velocity is in units of centimeters per second. | - | |This score type is not used in the | - | |multi-group :ref:`energy_mode`. | +----------------------+---------------------------------------------------+ |kappa-fission |The recoverable energy production rate due to | | |fission. The recoverable energy is defined as the | @@ -1882,8 +1876,6 @@ The ```` element accepts the following sub-elements: +----------------------+---------------------------------------------------+ |decay-rate |The delayed-nu-fission-weighted decay rate where | | |the decay rate is in units of inverse seconds. | - | |This score type is not used in the | - | |multi-group :ref:`energy_mode`. | +----------------------+---------------------------------------------------+ .. note:: diff --git a/examples/xml/pincell_multigroup/geometry.xml b/examples/xml/pincell_multigroup/geometry.xml index 94c70cb340..7f2fb36e82 100644 --- a/examples/xml/pincell_multigroup/geometry.xml +++ b/examples/xml/pincell_multigroup/geometry.xml @@ -1,16 +1,10 @@ + - - - - - - - - - - - - - - + + + + + + + diff --git a/examples/xml/pincell_multigroup/materials.xml b/examples/xml/pincell_multigroup/materials.xml index 61a5cc697a..f0225b1009 100644 --- a/examples/xml/pincell_multigroup/materials.xml +++ b/examples/xml/pincell_multigroup/materials.xml @@ -1,54 +1,12 @@ - + - - ./mgxs.h5 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + ./mgxs.h5 + + + + + + + + diff --git a/examples/xml/pincell_multigroup/mgxs.h5 b/examples/xml/pincell_multigroup/mgxs.h5 index 77f445c6d0..1d5561b000 100644 Binary files a/examples/xml/pincell_multigroup/mgxs.h5 and b/examples/xml/pincell_multigroup/mgxs.h5 differ diff --git a/examples/xml/pincell_multigroup/settings.xml b/examples/xml/pincell_multigroup/settings.xml index 9adc17c9d0..6917e8b1d7 100644 --- a/examples/xml/pincell_multigroup/settings.xml +++ b/examples/xml/pincell_multigroup/settings.xml @@ -1,37 +1,14 @@ - + - - multi-group - - - - 100 - 10 - 1000 - - - - - - - -0.63 -0.63 -1E50 - 0.63 0.63 1E50 - - - - - - true - true - - - false - + + 1000 + 100 + 10 + + + + -0.63 -0.63 -1 0.63 0.63 1 + + + multi-group diff --git a/examples/xml/pincell_multigroup/tallies.xml b/examples/xml/pincell_multigroup/tallies.xml index df65b461db..1591a556b4 100644 --- a/examples/xml/pincell_multigroup/tallies.xml +++ b/examples/xml/pincell_multigroup/tallies.xml @@ -1,4 +1,4 @@ - + 100 100 1 @@ -6,7 +6,7 @@ 0.63 0.63 1e+50 - + flux fission nu-fission diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index a8f54d16c0..f6bf4dab27 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -63,8 +63,8 @@ class Library(object): The highest legendre moment in the scattering matrices (default is 0) energy_groups : openmc.mgxs.EnergyGroups Energy group structure for energy condensation - delayed_groups : list of int - Delayed groups to filter out the xs + num_delayed_groups : int + Number of delayed groups estimator : str or None The tally estimator used to compute multi-group cross sections. If None, the default for each MGXS type is used. @@ -99,7 +99,7 @@ class Library(object): self._domain_type = None self._domains = 'all' self._energy_groups = None - self._delayed_groups = None + self._num_delayed_groups = 0 self._correction = 'P0' self._legendre_order = 0 self._tally_trigger = None @@ -132,7 +132,7 @@ class Library(object): clone._correction = self.correction clone._legendre_order = self.legendre_order clone._energy_groups = copy.deepcopy(self.energy_groups, memo) - clone._delayed_groups = copy.deepcopy(self.delayed_groups, memo) + clone._num_delayed_groups = self.num_delayed_groups clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo) clone._all_mgxs = copy.deepcopy(self.all_mgxs) clone._sp_filename = self._sp_filename @@ -202,8 +202,8 @@ class Library(object): return self._energy_groups @property - def delayed_groups(self): - return self._delayed_groups + def num_delayed_groups(self): + return self._num_delayed_groups @property def correction(self): @@ -225,13 +225,6 @@ class Library(object): def num_groups(self): return self.energy_groups.num_groups - @property - def num_delayed_groups(self): - if self.delayed_groups == None: - return 0 - else: - return len(self.delayed_groups) - @property def all_mgxs(self): return self._all_mgxs @@ -331,22 +324,14 @@ class Library(object): cv.check_type('energy groups', energy_groups, openmc.mgxs.EnergyGroups) self._energy_groups = energy_groups - @delayed_groups.setter - def delayed_groups(self, delayed_groups): + @num_delayed_groups.setter + def num_delayed_groups(self, num_delayed_groups): - if delayed_groups != None: - - cv.check_type('delayed groups', delayed_groups, list, int) - cv.check_greater_than('num delayed groups', len(delayed_groups), 0) - - # Check that the groups are within [1, MAX_DELAYED_GROUPS] - for group in delayed_groups: - cv.check_greater_than('delayed group', group, 0) - cv.check_less_than('delayed group', group, - openmc.mgxs.MAX_DELAYED_GROUPS, - equality=True) - - self._delayed_groups = delayed_groups + cv.check_less_than('num delayed groups', num_delayed_groups, + openmc.mgxs.MAX_DELAYED_GROUPS, equality=True) + cv.check_greater_than('num delayed groups', num_delayed_groups, 0, + equality=True) + self._num_delayed_groups = num_delayed_groups @correction.setter def correction(self, correction): @@ -431,7 +416,12 @@ class Library(object): mgxs.estimator = self.estimator if mgxs_type in openmc.mgxs.MDGXS_TYPES: - mgxs.delayed_groups = self.delayed_groups + if self.num_delayed_groups == 0: + mgxs.delayed_groups = None + else: + delayed_groups \ + = list(range(1,self.num_delayed_groups+1)) + mgxs.delayed_groups = delayed_groups # If a tally trigger was specified, add it to the MGXS if self.tally_trigger is not None: @@ -896,6 +886,7 @@ class Library(object): if nuclide != 'total': name += '_' + nuclide xsdata = openmc.XSdata(name, self.energy_groups) + xsdata.num_delayed_groups = self.num_delayed_groups if order is None: # Set the order to the Library's order (the defualt behavior) @@ -922,41 +913,79 @@ class Library(object): mymgxs = self.get_mgxs(domain, 'nu-transport') xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], subdomain=subdomain) + elif 'total' in self.mgxs_types: mymgxs = self.get_mgxs(domain, 'total') xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], subdomain=subdomain) + if 'absorption' in self.mgxs_types: mymgxs = self.get_mgxs(domain, 'absorption') xsdata.set_absorption_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], subdomain=subdomain) + if 'fission' in self.mgxs_types: mymgxs = self.get_mgxs(domain, 'fission') xsdata.set_fission_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], subdomain=subdomain) + if 'kappa-fission' in self.mgxs_types: mymgxs = self.get_mgxs(domain, 'kappa-fission') xsdata.set_kappa_fission_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], subdomain=subdomain) - # For chi and nu-fission we can either have only a nu-fission matrix - # provided, or vectors of chi and nu-fission provided + + if 'inverse-velocity' in self.mgxs_types: + mymgxs = self.get_mgxs(domain, 'inverse-velocity') + xsdata.set_inverse_velocity_mgxs(mymgxs, xs_type=xs_type, + nuclide=[nuclide], + subdomain=subdomain) + if 'nu-fission matrix' in self.mgxs_types: mymgxs = self.get_mgxs(domain, 'nu-fission matrix') xsdata.set_nu_fission_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], subdomain=subdomain) - else: - if 'chi' in self.mgxs_types: - mymgxs = self.get_mgxs(domain, 'chi') - xsdata.set_chi_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], - subdomain=subdomain) - if 'nu-fission' in self.mgxs_types: - mymgxs = self.get_mgxs(domain, 'nu-fission') - xsdata.set_nu_fission_mgxs(mymgxs, xs_type=xs_type, - nuclide=[nuclide], - subdomain=subdomain) + + if 'chi' in self.mgxs_types: + mymgxs = self.get_mgxs(domain, 'chi') + xsdata.set_chi_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], + subdomain=subdomain) + + if 'chi-prompt' in self.mgxs_types: + mymgxs = self.get_mgxs(domain, 'chi-prompt') + xsdata.set_chi_prompt_mgxs(mymgxs, xs_type=xs_type, + nuclide=[nuclide], subdomain=subdomain) + + if 'chi-delayed' in self.mgxs_types: + mymgxs = self.get_mgxs(domain, 'chi-delayed') + xsdata.set_chi_delayed_mgxs(mymgxs, xs_type=xs_type, + nuclide=[nuclide], subdomain=subdomain) + + if 'nu-fission' in self.mgxs_types: + mymgxs = self.get_mgxs(domain, 'nu-fission') + xsdata.set_nu_fission_mgxs(mymgxs, xs_type=xs_type, + nuclide=[nuclide], + subdomain=subdomain) + + if 'prompt-nu-fission' in self.mgxs_types: + mymgxs = self.get_mgxs(domain, 'prompt-nu-fission') + 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 'beta' in self.mgxs_types: + mymgxs = self.get_mgxs(domain, 'nu-fission') + xsdata.set_beta_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], + subdomain=subdomain) + # If multiplicity matrix is available, prefer that if 'multiplicity matrix' in self.mgxs_types: mymgxs = self.get_mgxs(domain, 'multiplicity matrix') @@ -964,6 +993,7 @@ class Library(object): nuclide=[nuclide], subdomain=subdomain) using_multiplicity = True + # multiplicity will fall back to using scatter and nu-scatter elif ((('scatter matrix' in self.mgxs_types) and ('nu-scatter matrix' in self.mgxs_types))): @@ -974,6 +1004,7 @@ class Library(object): nuclide=[nuclide], subdomain=subdomain) using_multiplicity = True + else: using_multiplicity = False @@ -995,8 +1026,8 @@ class Library(object): # absorption cross section. if 'total' in self.mgxs_types or 'transport' in self.mgxs_types: for i in range(len(xsdata.temperatures)): - xsdata._absorption[i] = \ - np.subtract(xsdata._total[i], np.sum( + xsdata._absorption[i] \ + = np.subtract(xsdata._total[i], np.sum( xsdata._scatter_matrix[i][0, :, :], axis=1)) return xsdata @@ -1051,7 +1082,9 @@ class Library(object): xs_type = 'macro' # Initialize file - mgxs_file = openmc.MGXSLibrary(self.energy_groups) + mgxs_file = openmc.MGXSLibrary(self.energy_groups, + num_delayed_groups=\ + self.num_delayed_groups) if self.domain_type == 'mesh': # Create the xsdata objects and add to the mgxs_file diff --git a/openmc/mgxs/mdgxs.py b/openmc/mgxs/mdgxs.py index 161d9e0f24..a99a99bee1 100644 --- a/openmc/mgxs/mdgxs.py +++ b/openmc/mgxs/mdgxs.py @@ -184,7 +184,7 @@ class MDGXS(MGXS): cv.check_less_than('delayed group', group, MAX_DELAYED_GROUPS, equality=True) - self._delayed_groups = delayed_groups + self._delayed_groups = delayed_groups @property def filters(self): @@ -1464,6 +1464,11 @@ class Beta(MDGXS): \beta_{d,g} &= \frac{\langle \nu^d \sigma_f \phi \rangle} {\langle \nu \sigma_f \phi \rangle} + NOTE: The Beta MGXS is the delayed neutron fraction computed directly from + the nuclear data. Often the delayed neutron fraction is + "importance-weighted" by the adjoint flux and called "beta-effective". It + is important to make clear that this Beta is not importance-weighted. + Parameters ---------- domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 2a8b0f6b76..adfabe022f 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -9,13 +9,14 @@ import h5py import openmc import openmc.mgxs from openmc.checkvalue import check_type, check_value, check_greater_than, \ - check_iterable_type + check_iterable_type, check_less_than # Supported incoming particle MGXS angular treatment representations _REPRESENTATIONS = ['isotropic', 'angle'] _SCATTER_TYPES = ['tabular', 'legendre', 'histogram'] -_SCATTER_SHAPES = ["[Order][G][G']"] +_XS_SHAPES = ["[Order][G][G']", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]", + "[G'][DG]"] class XSdata(object): @@ -34,6 +35,8 @@ class XSdata(object): temperatures : Iterable of float Temperatures (in units of Kelvin) of the provided datasets. Defaults to a single temperature at 294K. + num_delayed_groups : int + Number of delayed groups Attributes ---------- @@ -47,12 +50,12 @@ class XSdata(object): to a single temperature at 294K. energy_groups : openmc.mgxs.EnergyGroups Energy group structure + num_delayed_groups : int + Num delayed groups fissionable : bool Whether or not this is a fissionable data set. scatter_format : {'legendre', 'histogram', or 'tabular'} Angular distribution representation (legendre, histogram, or tabular) - scatter_shapes : {"[Order][G][G']"} - Dimensionality of the scattering and multiplicity matrices order : int Either the Legendre order, number of bins, or number of points used to describe the angular distribution associated with each group-to-group @@ -68,20 +71,6 @@ class XSdata(object): Number of equal width angular bins that the polar angular domain is subdivided into. This only applies when :attr:`XSdata.representation` is "angle". - use_chi : bool - Whether or not a chi vector or nu-fission matrix was used. - vector_shape : iterable of int - Dimensionality of vector multi-group cross sections (e.g., the total - cross section). The return result depends on the value of - :attr:`XSdata.representation`. - matrix_shape : iterable of int - Dimensionality of matrix multi-group cross sections (e.g., the - fission matrix cross section). The return result depends on the - value of :attr:`XSdata.representation`. - pn_matrix_shape : iterable of int - Dimensionality of scattering matrix data (e.g., the - scattering matrix cross section). The return result depends on the - value of :attr:`XSdata.representation`. total : dict of numpy.ndarray Group-wise total cross section. absorption : dict of numpy.ndarray @@ -98,7 +87,7 @@ class XSdata(object): fission : dict of numpy.ndarray Group-wise fission cross section. kappa_fission : dict of numpy.ndarray - Group-wise kappa-fission cross section. + Group-wise kappa_fission cross section. chi : dict of numpy.ndarray Group-wise fission spectra ordered by increasing group index (i.e., fast to thermal). This attribute should be used if making the common @@ -106,58 +95,93 @@ class XSdata(object): energy. If the user does not wish to make this approximation, then this should not be provided and this information included in the :attr:`XSdata.nu_fission` attribute instead. + chi_prompt : dict of numpy.ndarray + Group-wise prompt fission spectra ordered by increasing group index + (i.e., fast to thermal). This attribute should be used if chi from + prompt and delayed neutrons is being set separately. + chi_delayed : dict of numpy.ndarray + Group-wise delayed fission spectra ordered by increasing group index + (i.e., fast to thermal). This attribute should be used if chi from + prompt and delayed neutrons is being set separately. nu_fission : dict of numpy.ndarray - Group-wise fission production cross section vector (i.e., if - :attr:`XSdata.chi` is provided), or is the group-wise fission - production matrix. - inverse_velocities : dict of numpy.ndarray - Inverse of velocities, in units of sec/cm. + Group-wise fission production cross section vector (i.e., if ``chi`` is + provided), or is the group-wise fission production matrix. + prompt_nu_fission : dict of numpy.ndarray + Group-wise prompt fission production cross section vector. + delayed_nu_fission : dict of numpy.ndarray + Group-wise delayed fission production cross section vector. + beta : dict of numpy.ndarray + Delayed-group-wise delayed neutron fraction cross section vector. + decay_rate : dict of numpy.ndarray + Delayed-group-wise decay rate vector. + inverse_velocity : dict of numpy.ndarray + Inverse of velocity, in units of sec/cm. + xs_shapes : dict of iterable of int + Dictionary with keys of _XS_SHAPES and iterable of int values with the + corresponding shapes where "Order" corresponds to the pn scattering + order, "G" corresponds to incoming energy group, "G'" corresponds to + outgoing energy group, and "DG" corresponds to delayed group. Notes ----- The parameters containing cross section data have dimensionalities which depend upon the value of :attr:`XSdata.representation` as well as the number of Legendre or other angular dimensions as described by - :attr:`XSdata.order`. The :attr:`XSdata.vector_shape`, - :attr:`XSdata.matrix_shape`, and :attr:`XSdata.pn_matrix_shape` properties - are provided to obtain the dimensionality of the data for each temperature. + :attr:`XSdata.order`. The :attr:`XSdata.xs_shapes` are provided to obtain + the dimensionality of the data for each temperature. - The following are cross sections which should use each of these properties: + The following are cross sections which should use each of the properties. + Note that some cross sections can be input in more than one shape so they + are listed multiple times: - vector_shape : - total, absorption, fission, kappa_fission, chi, nu_fission (if chi is - provided), inverse_velocities - matrix_shape : - multiplicity_matrix, nu_fission (if chi is not provided) - pn_matrix_shape : - scatter_matrix + [Order][G][G']: scatter_matrix + + [G]: total, absorption, fission, kappa_fission, nu_fission, + prompt_nu_fission, inverse_velocity + + [G']: chi, chi_prompt, chi_delayed + + [G][G']: multiplicity_matrix, nu_fission + + [DG]: beta, decay_rate + + [G][DG]: delayed_nu_fission, beta, decay_rate + + [G'][DG]: chi_delayed """ def __init__(self, name, energy_groups, temperatures=[294.], - representation='isotropic'): + representation='isotropic', num_delayed_groups=0): + # Initialize class attributes self.name = name self.energy_groups = energy_groups + self.num_delayed_groups = num_delayed_groups self.temperatures = temperatures self.representation = representation self._atomic_weight_ratio = None self._fissionable = False self._scatter_format = 'legendre' - self._scatter_shape = "[Order][G][G']" self._order = None self._num_polar = None self._num_azimuthal = None - self._use_chi = None self._total = len(temperatures) * [None] self._absorption = len(temperatures) * [None] self._scatter_matrix = len(temperatures) * [None] self._multiplicity_matrix = len(temperatures) * [None] self._fission = len(temperatures) * [None] self._nu_fission = len(temperatures) * [None] + self._prompt_nu_fission = len(temperatures) * [None] + self._delayed_nu_fission = len(temperatures) * [None] self._kappa_fission = len(temperatures) * [None] self._chi = len(temperatures) * [None] - self._inverse_velocities = len(temperatures) * [None] + self._chi_prompt = len(temperatures) * [None] + self._chi_delayed = len(temperatures) * [None] + self._beta = len(temperatures) * [None] + self._decay_rate = len(temperatures) * [None] + self._inverse_velocity = len(temperatures) * [None] + self._xs_shapes = None @property def name(self): @@ -167,6 +191,10 @@ class XSdata(object): def energy_groups(self): return self._energy_groups + @property + def num_delayed_groups(self): + return self._num_delayed_groups + @property def representation(self): return self._representation @@ -187,10 +215,6 @@ class XSdata(object): def scatter_format(self): return self._scatter_format - @property - def scatter_shape(self): - return self._scatter_shape - @property def order(self): return self._order @@ -203,10 +227,6 @@ class XSdata(object): def num_azimuthal(self): return self._num_azimuthal - @property - def use_chi(self): - return self._use_chi - @property def total(self): return self._total @@ -231,6 +251,14 @@ class XSdata(object): def nu_fission(self): return self._nu_fission + @property + def prompt_nu_fission(self): + return self._prompt_nu_fission + + @property + def delayed_nu_fission(self): + return self._delayed_nu_fission + @property def kappa_fission(self): return self._kappa_fission @@ -239,6 +267,14 @@ class XSdata(object): def chi(self): return self._chi + @property + def chi_prompt(self): + return self._chi_prompt + + @property + def chi_delayed(self): + return self._chi_delayed + @property def num_orders(self): if self._order is not None: @@ -248,32 +284,31 @@ class XSdata(object): return self._order @property - def vector_shape(self): - if self.representation == 'isotropic': - return (self.energy_groups.num_groups,) - elif self.representation == 'angle': - return (self.num_polar, self.num_azimuthal, - self.energy_groups.num_groups) + def xs_shapes(self): - @property - def matrix_shape(self): - if self.representation == 'isotropic': - return (self.energy_groups.num_groups, - self.energy_groups.num_groups) - elif self.representation == 'angle': - return (self.num_polar, self.num_azimuthal, - self.energy_groups.num_groups, - self.energy_groups.num_groups) + if self._xs_shapes is None: - @property - def pn_matrix_shape(self): - if self.representation == 'isotropic': - return (self.num_orders, self.energy_groups.num_groups, - self.energy_groups.num_groups) - elif self.representation == 'angle': - return (self.num_polar, self.num_azimuthal, self.num_orders, - self.energy_groups.num_groups, - self.energy_groups.num_groups) + self._xs_shapes = {} + self._xs_shapes["[G]"] = (self.energy_groups.num_groups,) + self._xs_shapes["[G']"] = (self.energy_groups.num_groups,) + self._xs_shapes["[G][G']"] = (self.energy_groups.num_groups, + self.energy_groups.num_groups) + self._xs_shapes["[DG]"] = (self.num_delayed_groups,) + self._xs_shapes["[G][DG]"] = (self.energy_groups.num_groups, + self.num_delayed_groups) + self._xs_shapes["[G'][DG]"] = (self.energy_groups.num_groups, + self.num_delayed_groups) + self._xs_shapes["[Order][G][G']"] \ + = (self.num_orders, self.energy_groups.num_groups, + self.energy_groups.num_groups) + + # If representation is by angle prepend num polar and num azim + if self.representation == 'angle': + for key,shapes in self._xs_shapes.items(): + self._xs_shapes[key] \ + = (self.num_polar, self.num_azimuthal) + shapes + + return self._xs_shapes @name.setter def name(self, name): @@ -282,6 +317,7 @@ class XSdata(object): @energy_groups.setter def energy_groups(self, energy_groups): + # Check validity of energy_groups check_type('energy_groups', energy_groups, openmc.mgxs.EnergyGroups) @@ -289,16 +325,30 @@ class XSdata(object): msg = 'Unable to assign an EnergyGroups object ' \ 'with uninitialized group edges' raise ValueError(msg) + self._energy_groups = energy_groups + @num_delayed_groups.setter + def num_delayed_groups(self, num_delayed_groups): + + # Check validity of num_delayed_groups + check_type('num_delayed_groups', num_delayed_groups, int) + check_less_than('num_delayed_groups', num_delayed_groups, + openmc.mgxs.MAX_DELAYED_GROUPS, equality=True) + check_greater_than('num_delayed_groups', num_delayed_groups, 0, + equality=True) + self._num_delayed_groups = num_delayed_groups + @representation.setter def representation(self, representation): + # Check it is of valid type. check_value('representation', representation, _REPRESENTATIONS) self._representation = representation @atomic_weight_ratio.setter def atomic_weight_ratio(self, atomic_weight_ratio): + # Check validity of type and that the atomic_weight_ratio value is > 0 check_type('atomic_weight_ratio', atomic_weight_ratio, Real) check_greater_than('atomic_weight_ratio', atomic_weight_ratio, 0.0) @@ -306,24 +356,20 @@ class XSdata(object): @temperatures.setter def temperatures(self, temperatures): - check_iterable_type('temperatures', temperatures, Real) + check_iterable_type('temperatures', temperatures, Real) self._temperatures = np.array(temperatures) @scatter_format.setter def scatter_format(self, scatter_format): + # check to see it is of a valid type and value check_value('scatter_format', scatter_format, _SCATTER_TYPES) self._scatter_format = scatter_format - @scatter_shape.setter - def scatter_shape(self, scatter_shape): - # check to see it is of a valid type and value - check_value('scatter_shape', scatter_shape, _SCATTER_SHAPES) - self._scatter_shape = scatter_shape - @order.setter def order(self, order): + # Check type and value check_type('order', order, Integral) check_greater_than('order', order, 0, equality=True) @@ -331,6 +377,7 @@ class XSdata(object): @num_polar.setter def num_polar(self, num_polar): + # Make sure we have positive ints check_type('num_polar', num_polar, Integral) check_greater_than('num_polar', num_polar, 0) @@ -338,15 +385,11 @@ class XSdata(object): @num_azimuthal.setter def num_azimuthal(self, num_azimuthal): + check_type('num_azimuthal', num_azimuthal, Integral) check_greater_than('num_azimuthal', num_azimuthal, 0) self._num_azimuthal = num_azimuthal - @use_chi.setter - def use_chi(self, use_chi): - check_type('use_chi', use_chi, bool) - self._use_chi = use_chi - def add_temperature(self, temperature): """This method re-sizes the attributes of this XSdata object so that it can accomodate an additional temperature. Note that the set_* methods @@ -358,6 +401,7 @@ class XSdata(object): Temperature (in units of Kelvin) of the provided dataset. """ + check_type('temperature', temperature, Real) temp_store = self.temperatures.tolist().append(temperature) @@ -369,9 +413,15 @@ class XSdata(object): self._multiplicity_matrix.append(None) self._fission.append(None) self._nu_fission.append(None) + self._prompt_nu_fission.append(None) + self._delayed_nu_fission.append(None) self._kappa_fission.append(None) self._chi.append(None) - self._inverse_velocities.append(None) + self._chi_prompt.append(None) + self._chi_delayed.append(None) + self._beta.append(None) + self._decay_rate.append(None) + self._inverse_velocity.append(None) def set_total(self, total, temperature=294.): """This method sets the cross section for this XSdata object at the @@ -390,10 +440,13 @@ class XSdata(object): openmc.mgxs_library.set_total_mgxs() """ - check_type('total', total, Iterable, expected_iter_type=Real) + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G]"]] + # Convert to a numpy array so we can easily get the shape for checking total = np.asarray(total) - check_value('total shape', total.shape, [self.vector_shape]) + check_value('total shape', total.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -417,10 +470,13 @@ class XSdata(object): openmc.mgxs_library.set_absorption_mgxs() """ - check_type('absorption', absorption, Iterable, expected_iter_type=Real) + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G]"]] + # Convert to a numpy array so we can easily get the shape for checking absorption = np.asarray(absorption) - check_value('absorption shape', absorption.shape, [self.vector_shape]) + check_value('absorption shape', absorption.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -444,10 +500,13 @@ class XSdata(object): openmc.mgxs_library.set_fission_mgxs() """ - check_type('fission', fission, Iterable, expected_iter_type=Real) + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G]"]] + # Convert to a numpy array so we can easily get the shape for checking fission = np.asarray(fission) - check_value('fission shape', fission.shape, [self.vector_shape]) + check_value('fission shape', fission.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -474,12 +533,13 @@ class XSdata(object): openmc.mgxs_library.set_kappa_fission_mgxs() """ - check_type('kappa_fission', kappa_fission, Iterable, - expected_iter_type=Real) + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G]"]] + # Convert to a numpy array so we can easily get the shape for checking kappa_fission = np.asarray(kappa_fission) - check_value('kappa fission shape', kappa_fission.shape, - [self.vector_shape]) + check_value('kappa fission shape', kappa_fission.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -506,27 +566,138 @@ class XSdata(object): openmc.mgxs_library.set_chi_mgxs() """ - if self.use_chi is not None: - if not self.use_chi: - msg = 'Providing "chi" when "nu-fission" already provided ' \ - 'as a matrix' - raise ValueError(msg) - check_type('chi', chi, Iterable, expected_iter_type=Real) + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G']"]] + # Convert to a numpy array so we can easily get the shape for checking chi = np.asarray(chi) - # Check the shape - if chi.shape != self.vector_shape: - msg = 'Provided chi iterable does not have the expected shape.' - raise ValueError(msg) + check_value('chi shape', chi.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) i = np.where(self.temperatures == temperature)[0][0] self._chi[i] = chi - if self.use_chi is not None: - self.use_chi = True + def set_chi_prompt(self, chi_prompt, temperature=294.): + """This method sets the cross section for this XSdata object at the + provided temperature. + + Parameters + ---------- + chi_prompt : np.ndarray + Prompt fission Spectrum + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + + See also + -------- + openmc.mgxs_library.set_chi_prompt_mgxs() + + """ + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G']"]] + + # Convert to a numpy array so we can easily get the shape for checking + chi_prompt = np.asarray(chi_prompt) + check_value('chi prompt shape', chi_prompt.shape, shapes) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + self._chi_prompt[i] = chi_prompt + + def set_chi_delayed(self, chi_delayed, temperature=294.): + """This method sets the cross section for this XSdata object at the + provided temperature. + + Parameters + ---------- + chi_delayed : np.ndarray + Delayed fission Spectrum + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + + See also + -------- + openmc.mgxs_library.set_chi_delayed_mgxs() + + """ + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G']"], self.xs_shapes["[G'][DG]"]] + + # Convert to a numpy array so we can easily get the shape for checking + chi_delayed = np.asarray(chi_delayed) + check_value('chi delayed shape', chi_delayed.shape, shapes) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + self._chi_delayed[i] = chi_delayed + + def set_beta(self, beta, temperature=294.): + """This method sets the cross section for this XSdata object at the + provided temperature. + + Parameters + ---------- + beta : np.ndarray + Delayed fission spectrum + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + + See also + -------- + openmc.mgxs_library.set_beta_mgxs() + + """ + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[DG]"], self.xs_shapes["[G][DG]"]] + + # Convert to a numpy array so we can easily get the shape for checking + beta = np.asarray(beta) + check_value('beta shape', beta.shape, shapes) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + self._beta[i] = beta + + def set_decay_rate(self, decay_rate, temperature=294.): + """This method sets the cross section for this XSdata object at the + provided temperature. + + Parameters + ---------- + decay_rate : np.ndarray + Delayed neutron precursor decay rate + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + + See also + -------- + openmc.mgxs_library.set_decay_rate_mgxs() + + """ + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[DG]"], self.xs_shapes["[G][DG]"]] + + # Convert to a numpy array so we can easily get the shape for checking + decay_rate = np.asarray(decay_rate) + check_value('decay rate shape', decay_rate.shape, shapes) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + self._decay_rate[i] = decay_rate def set_scatter_matrix(self, scatter, temperature=294.): """This method sets the cross section for this XSdata object at the @@ -545,11 +716,15 @@ class XSdata(object): openmc.mgxs_library.set_scatter_matrix_mgxs() """ + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[Order][G][G']"]] + # Convert to a numpy array so we can easily get the shape for checking scatter = np.asarray(scatter) check_iterable_type('scatter', scatter, Real, max_depth=len(scatter.shape)) - check_value('scatter shape', scatter.shape, [self.pn_matrix_shape]) + check_value('scatter shape', scatter.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -573,12 +748,15 @@ class XSdata(object): openmc.mgxs_library.set_multiplicity_matrix_mgxs() """ + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G][G']"]] + # Convert to a numpy array so we can easily get the shape for checking multiplicity = np.asarray(multiplicity) check_iterable_type('multiplicity', multiplicity, Real, max_depth=len(multiplicity.shape)) - check_value('multiplicity shape', multiplicity.shape, - [self.matrix_shape]) + check_value('multiplicity shape', multiplicity.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -602,82 +780,129 @@ class XSdata(object): openmc.mgxs_library.set_nu_fission_mgxs() """ - # The NuFissionXS class does not have the capability to produce - # a fission matrix and therefore if this path is pursued, we know - # chi must be used. - # nu_fission can be given as a vector or a matrix - # Vector is used when chi also exists. - # Matrix is used when chi does not exist. - # We have to check that the correct form is given, but only if - # chi already has been set. If not, we just check that this is OK - # and set the use_chi flag accordingly + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G]"], self.xs_shapes["[G][G']"]] # Convert to a numpy array so we can easily get the shape for checking nu_fission = np.asarray(nu_fission) - + check_value('nu_fission shape', nu_fission.shape, shapes) check_iterable_type('nu_fission', nu_fission, Real, max_depth=len(nu_fission.shape)) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) - if self.use_chi is not None: - if self.use_chi: - check_value('nu_fission shape', nu_fission.shape, - [self.vector_shape]) - else: - check_value('nu_fission shape', nu_fission.shape, - [self.matrix_shape]) - else: - check_value('nu_fission shape', nu_fission.shape, - [self.vector_shape, self.matrix_shape]) - # Find out if we have a nu-fission matrix or vector - # and set a flag to allow other methods to check this later. - if nu_fission.shape == self.vector_shape: - self.use_chi = True - else: - self.use_chi = False - i = np.where(self.temperatures == temperature)[0][0] self._nu_fission[i] = nu_fission if np.sum(nu_fission) > 0.0: self._fissionable = True - def set_inverse_velocities(self, inv_vel, temperature=294.): - """This method sets the inverse velocities for this XSdata object at the + def set_prompt_nu_fission(self, prompt_nu_fission, temperature=294.): + """This method sets the cross section for this XSdata object at the + provided temperature. + + Parameters + ---------- + prompt_nu_fission: np.ndarray + Prompt-nu-fission Cross Section + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + + See also + -------- + openmc.mgxs_library.set_prompt_nu_fission_mgxs() + + """ + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G]"]] + + # Convert to a numpy array so we can easily get the shape for checking + prompt_nu_fission = np.asarray(prompt_nu_fission) + check_value('prompt_nu_fission shape', prompt_nu_fission.shape, shapes) + check_iterable_type('prompt_nu_fission', prompt_nu_fission, Real, + max_depth=len(prompt_nu_fission.shape)) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + self._prompt_nu_fission[i] = prompt_nu_fission + if np.sum(prompt_nu_fission) > 0.0: + self._fissionable = True + + def set_delayed_nu_fission(self, delayed_nu_fission, temperature=294.): + """This method sets the cross section for this XSdata object at the + provided temperature. + + Parameters + ---------- + delayed_nu_fission: np.ndarray + Delayed-nu-fission Cross Section + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + + See also + -------- + openmc.mgxs_library.set_delayed_nu_fission_mgxs() + + """ + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G][DG]"]] + + # Convert to a numpy array so we can easily get the shape for checking + delayed_nu_fission = np.asarray(delayed_nu_fission) + check_value('delayed_nu_fission shape', delayed_nu_fission.shape, + shapes) + check_iterable_type('delayed_nu_fission', delayed_nu_fission, Real, + max_depth=len(delayed_nu_fission.shape)) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + self._delayed_nu_fission[i] = delayed_nu_fission + if np.sum(delayed_nu_fission) > 0.0: + self._fissionable = True + + def set_inverse_velocity(self, inv_vel, temperature=294.): + """This method sets the inverse velocity for this XSdata object at the provided temperature. Parameters ---------- inv_vel: np.ndarray - Inverse velocities in units of sec/cm. + Inverse velocity in units of sec/cm. temperature : float Temperature (in Kelvin) of the data. Defaults to room temperature (294K). """ - check_type('inverse velocities', inv_vel, Iterable, - expected_iter_type=Real) + + # Get the accepted shapes for this xs + shapes = [self.xs_shapes["[G]"]] + # Convert to a numpy array so we can easily get the shape for checking inv_vel = np.asarray(inv_vel) - check_value('inverse velocities shape', inv_vel.shape, - [self.vector_shape]) + check_value('inverse_velocity shape', inv_vel.shape, shapes) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) i = np.where(self.temperatures == temperature)[0][0] - self._inverse_velocities[i] = inv_vel + self._inverse_velocity[i] = inv_vel def set_total_mgxs(self, total, temperature=294., nuclide='total', xs_type='macro', subdomain=None): """This method allows for an openmc.mgxs.TotalXS or - openmc.mgxs.TransportXS to be used to set the total cross section - for this XSdata object. + openmc.mgxs.TransportXS to be used to set the total cross section for + this XSdata object. Parameters ---------- total: openmc.mgxs.TotalXS or openmc.mgxs.TransportXS - MGXS Object containing the total or transport cross section - for the domain of interest. + MGXS Object containing the total, transport or nu-transport cross + section for the domain of interest. temperature : float Temperature (in Kelvin) of the data. Defaults to room temperature (294K). @@ -701,8 +926,7 @@ class XSdata(object): check_type('total', total, (openmc.mgxs.TotalXS, openmc.mgxs.TransportXS)) check_value('energy_groups', total.energy_groups, [self.energy_groups]) - check_value('domain_type', total.domain_type, - ['universe', 'cell', 'material', 'mesh']) + check_value('domain_type', total.domain_type, openmc.mgxs.DOMAIN_TYPES) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -748,7 +972,7 @@ class XSdata(object): check_value('energy_groups', absorption.energy_groups, [self.energy_groups]) check_value('domain_type', absorption.domain_type, - ['universe', 'cell', 'material', 'mesh']) + openmc.mgxs.DOMAIN_TYPES) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -795,7 +1019,7 @@ class XSdata(object): check_value('energy_groups', fission.energy_groups, [self.energy_groups]) check_value('domain_type', fission.domain_type, - ['universe', 'cell', 'material', 'mesh']) + openmc.mgxs.DOMAIN_TYPES) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -843,7 +1067,7 @@ class XSdata(object): check_value('energy_groups', nu_fission.energy_groups, [self.energy_groups]) check_value('domain_type', nu_fission.domain_type, - ['universe', 'cell', 'material', 'mesh']) + openmc.mgxs.DOMAIN_TYPES) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -856,14 +1080,117 @@ class XSdata(object): msg = 'Angular-Dependent MGXS have not yet been implemented' raise ValueError(msg) - if isinstance(nu_fission, openmc.mgxs.NuFissionMatrixXS): - self.use_chi = False - else: - self.use_chi = True - if np.sum(self._nu_fission) > 0.0: self._fissionable = True + 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. + + Parameters + ---------- + prompt_nu_fission: openmc.mgxs.PromptNuFissionXS + MGXS Object containing the prompt-nu-fission cross section + for the domain of interest. + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + nuclide : str + Individual nuclide (or 'total' if obtaining material-wise data) + to gather data for. Defaults to 'total'. + xs_type: {'macro', 'micro'} + Provide the macro or micro cross section in units of cm^-1 or + barns. Defaults to 'macro'. + subdomain : iterable of int + If the MGXS contains a mesh domain type, the subdomain parameter + specifies which mesh cell (i.e., [i, j, k] index) to use. + + See also + -------- + openmc.mgxs.Library.create_mg_library() + openmc.mgxs.Library.get_xsdata + + """ + + check_type('prompt_nu_fission', prompt_nu_fission, + (openmc.mgxs.PromptNuFissionXS,)) + check_value('energy_groups', prompt_nu_fission.energy_groups, + [self.energy_groups]) + check_value('domain_type', prompt_nu_fission.domain_type, + openmc.mgxs.DOMAIN_TYPES) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + if self.representation is 'isotropic': + self._prompt_nu_fission[i] = prompt_nu_fission.get_xs\ + (nuclides=nuclide, + xs_type=xs_type, + subdomains=subdomain) + elif self.representation is 'angle': + msg = 'Angular-Dependent MGXS have not yet been implemented' + raise ValueError(msg) + + if np.sum(self._prompt_nu_fission) > 0.0: + self._fissionable = True + + 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. + + Parameters + ---------- + delayed_nu_fission: openmc.mgxs.DelayedNuFissionXS + MGXS Object containing the delayed-nu-fission cross section + for the domain of interest. + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + nuclide : str + Individual nuclide (or 'total' if obtaining material-wise data) + to gather data for. Defaults to 'total'. + xs_type: {'macro', 'micro'} + Provide the macro or micro cross section in units of cm^-1 or + barns. Defaults to 'macro'. + subdomain : iterable of int + If the MGXS contains a mesh domain type, the subdomain parameter + specifies which mesh cell (i.e., [i, j, k] index) to use. + + See also + -------- + openmc.mgxs.Library.create_mg_library() + openmc.mgxs.Library.get_xsdata + + """ + + check_type('delayed_nu_fission', delayed_nu_fission, + (openmc.mgxs.DelayedNuFissionXS,)) + check_value('energy_groups', delayed_nu_fission.energy_groups, + [self.energy_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) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + if self.representation is 'isotropic': + self._delayed_nu_fission[i] = delayed_nu_fission.get_xs\ + (nuclides=nuclide, + xs_type=xs_type, + subdomains=subdomain) + elif self.representation is 'angle': + msg = 'Angular-Dependent MGXS have not yet been implemented' + raise ValueError(msg) + + if np.sum(self._delayed_nu_fission) > 0.0: + self._fissionable = True + def set_kappa_fission_mgxs(self, k_fission, temperature=294., nuclide='total', xs_type='macro', subdomain=None): @@ -896,11 +1223,11 @@ class XSdata(object): """ - check_type('k_fission', k_fission, openmc.mgxs.KappaFissionXS) + check_type('kappa_fission', k_fission, openmc.mgxs.KappaFissionXS) check_value('energy_groups', k_fission.energy_groups, [self.energy_groups]) check_value('domain_type', k_fission.domain_type, - ['universe', 'cell', 'material', 'mesh']) + openmc.mgxs.DOMAIN_TYPES) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -942,16 +1269,9 @@ class XSdata(object): """ - if self.use_chi is not None: - if not self.use_chi: - msg = 'Providing chi when nu_fission already provided as a ' \ - 'matrix!' - raise ValueError(msg) - check_type('chi', chi, openmc.mgxs.Chi) check_value('energy_groups', chi.energy_groups, [self.energy_groups]) - check_value('domain_type', chi.domain_type, - ['universe', 'cell', 'material', 'mesh']) + check_value('domain_type', chi.domain_type, openmc.mgxs.DOMAIN_TYPES) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -963,8 +1283,190 @@ class XSdata(object): msg = 'Angular-Dependent MGXS have not yet been implemented' raise ValueError(msg) - if self.use_chi is not None: - self.use_chi = True + def set_chi_prompt_mgxs(self, chi_prompt, temperature=294., nuclide='total', + xs_type='macro', subdomain=None): + """This method allows for an openmc.mgxs.ChiPrompt + to be used to set chi-prompt for this XSdata object. + + Parameters + ---------- + chi_prompt: openmc.mgxs.ChiPrompt + MGXS Object containing chi-prompt for the domain of interest. + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + nuclide : str + Individual nuclide (or 'total' if obtaining material-wise data) + to gather data for. Defaults to 'total'. + xs_type: {'macro', 'micro'} + Provide the macro or micro cross section in units of cm^-1 or + barns. Defaults to 'macro'. + subdomain : iterable of int + If the MGXS contains a mesh domain type, the subdomain parameter + specifies which mesh cell (i.e., [i, j, k] index) to use. + + See also + -------- + openmc.mgxs.Library.create_mg_library() + openmc.mgxs.Library.get_xsdata + + """ + + check_type('chi_prompt', chi_prompt, openmc.mgxs.ChiPrompt) + check_value('energy_groups', chi_prompt.energy_groups, + [self.energy_groups]) + check_value('domain_type', chi_prompt.domain_type, + openmc.mgxs.DOMAIN_TYPES) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + if self.representation is 'isotropic': + self._chi_prompt[i] = chi_prompt.get_xs(nuclides=nuclide, + xs_type=xs_type, + subdomains=subdomain) + elif self.representation is 'angle': + msg = 'Angular-Dependent MGXS have not yet been implemented' + raise ValueError(msg) + + def set_chi_delayed_mgxs(self, chi_delayed, temperature=294., + nuclide='total', xs_type='macro', subdomain=None): + """This method allows for an openmc.mgxs.ChiDelayed + to be used to set chi-delayed for this XSdata object. + + Parameters + ---------- + chi_delayed: openmc.mgxs.ChiDelayed + MGXS Object containing chi-delayed for the domain of interest. + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + nuclide : str + Individual nuclide (or 'total' if obtaining material-wise data) + to gather data for. Defaults to 'total'. + xs_type: {'macro', 'micro'} + Provide the macro or micro cross section in units of cm^-1 or + barns. Defaults to 'macro'. + subdomain : iterable of int + If the MGXS contains a mesh domain type, the subdomain parameter + specifies which mesh cell (i.e., [i, j, k] index) to use. + + See also + -------- + openmc.mgxs.Library.create_mg_library() + openmc.mgxs.Library.get_xsdata + + """ + + check_type('chi_delayed', chi_delayed, openmc.mgxs.ChiDelayed) + check_value('energy_groups', chi_delayed.energy_groups, + [self.energy_groups]) + check_value('num_delayed_groups', chi_delayed.num_delayed_groups, + [self.num_delayed_groups]) + check_value('domain_type', chi_delayed.domain_type, + openmc.mgxs.DOMAIN_TYPES) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + if self.representation is 'isotropic': + self._chi_delayed[i] = chi_delayed.get_xs(nuclides=nuclide, + xs_type=xs_type, + subdomains=subdomain) + elif self.representation is 'angle': + msg = 'Angular-Dependent MGXS have not yet been implemented' + raise ValueError(msg) + + def set_beta_mgxs(self, beta, temperature=294., + nuclide='total', xs_type='macro', subdomain=None): + """This method allows for an openmc.mgxs.Beta + to be used to set beta for this XSdata object. + + Parameters + ---------- + beta : openmc.mgxs.Beta + MGXS Object containing beta for the domain of interest. + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + nuclide : str + Individual nuclide (or 'total' if obtaining material-wise data) + to gather data for. Defaults to 'total'. + xs_type: {'macro', 'micro'} + Provide the macro or micro cross section in units of cm^-1 or + barns. Defaults to 'macro'. + subdomain : iterable of int + If the MGXS contains a mesh domain type, the subdomain parameter + specifies which mesh cell (i.e., [i, j, k] index) to use. + + See also + -------- + openmc.mgxs.Library.create_mg_library() + openmc.mgxs.Library.get_xsdata + + """ + + check_type('beta', beta, openmc.mgxs.Beta) + check_value('num_delayed_groups', beta.num_delayed_groups, + [self.num_delayed_groups]) + check_value('domain_type', beta.domain_type, openmc.mgxs.DOMAIN_TYPES) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + if self.representation is 'isotropic': + self._beta[i] = beta.get_xs(nuclides=nuclide, + xs_type=xs_type, + subdomains=subdomain) + elif self.representation is 'angle': + msg = 'Angular-Dependent MGXS have not yet been implemented' + raise ValueError(msg) + + def set_decay_rate_mgxs(self, decay_rate, temperature=294., + nuclide='total', xs_type='macro', subdomain=None): + """This method allows for an openmc.mgxs.DecayRate + to be used to set decay rate for this XSdata object. + + Parameters + ---------- + decay_rate : openmc.mgxs.DecayRate + MGXS Object containing decay rate for the domain of interest. + temperature : float + Temperature (in units of Kelvin) of the provided dataset. Defaults + to 294K + nuclide : str + Individual nuclide (or 'total' if obtaining material-wise data) + to gather data for. Defaults to 'total'. + xs_type: {'macro', 'micro'} + Provide the macro or micro cross section in units of cm^-1 or + barns. Defaults to 'macro'. + subdomain : iterable of int + If the MGXS contains a mesh domain type, the subdomain parameter + specifies which mesh cell (i.e., [i, j, k] index) to use. + + See also + -------- + openmc.mgxs.Library.create_mg_library() + openmc.mgxs.Library.get_xsdata + + """ + + check_type('decay_rate', decay_rate, openmc.mgxs.DecayRate) + check_value('num_delayed_groups', decay_rate.num_delayed_groups, + [self.num_delayed_groups]) + check_value('domain_type', decay_rate.domain_type, + openmc.mgxs.DOMAIN_TYPES) + check_type('temperature', temperature, Real) + check_value('temperature', temperature, self.temperatures) + + i = np.where(self.temperatures == temperature)[0][0] + if self.representation is 'isotropic': + self._decay_rate[i] = decay_rate.get_xs(nuclides=nuclide, + xs_type=xs_type, + subdomains=subdomain) + elif self.representation is 'angle': + msg = 'Angular-Dependent MGXS have not yet been implemented' + raise ValueError(msg) def set_scatter_matrix_mgxs(self, scatter, temperature=294., nuclide='total', xs_type='macro', @@ -1003,7 +1505,7 @@ class XSdata(object): check_value('energy_groups', scatter.energy_groups, [self.energy_groups]) check_value('domain_type', scatter.domain_type, - ['universe', 'cell', 'material', 'mesh']) + openmc.mgxs.DOMAIN_TYPES) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -1085,7 +1587,7 @@ class XSdata(object): check_value('energy_groups', nuscatter.energy_groups, [self.energy_groups]) check_value('domain_type', nuscatter.domain_type, - ['universe', 'cell', 'material', 'mesh']) + openmc.mgxs.DOMAIN_TYPES) check_type('temperature', temperature, Real) check_value('temperature', temperature, self.temperatures) @@ -1099,8 +1601,7 @@ class XSdata(object): check_value('energy_groups', scatter.energy_groups, [self.energy_groups]) check_value('domain_type', scatter.domain_type, - ['universe', 'cell', 'material', 'mesh']) - + openmc.mgxs.DOMAIN_TYPES) i = np.where(self.temperatures == temperature)[0][0] if self.representation == 'isotropic': nuscatt = nuscatter.get_xs(nuclides=nuclide, @@ -1133,17 +1634,19 @@ class XSdata(object): grp.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio if self.fissionable is not None: grp.attrs['fissionable'] = self.fissionable + if self.representation is not None: grp.attrs['representation'] = np.string_(self.representation) if self.representation == 'angle': if self.num_azimuthal is not None: grp.attrs['num_azimuthal'] = self.num_azimuthal + if self.num_polar is not None: grp.attrs['num_polar'] = self.num_polar + + grp.attrs['scatter_shape'] = np.string_("[Order][G][G']") if self.scatter_format is not None: grp.attrs['scatter_format'] = np.string_(self.scatter_format) - if self.scatter_shape is not None: - grp.attrs['scatter_shape'] = np.string_(self.scatter_shape) if self.order is not None: grp.attrs['order'] = self.order @@ -1155,27 +1658,63 @@ class XSdata(object): # Create the temperature datasets for i, temperature in enumerate(self.temperatures): + xs_grp = grp.create_group(str(int(np.round(temperature))) + "K") + if self._total[i] is None: raise ValueError('total data must be provided when writing ' 'the HDF5 library') + xs_grp.create_dataset("total", data=self._total[i]) + if self._absorption[i] is None: raise ValueError('absorption data must be provided when ' 'writing the HDF5 library') + xs_grp.create_dataset("absorption", data=self._absorption[i]) + if self.fissionable: if self._fission[i] is not None: xs_grp.create_dataset("fission", data=self._fission[i]) + if self._kappa_fission[i] is not None: xs_grp.create_dataset("kappa-fission", data=self._kappa_fission[i]) + if self._chi[i] is not None: xs_grp.create_dataset("chi", data=self._chi[i]) - if self._nu_fission[i] is None: - raise ValueError('nu-fission data must be provided when ' - 'writing the HDF5 library') - xs_grp.create_dataset("nu-fission", data=self._nu_fission[i]) + + if self._chi_prompt[i] is not None: + xs_grp.create_dataset("chi-prompt", + data=self._chi_prompt[i]) + + if self._chi_delayed[i] is not None: + xs_grp.create_dataset("chi-delayed", + data=self._chi_delayed[i]) + + if self._nu_fission[i] is None and \ + (self._delayed_nu_fission[i] is None or \ + self._prompt_nu_fission[i] is None): + raise ValueError('nu-fission or prompt-nu-fission and ' + 'delayed-nu-fission data must be provided ' + 'when writing the HDF5 library') + + if self._nu_fission[i] is not None: + xs_grp.create_dataset("nu-fission", data=self._nu_fission[i]) + + if self._prompt_nu_fission[i] is not None: + xs_grp.create_dataset("prompt-nu-fission", + data=self._prompt_nu_fission[i]) + + if self._delayed_nu_fission[i] is not None: + xs_grp.create_dataset("delayed-nu-fission", + data=self._delayed_nu_fission[i]) + + if self._beta[i] is not None: + xs_grp.create_dataset("beta", data=self._beta[i]) + + if self._decay_rate[i] is not None: + xs_grp.create_dataset("decay rate", data=self._decay_rate[i]) if self._scatter_matrix[i] is None: raise ValueError('Scatter matrix must be provided when ' @@ -1184,11 +1723,14 @@ class XSdata(object): # Get the sparse scattering data to print to the library G = self.energy_groups.num_groups if self.representation == 'isotropic': + g_out_bounds = np.zeros((G, 2), dtype=np.int) + for g_in in range(G): nz = np.nonzero(self._scatter_matrix[i][0, g_in, :]) g_out_bounds[g_in, 0] = nz[0][0] g_out_bounds[g_in, 1] = nz[0][-1] + # Now create the flattened scatter matrix array matrix = self._scatter_matrix[i] flat_scatt = [] @@ -1197,12 +1739,15 @@ class XSdata(object): g_out_bounds[g_in, 1] + 1): for l in range(len(matrix[:, g_in, g_out])): flat_scatt.append(matrix[l, g_in, g_out]) + # And write it. scatt_grp = xs_grp.create_group('scatter_data') scatt_grp.create_dataset("scatter_matrix", data=np.array(flat_scatt)) + # Repeat for multiplicity if self._multiplicity_matrix[i] is not None: + # Now create the flattened scatter matrix array matrix = self._multiplicity_matrix[i][:, :] flat_mult = [] @@ -1210,7 +1755,8 @@ class XSdata(object): for g_out in range(g_out_bounds[g_in, 0], g_out_bounds[g_in, 1] + 1): flat_mult.append(matrix[g_in, g_out]) - scatt_grp.create_dataset("multiplicity_matrix", + + scatt_grp.create_dataset("multiplicity matrix", data=np.array(flat_mult)) # And finally, adjust g_out_bounds for 1-based group counting @@ -1223,6 +1769,7 @@ class XSdata(object): Np = self.num_polar Na = self.num_azimuthal g_out_bounds = np.zeros((Np, Na, G, 2), dtype=np.int) + for p in range(Np): for a in range(Na): for g_in in range(G): @@ -1230,6 +1777,7 @@ class XSdata(object): nz = np.nonzero(matrix) g_out_bounds[p, a, g_in, 0] = nz[0][0] g_out_bounds[p, a, g_in, 1] = nz[0][-1] + # Now create the flattened scatter matrix array flat_scatt = [] for p in range(Np): @@ -1240,12 +1788,15 @@ class XSdata(object): g_out_bounds[p, a, g_in, 1] + 1): for l in range(len(matrix[:, g_in, g_out])): flat_scatt.append(matrix[l, g_in, g_out]) + # And write it. scatt_grp = xs_grp.create_group('scatter_data') scatt_grp.create_dataset("scatter_matrix", data=np.array(flat_scatt)) + # Repeat for multiplicity if self._multiplicity_matrix[i] is not None: + # Now create the flattened scatter matrix array flat_mult = [] for p in range(Np): @@ -1255,9 +1806,10 @@ class XSdata(object): for g_out in range(g_out_bounds[p, a, g_in, 0], g_out_bounds[p, a, g_in, 1] + 1): flat_mult.append(matrix[g_in, g_out]) - # And write it. - scatt_grp.create_dataset("multiplicity_matrix", - data=np.array(flat_mult)) + + # And write it. + scatt_grp.create_dataset("multiplicity_matrix", + data=np.array(flat_mult)) # And finally, adjust g_out_bounds for 1-based group counting # and write it. @@ -1266,9 +1818,9 @@ class XSdata(object): scatt_grp.create_dataset("g_max", data=g_out_bounds[:, :, :, 1]) # Add the kinetics data - if self._inverse_velocities[i] is not None: - xs_grp.create_dataset("inverse_velocities", - data=self._inverse_velocities[i]) + if self._inverse_velocity[i] is not None: + xs_grp.create_dataset("inverse-velocity", + data=self._inverse_velocity[i]) class MGXSLibrary(object): @@ -1280,23 +1832,32 @@ class MGXSLibrary(object): ---------- energy_groups : openmc.mgxs.EnergyGroups Energy group structure + num_delayed_groups : int + Num delayed groups Attributes ---------- energy_groups : openmc.mgxs.EnergyGroups Energy group structure. + num_delayed_groups : int + Num delayed groups xsdatas : Iterable of openmc.XSdata Iterable of multi-Group cross section data objects """ - def __init__(self, energy_groups): + def __init__(self, energy_groups, num_delayed_groups=0): self.energy_groups = energy_groups + self.num_delayed_groups = num_delayed_groups self._xsdatas = [] @property def energy_groups(self): return self._energy_groups + @property + def num_delayed_groups(self): + return self._num_delayed_groups + @property def temperatures(self): return self._temperatures @@ -1310,6 +1871,15 @@ class MGXSLibrary(object): check_type('energy groups', energy_groups, openmc.mgxs.EnergyGroups) self._energy_groups = energy_groups + @num_delayed_groups.setter + def num_delayed_groups(self, num_delayed_groups): + check_type('num_delayed_groups', num_delayed_groups, int) + check_greater_than('num_delayed_groups', num_delayed_groups, 0, + equality=True) + check_less_than('num_delayed_groups', num_delayed_groups, + openmc.mgxs.MAX_DELAYED_GROUPS, equality=True) + self._num_delayed_groups = num_delayed_groups + def add_xsdata(self, xsdata): """Add an XSdata entry to the file. @@ -1319,10 +1889,12 @@ class MGXSLibrary(object): MGXS information to add """ + if not isinstance(xsdata, XSdata): msg = 'Unable to add a non-XSdata "{0}" to the ' \ 'MGXSLibrary instance'.format(xsdata) raise ValueError(msg) + if xsdata.energy_groups != self._energy_groups: msg = 'Energy groups of XSdata do not match that of MGXSLibrary.' raise ValueError(msg) @@ -1338,6 +1910,7 @@ class MGXSLibrary(object): XSdatas to add """ + check_iterable_type('xsdatas', xsdatas, XSdata) for xsdata in xsdatas: @@ -1374,7 +1947,8 @@ class MGXSLibrary(object): # Create and write to the HDF5 file file = h5py.File(filename, "w") - file.attrs['groups'] = self.energy_groups.num_groups + file.attrs['energy_groups'] = self.energy_groups.num_groups + file.attrs['delayed_groups'] = self.num_delayed_groups file.attrs['group structure'] = self.energy_groups.group_edges for xsdata in self._xsdatas: diff --git a/scripts/openmc-update-mgxs b/scripts/openmc-update-mgxs index 658be9fef7..cad48d0e31 100755 --- a/scripts/openmc-update-mgxs +++ b/scripts/openmc-update-mgxs @@ -85,13 +85,13 @@ if __name__ == '__main__': temp = np.array(temp.split()) group_structure = temp.astype(np.float) energy_groups = openmc.mgxs.EnergyGroups(group_structure) - temp = tree.find('inverse_velocities') + temp = tree.find('inverse-velocity') if temp is not None: temp = temp.text.strip() temp = np.array(temp.split()) - inverse_velocities = temp.astype(np.float) + inverse_velocity = temp.astype(np.float) else: - inverse_velocities = None + inverse_velocity = None xsd = [] names = [] @@ -166,8 +166,8 @@ if __name__ == '__main__': total.shape = xsd[i].vector_shape xsd[i].set_total(total, temperature) - if inverse_velocities is not None: - xsd[i].set_inverse_velocities(inverse_velocities, temperature) + if inverse_velocity is not None: + xsd[i].set_inverse_velocity(inverse_velocity, temperature) temp = get_data(xsdata_elem, 'absorption') temp = np.array(temp.split()) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 76aee23695..f107627503 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -107,7 +107,7 @@ contains if (.not. run_CE) then do i = 1, ng found = .false. - do g = 1, energy_groups + 1 + do g = 1, num_energy_groups + 1 if (cmfd % egrid(i) == energy_bins(g)) then found = .true. exit diff --git a/src/global.F90 b/src/global.F90 index 09a8422a11..8111fd5144 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -112,7 +112,10 @@ module global type(MgxsContainer), target, allocatable :: macro_xs(:) ! Number of energy groups - integer :: energy_groups + integer :: num_energy_groups + + ! Number of delayed groups + integer :: num_delayed_groups ! Energy group structure real(8), allocatable :: energy_bins(:) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index ce689afa3b..c8892c53c4 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -98,6 +98,7 @@ module hdf5_interface public :: open_dataset public :: close_dataset public :: get_shape + public :: get_ndims public :: write_attribute_string public :: get_groups public :: get_datasets @@ -2586,6 +2587,26 @@ contains end if end subroutine get_shape + subroutine get_ndims(obj_id, ndims) + integer(HID_T), intent(in) :: obj_id + integer(HID_T), intent(out) :: ndims + + integer :: hdf5_err + integer :: type + integer(HID_T) :: space_id + + call h5iget_type_f(obj_id, type, hdf5_err) + if (type == H5I_DATASET_F) then + call h5dget_space_f(obj_id, space_id, hdf5_err) + call h5sget_simple_extent_ndims_f(space_id, ndims, hdf5_err) + call h5sclose_f(space_id, hdf5_err) + elseif (type == H5I_ATTR_F) then + call h5aget_space_f(obj_id, space_id, hdf5_err) + call h5sget_simple_extent_ndims_f(space_id, ndims, hdf5_err) + call h5sclose_f(space_id, hdf5_err) + end if + end subroutine get_ndims + function using_mpio_device(obj_id) result(mpio) integer(HID_T), intent(in) :: obj_id logical :: mpio diff --git a/src/input_xml.F90 b/src/input_xml.F90 index aac22b9e88..ea798b9e98 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1062,19 +1062,6 @@ contains end select end if - ! Check whether create fission sites - if (run_mode == MODE_FIXEDSOURCE) then - if (check_for_node(doc, "create_fission_neutrons")) then - call get_node_value(doc, "create_fission_neutrons", temp_str) - temp_str = to_lower(temp_str) - if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') then - create_fission_neutrons = .true. - else if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') then - create_fission_neutrons = .false. - end if - end if - end if - ! Check for tabular_legendre options if (check_for_node(doc, "tabular_legendre")) then @@ -1100,6 +1087,19 @@ contains end if end if + ! Check whether create fission sites + if (run_mode == MODE_FIXEDSOURCE) then + if (check_for_node(doc, "create_fission_neutrons")) then + call get_node_value(doc, "create_fission_neutrons", temp_str) + temp_str = to_lower(temp_str) + if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') then + create_fission_neutrons = .true. + else if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') then + create_fission_neutrons = .false. + end if + end if + end if + ! Close settings XML file call close_xmldoc(doc) @@ -3017,8 +3017,8 @@ contains ! index is simply the group (after flipping for the different ! ordering of the library and tallying systems). if (.not. run_CE) then - if (n_words == energy_groups + 1) then - if (all(filt % bins == energy_bins(energy_groups + 1:1:-1))) & + if (n_words == num_energy_groups + 1) then + if (all(filt % bins == energy_bins(num_energy_groups + 1:1:-1))) & then filt % matches_transport_groups = .true. end if @@ -3043,8 +3043,8 @@ contains ! index is simply the group (after flipping for the different ! ordering of the library and tallying systems). if (.not. run_CE) then - if (n_words == energy_groups + 1) then - if (all(filt % bins == energy_bins(energy_groups + 1:1:-1))) & + if (n_words == num_energy_groups + 1) then + if (all(filt % bins == energy_bins(num_energy_groups + 1:1:-1))) & then filt % matches_transport_groups = .true. end if @@ -3058,14 +3058,6 @@ contains t % estimator = ESTIMATOR_ANALOG case ('delayedgroup') - ! Check to see if running in MG mode, because if so, the current - ! system isnt set up yet to support delayed group data and thus - ! these tallies - if (.not. run_CE) then - call fatal_error("delayedgroup filter on tally " & - // trim(to_str(t % id)) // " not yet supported& - & for multi-group mode.") - end if ! Allocate and declare the filter type allocate(DelayedGroupFilter::t % filters(j) % obj) @@ -3556,7 +3548,13 @@ contains end if case ('decay-rate') t % score_bins(j) = SCORE_DECAY_RATE - t % estimator = ESTIMATOR_ANALOG + + ! Set tally estimator to analog for CE mode + ! (MG mode has all data available without a collision being + ! necessary) + if (run_CE) then + t % estimator = ESTIMATOR_ANALOG + end if case ('delayed-nu-fission') t % score_bins(j) = SCORE_DELAYED_NU_FISSION if (t % find_filter(FILTER_ENERGYOUT) > 0) then @@ -3569,12 +3567,6 @@ contains ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG end if - - ! Disallow for MG mode since data not present - if (.not. run_CE) then - call fatal_error("Cannot tally delayed nu-fission rate in & - &multi-group mode") - end if case ('kappa-fission') t % score_bins(j) = SCORE_KAPPA_FISSION case ('inverse-velocity') @@ -4611,15 +4603,23 @@ contains ! Open file for reading file_id = file_open(path_cross_sections, 'r', parallel=.true.) - if (attribute_exists(file_id, "groups")) then - ! Get neutron group count - call read_attribute(energy_groups, file_id, "groups") + if (attribute_exists(file_id, "energy_groups")) then + ! Get neutron energy group count + call read_attribute(num_energy_groups, file_id, "energy_groups") else - call fatal_error("'groups' attribute must exist!") + call fatal_error("'energy_groups' attribute must exist!") end if - allocate(rev_energy_bins(energy_groups + 1)) - allocate(energy_bins(energy_groups + 1)) + if (attribute_exists(file_id, "delayed_groups")) then + ! Get neutron delayed group count + call read_attribute(num_delayed_groups, file_id, "delayed_groups") + else + num_delayed_groups = 0 + end if + + allocate(rev_energy_bins(num_energy_groups + 1)) + allocate(energy_bins(num_energy_groups + 1)) + if (attribute_exists(file_id, "group structure")) then ! Get neutron group structure call read_attribute(energy_bins, file_id, "group structure") @@ -4628,10 +4628,10 @@ contains end if ! First reverse the order of energy_groups - energy_bins = energy_bins(energy_groups + 1:1:-1) + energy_bins = energy_bins(num_energy_groups + 1:1:-1) - allocate(energy_bin_avg(energy_groups)) - do i = 1, energy_groups + allocate(energy_bin_avg(num_energy_groups)) + do i = 1, num_energy_groups energy_bin_avg(i) = HALF * (energy_bins(i) + energy_bins(i + 1)) end do diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index 0535adf6c9..f823304c0f 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -22,26 +22,25 @@ contains !=============================================================================== subroutine read_mgxs() - integer :: i ! index in materials array - integer :: j ! index over nuclides in material - integer :: i_xsdata ! index in xsdata_dict - integer :: i_nuclide ! index in nuclides array - character(20) :: name ! name of library to load - integer :: representation ! Data representation + integer :: i ! index in materials array + integer :: j ! index over nuclides in material + integer :: i_xsdata ! index in list + integer :: i_nuclide ! index in nuclides array + character(20) :: name ! name of library to load + integer :: representation ! Data representation character(MAX_LINE_LEN) :: temp_str type(Material), pointer :: mat - type(SetChar) :: already_read - integer(HID_T) :: file_id - integer(HID_T) :: xsdata_group - logical :: file_exists - logical :: get_kfiss, get_fiss - integer :: l - type(DictCharInt) :: xsdata_dict + type(SetChar) :: already_read + integer(HID_T) :: file_id + integer(HID_T) :: xsdata_group + logical :: file_exists + type(DictCharInt) :: xsdata_dict type(VectorReal), allocatable :: temps(:) ! Check if MGXS Library exists inquire(FILE=path_cross_sections, EXIST=file_exists) if (.not. file_exists) then + ! Could not find MGXS Library file call fatal_error("Cross sections HDF5 file '" & &// trim(path_cross_sections) // "' does not exist!") @@ -62,23 +61,6 @@ contains allocate(micro_xs(n_nuclides_total)) !$omp end parallel - ! Find out if we need fission & kappa fission - ! (i.e., are there any SCORE_FISSION or SCORE_KAPPA_FISSION tallies?) - get_kfiss = .false. - get_fiss = .false. - do i = 1, n_tallies - do l = 1, tallies(i) % n_score_bins - if (tallies(i) % score_bins(l) == SCORE_KAPPA_FISSION) then - get_kfiss = .true. - end if - if (tallies(i) % score_bins(l) == SCORE_FISSION .or. & - tallies(i) % score_bins(l) == SCORE_NU_FISSION) then - get_fiss = .true. - end if - end do - if (get_kfiss .and. get_fiss) exit - end do - ! ========================================================================== ! READ ALL MGXS CROSS SECTION TABLES @@ -105,7 +87,9 @@ contains ! First find out the data representation if (attribute_exists(xsdata_group, "representation")) then + call read_attribute(temp_str, xsdata_group, "representation") + if (trim(temp_str) == 'isotropic') then representation = MGXS_ISOTROPIC else if (trim(temp_str) == 'angle') then @@ -120,20 +104,24 @@ contains ! Now allocate accordingly select case(representation) + case(MGXS_ISOTROPIC) allocate(MgxsIso :: nuclides_MG(i_nuclide) % obj) + case(MGXS_ANGLE) allocate(MgxsAngle :: nuclides_MG(i_nuclide) % obj) + end select ! Now read in the data specific to the type we just declared call nuclides_MG(i_nuclide) % obj % from_hdf5(xsdata_group, & - energy_groups, temps(i_nuclide), temperature_method, & - temperature_tolerance, get_kfiss, get_fiss, max_order, & + num_energy_groups, num_delayed_groups, temps(i_nuclide), & + temperature_method, temperature_tolerance, max_order, & legendre_to_tabular, legendre_to_tabular_points) ! Add name to dictionary call already_read % add(name) + end if end do NUCLIDE_LOOP end do MATERIAL_LOOP @@ -149,6 +137,7 @@ contains ! Loop around nuclides in material NUCLIDE_LOOP2: do j = 1, mat % n_nuclides + ! Is this fissionable? if (nuclides_MG(mat % nuclide(j)) % obj % fissionable) then mat % fissionable = .true. @@ -163,12 +152,12 @@ contains end subroutine read_mgxs !=============================================================================== -! CREATE_MACRO_XS generates the macroscopic x/s from the microscopic input data +! CREATE_MACRO_XS generates the macroscopic xs from the microscopic input data !=============================================================================== subroutine create_macro_xs() - integer :: i_mat ! index in materials array - type(Material), pointer :: mat ! current material + integer :: i_mat ! index in materials array + type(Material), pointer :: mat ! current material type(VectorReal), allocatable :: kTs(:) allocate(macro_xs(n_materials)) @@ -176,27 +165,32 @@ contains ! Get temperatures to read for each material call get_mat_kTs(kTs) + ! Force all nuclides in a material to be the same representation. + ! Therefore type(nuclides(mat % nuclide(1)) % obj) dictates type(macroxs). + ! At the same time, we will find the scattering type, as that will dictate + ! how we allocate the scatter object within macroxs.allocate(macro_xs(n_materials)) do i_mat = 1, n_materials + + ! Get the material mat => materials(i_mat) - ! Check to see how our nuclides are represented - ! Force all to be the same type - ! Therefore type(nuclides(mat % nuclide(1)) % obj) dictates type(macroxs) + ! Get the scattering type for the first nuclide select type(nuc => nuclides_MG(mat % nuclide(1)) % obj) type is (MgxsIso) allocate(MgxsIso :: macro_xs(i_mat) % obj) type is (MgxsAngle) allocate(MgxsAngle :: macro_xs(i_mat) % obj) end select + ! Do not read materials which we do not actually use in the problem to ! save space if (allocated(kTs(i_mat) % data)) then call macro_xs(i_mat) % obj % combine(kTs(i_mat), mat, nuclides_MG, & - energy_groups, max_order, & - temperature_tolerance, & - temperature_method) + num_energy_groups, num_delayed_groups, max_order, & + temperature_tolerance, temperature_method) end if end do + end subroutine create_macro_xs !=============================================================================== @@ -205,16 +199,17 @@ contains !=============================================================================== subroutine get_mat_kTs(kTs) - type(VectorReal), allocatable, intent(out) :: kTs(:) - integer :: i, j + type(VectorReal), allocatable, intent(out) :: kTs(:) + integer :: i, j ! Cell and material index integer :: i_material ! Index in materials array - real(8) :: kT ! temperature in eV + real(8) :: kT ! temperature in eV allocate(kTs(size(materials))) do i = 1, size(cells) do j = 1, size(cells(i) % material) + ! Skip any non-material cells and void materials if (cells(i) % material(j) == NONE .or. & cells(i) % material(j) == MATERIAL_VOID) cycle @@ -238,5 +233,4 @@ contains end subroutine get_mat_kTs - end module mgxs_data diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index 98b2565c1f..fc5e8cb3e1 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -26,11 +26,14 @@ module mgxs_header real(8), allocatable :: total(:) ! total cross section real(8), allocatable :: absorption(:) ! absorption cross section class(ScattData), allocatable :: scatter ! scattering info - real(8), allocatable :: nu_fission(:) ! fission production - real(8), allocatable :: k_fission(:) ! kappa-fission - real(8), allocatable :: fission(:) ! fission - real(8), allocatable :: chi(:, :) ! Fission Spectra - real(8), allocatable :: inv_vel(:) ! Inverse velocities + real(8), allocatable :: delayed_nu_fission(:,:) ! Delayed fission matrix (Dg x Gin) + real(8), allocatable :: prompt_nu_fission(:) ! Prompt fission vector (Gin) + real(8), allocatable :: kappa_fission(:) ! Kappa fission + real(8), allocatable :: fission(:) ! Neutron production + real(8), allocatable :: decay_rate(:) ! Delayed neutron precursor decay rate + real(8), allocatable :: inverse_velocity(:) ! Inverse neutron velocity + real(8), allocatable :: chi_delayed(:, :, :) ! Delayed fission spectra + real(8), allocatable :: chi_prompt(:, :) ! Prompt fission spectra end type XsDataIso type :: XsDataAngle @@ -39,11 +42,14 @@ module mgxs_header real(8), allocatable :: total(:, :, :) ! total cross section real(8), allocatable :: absorption(:, :, :) ! absorption cross section type(ScattDataContainer), allocatable :: scatter(:, :) ! scattering info - real(8), allocatable :: nu_fission(:, :, :) ! fission production - real(8), allocatable :: k_fission(:, :, :) ! kappa-fission - real(8), allocatable :: fission(:, :, :) ! fission - real(8), allocatable :: chi(:, :, :, :) ! Fission Spectra - real(8), allocatable :: inv_vel(:, :, :) ! Inverse velocities + real(8), allocatable :: delayed_nu_fission(:, :, :, :) ! Delayed fission matrix (Gout x Gin) + real(8), allocatable :: prompt_nu_fission(:, :, :) ! Prompt fission matrix (Gout x Gin) + real(8), allocatable :: kappa_fission(:, :, :) ! Kappa fission + real(8), allocatable :: fission(:, :, :) ! Neutron production + real(8), allocatable :: decay_rate(:, :, :) ! Delayed neutron precursor decay rate + real(8), allocatable :: inverse_velocity(:, :, :) ! Inverse neutron velocity + real(8), allocatable :: chi_delayed(:, :, :, :, :) ! Delayed fission spectra + real(8), allocatable :: chi_prompt(:, :, :, :) ! Prompt fission spectra end type XsDataAngle !=============================================================================== @@ -58,6 +64,7 @@ module mgxs_header ! Fission information logical :: fissionable ! mgxs object is fissionable? integer :: scatter_format ! either legendre, histogram, or tabular. + integer :: num_delayed_groups ! Num delayed groups ! Caching information integer :: index_temp ! temperature index for nuclide @@ -66,12 +73,17 @@ module mgxs_header procedure(mgxs_from_hdf5_), deferred :: from_hdf5 ! Load the data procedure(mgxs_combine_), deferred :: combine ! initializes object procedure(mgxs_get_xs_), deferred :: get_xs ! Get the requested xs + ! Sample the outgoing energy from a fission event procedure(mgxs_sample_fission_), deferred :: sample_fission_energy + ! Sample the outgoing energy and angle from a scatter event procedure(mgxs_sample_scatter_), deferred :: sample_scatter + ! Calculate the material specific MGXS data from the nuclides procedure(mgxs_calculate_xs_), deferred :: calculate_xs + + ! Find the temperature procedure :: find_temperature => mgxs_find_temperature end type Mgxs @@ -88,38 +100,38 @@ module mgxs_header !=============================================================================== abstract interface - subroutine mgxs_from_hdf5_(this, xs_id, groups, temperature, method, & - tolerance, get_kfiss, get_fiss, max_order, & - legendre_to_tabular, legendre_to_tabular_points) + subroutine mgxs_from_hdf5_(this, xs_id, energy_groups, delayed_groups, & + temperature, method, tolerance, max_order, legendre_to_tabular, & + legendre_to_tabular_points) import Mgxs, HID_T, VectorReal class(Mgxs), intent(inout) :: this ! Working Object integer(HID_T), intent(in) :: xs_id ! Library data - integer, intent(in) :: groups ! Number of Energy groups + integer, intent(in) :: energy_groups ! Number of energy groups + integer, intent(in) :: delayed_groups ! Number of delayed groups type(VectorReal), intent(in) :: temperature ! list of desired temperatures integer, intent(inout) :: method ! Type of temperature access real(8), intent(in) :: tolerance ! Tolerance on method - logical, intent(in) :: get_kfiss ! Need Kappa-Fission? - logical, intent(in) :: get_fiss ! Should we get fiss data? integer, intent(in) :: max_order ! Maximum requested order logical, intent(in) :: legendre_to_tabular ! Convert Legendres to Tabular? integer, intent(in) :: legendre_to_tabular_points ! Number of points to use ! in that conversion end subroutine mgxs_from_hdf5_ - subroutine mgxs_combine_(this, temps, mat, nuclides, groups, max_order, & - tolerance, method) + subroutine mgxs_combine_(this, temps, mat, nuclides, energy_groups, & + delayed_groups, max_order, tolerance, method) import Mgxs, Material, MgxsContainer, VectorReal class(Mgxs), intent(inout) :: this ! The Mgxs to initialize type(VectorReal), intent(in) :: temps ! Temperatures to obtain type(Material), pointer, intent(in) :: mat ! base material type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from - integer, intent(in) :: groups ! Number of E groups + integer, intent(in) :: energy_groups ! Number of energy groups + integer, intent(in) :: delayed_groups ! Number of delayed groups integer, intent(in) :: max_order ! Maximum requested order real(8), intent(in) :: tolerance ! Tolerance on method integer, intent(in) :: method ! Type of temperature access end subroutine mgxs_combine_ - pure function mgxs_get_xs_(this, xstype, gin, gout, uvw, mu) result(xs_val) + pure function mgxs_get_xs_(this, xstype, gin, gout, uvw, mu, dg) result(xs_val) import Mgxs class(Mgxs), intent(in) :: this character(*), intent(in) :: xstype ! Cross Section Type @@ -127,17 +139,19 @@ module mgxs_header integer, optional, intent(in) :: gout ! Outgoing Group real(8), optional, intent(in) :: uvw(3) ! Requested Angle real(8), optional, intent(in) :: mu ! Change in angle + integer, optional, intent(in) :: dg ! Delayed group real(8) :: xs_val ! Resultant xs end function mgxs_get_xs_ - function mgxs_sample_fission_(this, gin, uvw) result(gout) + subroutine mgxs_sample_fission_(this, gin, uvw, dg, gout) import Mgxs class(Mgxs), intent(in) :: this integer, intent(in) :: gin ! Incoming energy group real(8), intent(in) :: uvw(3) ! Particle Direction - integer :: gout ! Sampled outgoing group + integer, intent(out) :: dg ! Delayed group + integer, intent(out) :: gout ! Sampled outgoing group - end function mgxs_sample_fission_ + end subroutine mgxs_sample_fission_ subroutine mgxs_sample_scatter_(this, uvw, gin, gout, mu, wgt) import Mgxs @@ -400,29 +414,30 @@ module mgxs_header end subroutine mgxs_from_hdf5 - subroutine mgxsiso_from_hdf5(this, xs_id, groups, temperature, method, & - tolerance, get_kfiss, get_fiss, max_order, & - legendre_to_tabular, legendre_to_tabular_points) + subroutine mgxsiso_from_hdf5(this, xs_id, energy_groups, delayed_groups, & + temperature, method, tolerance, max_order, & + legendre_to_tabular, legendre_to_tabular_points) class(MgxsIso), intent(inout) :: this ! Working Object integer(HID_T), intent(in) :: xs_id ! Group in H5 file - integer, intent(in) :: groups ! Number of Energy groups + integer, intent(in) :: energy_groups ! Number of energy groups + integer, intent(in) :: delayed_groups ! Number of delayed groups type(VectorReal), intent(in) :: temperature ! list of desired temperatures integer, intent(inout) :: method ! Type of temperature access real(8), intent(in) :: tolerance ! Tolerance on method - logical, intent(in) :: get_kfiss ! Need Kappa-Fission? - logical, intent(in) :: get_fiss ! Should we get fiss data? integer, intent(in) :: max_order ! Maximum requested order logical, intent(in) :: legendre_to_tabular ! Convert Legendres to Tabular? integer, intent(in) :: legendre_to_tabular_points ! Number of points to use ! in that conversion character(MAX_LINE_LEN) :: temp_str - integer(HID_T) :: xsdata_grp, scatt_grp + integer(HID_T) :: xsdata, xsdata_grp, scatt_grp, ndims + integer(HSIZE_T) :: dims(2) real(8), allocatable :: temp_arr(:), temp_2d(:, :) - real(8) :: dmu, mu, norm + real(8), allocatable :: temp_beta(:, :) + real(8) :: dmu, mu, norm, chi_sum integer :: order, order_dim, gin, gout, l, imu, length type(VectorInt) :: temps_to_read - integer :: t + integer :: t, dg type(Jagged2D), allocatable :: input_scatt(:), scatt_coeffs(:) type(Jagged1D), allocatable :: temp_mult(:) integer, allocatable :: gmin(:), gmax(:) @@ -431,105 +446,467 @@ module mgxs_header call mgxs_from_hdf5(this, xs_id, temperature, method, tolerance, & temps_to_read, order_dim) + ! Set the number of delayed groups + this % num_delayed_groups = delayed_groups + ! Load the more specific data do t = 1, temps_to_read % size() associate(xs => this % xs(t)) + ! Get temperature as a string temp_str = trim(to_str(temps_to_read % data(t))) // "K" xsdata_grp = open_group(xs_id, trim(temp_str)) - allocate(xs % nu_fission(groups)) - allocate(xs % chi(groups, groups)) + + ! Allocate data for all the cross sections + allocate(xs % total(energy_groups)) + allocate(xs % absorption(energy_groups)) + allocate(xs % delayed_nu_fission(delayed_groups, energy_groups)) + allocate(xs % prompt_nu_fission(energy_groups)) + allocate(xs % fission(energy_groups)) + allocate(xs % kappa_fission(energy_groups)) + allocate(xs % decay_rate(delayed_groups)) + allocate(xs % inverse_velocity(energy_groups)) + allocate(xs % chi_delayed(delayed_groups, energy_groups, & + energy_groups)) + allocate(xs % chi_prompt(energy_groups, energy_groups)) + + ! Set all fissionable terms to zero + xs % delayed_nu_fission = ZERO + xs % prompt_nu_fission = ZERO + xs % fission = ZERO + xs % kappa_fission = ZERO + xs % chi_delayed = ZERO + xs % chi_prompt = ZERO + xs % decay_rate = ZERO + xs % inverse_velocity = ZERO + if (this % fissionable) then - if (object_exists(xsdata_grp, "chi")) then - ! Chi was provided, that means we need chi and nu-fission vectors - ! Get chi - allocate(temp_arr(groups)) - call read_dataset(temp_arr, xsdata_grp, "chi") - do gin = 1, groups - do gout = 1, groups - xs % chi(gout, gin) = temp_arr(gout) + + ! Allocate temporary array for beta + allocate(temp_beta(delayed_groups, energy_groups)) + + ! Set beta + if (object_exists(xsdata_grp, "beta")) then + + ! Get the dimensions of the beta dataset + xsdata = open_dataset(xsdata_grp, "beta") + call get_ndims(xsdata, ndims) + + ! Beta is input as (delayed_groups) + if (ndims == 1) then + + ! Allocate temporary array for beta + allocate(temp_arr(delayed_groups)) + + ! Read beta + call read_dataset(temp_arr, xsdata_grp, "beta") + + do dg = 1, delayed_groups + do gin = 1, energy_groups + temp_beta(dg, gin) = temp_arr(dg) + end do end do - ! Normalize chi so its CDF goes to 1 - xs % chi(:, gin) = xs % chi(:, gin) / sum(xs % chi(:, gin)) - end do - deallocate(temp_arr) - ! Get nu_fission (as a vector) - if (object_exists(xsdata_grp, "nu-fission")) then - call read_dataset(xs % nu_fission, xsdata_grp, "nu-fission") + ! Deallocate temporary beta array + deallocate(temp_arr) + + ! Beta is input as (delayed_groups, energy_groups) + else if (ndims == 2) then + + ! Allocate temporary array for beta + allocate(temp_arr(delayed_groups * energy_groups)) + + ! Read beta + call read_dataset(temp_arr, xsdata_grp, "beta") + + ! Reshape array and set to dedicated beta array + temp_beta = reshape(temp_arr, (/delayed_groups, energy_groups/)) + + ! Deallocate temporary beta array + deallocate(temp_arr) + else - call fatal_error("If fissionable, must provide nu-fission!") + call fatal_error("beta must be provided as a 1D or 2D array") end if - else - ! chi isnt provided but is within nu_fission, existing as a matrix - ! So, get nu_fission (as a matrix) - if (object_exists(xsdata_grp, "nu-fission")) then - allocate(temp_2d(groups, groups)) - call read_dataset(temp_2d, xsdata_grp, "nu-fission") - else - call fatal_error("If fissionable, must provide nu-fission!") - end if + temp_beta = ZERO + end if - ! Set the vector nu-fission from the matrix nu-fission - do gin = 1, groups - xs % nu_fission(gin) = sum(temp_2d(:, gin)) + ! If chi provided, set chi-prompt and chi-delayed + if (object_exists(xsdata_grp, "chi")) then + + ! Allocate temporary array for chi + allocate(temp_arr(energy_groups)) + + ! Read chi + call read_dataset(temp_arr, xsdata_grp, "chi") + + do gin = 1, energy_groups + do gout = 1, energy_groups + xs % chi_prompt(gout, gin) = temp_arr(gout) + end do + + ! Normalize chi-prompt so its CDF goes to 1 + chi_sum =sum(xs % chi_prompt(:, gin)) + if (chi_sum == ZERO) then + call fatal_error("Encountered chi for a group that sums to & + &zero") + else + xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / chi_sum + end if end do - ! Now pull out information needed for chi - xs % chi(:, :) = temp_2d - ! Normalize chi so its CDF goes to 1 - do gin = 1, groups - xs % chi(:, gin) = xs % chi(:, gin) / sum(xs % chi(:, gin)) + ! Set chi-delayed to chi-prompt + do dg = 1, delayed_groups + xs % chi_delayed(dg, :, :) = xs % chi_prompt(:, :) end do - deallocate(temp_2d) + + ! Deallocate temporary chi array + deallocate(temp_arr) end if - ! If we have a need* for the fission and kappa-fission x/s, get them - ! (*Need is defined as will be using it to tally) - if (get_fiss) then - if (object_exists(xsdata_grp, "fission")) then - allocate(xs % fission(groups)) - call read_dataset(xs % fission, xsdata_grp, "fission") + ! If nu-fission provided, set prompt-nu_-ission and + ! delayed-nu-fission. If nu fission is a matrix, set chi-prompt and + ! chi-delayed. + if (object_exists(xsdata_grp, "nu-fission")) then + + ! Get the dimensions of the nu-fission dataset + xsdata = open_dataset(xsdata_grp, "nu-fission") + call get_ndims(xsdata, ndims) + + ! If nu-fission is a vector + if (ndims == 1) then + + ! Get nu-fission + call read_dataset(xs % prompt_nu_fission, xsdata_grp, & + "nu-fission") + + ! Set delayed-nu-fission and correct prompt-nu-fission with + ! beta + do gin = 1, energy_groups + do dg = 1, delayed_groups + + ! Set delayed-nu-fission using delayed neutron fraction + xs % delayed_nu_fission(dg, gin) = temp_beta(dg, gin) * & + xs % prompt_nu_fission(gin) + end do + + ! Correct prompt-nu-fission using delayed neutron fraction + if (delayed_groups > 0) then + xs % prompt_nu_fission(gin) = (1 - sum(temp_beta(:, gin))) & + * xs % prompt_nu_fission(gin) + end if + end do + + ! If nu-fission is a matrix, set prompt-nu-fission, + ! delayed-nu-fission, chi-prompt, and chi-delayed. + else if (ndims == 2) then + + ! chi is embedded in nu-fission -> extract chi + allocate(temp_arr(energy_groups * energy_groups)) + call read_dataset(temp_arr, xsdata_grp, "nu-fission") + allocate(temp_2d(energy_groups, energy_groups)) + temp_2d = reshape(temp_arr, (/energy_groups, energy_groups/)) + + ! Deallocate temporary 1D array for nu-fission matrix + deallocate(temp_arr) + + ! Set the vector nu-fission from the matrix nu-fission + do gin = 1, energy_groups + xs % prompt_nu_fission(gin) = sum(temp_2d(:, gin)) + end do + + ! Set delayed-nu-fission and correct prompt-nu-fission with + ! beta + do gin = 1, energy_groups + do dg = 1, delayed_groups + + ! Set delayed-nu-fission using delayed neutron fraction + xs % delayed_nu_fission(dg, gin) = temp_beta(dg, gin) * & + xs % prompt_nu_fission(gin) + end do + + ! Correct prompt-nu-fission using delayed neutron fraction + if (delayed_groups > 0) then + xs % prompt_nu_fission(gin) = (1 - sum(temp_beta(:, gin))) & + * xs % prompt_nu_fission(gin) + end if + 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 for a group that sums to & + &zero") + else + xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / chi_sum + end if + end do + + ! Set chi-delayed to chi-prompt + do dg = 1, delayed_groups + xs % chi_delayed(dg, :, :) = xs % chi_prompt(:, :) + end do else - call fatal_error("Fission data missing, required due to fission& - & tallies in tallies.xml file!") + call fatal_error("nu-fission must be provided as a 1D or 2D & + &array") end if end if - if (get_kfiss) then - if (object_exists(xsdata_grp, "kappa-fission")) then - allocate(xs % k_fission(groups)) - call read_dataset(xs % k_fission, xsdata_grp, "kappa-fission") + + ! If chi-prompt provided, set chi-prompt + if (object_exists(xsdata_grp, "chi-prompt")) then + + ! Allocate temporary array for chi-prompt + allocate(temp_arr(energy_groups)) + + ! Get array with chi-prompt + call read_dataset(temp_arr, xsdata_grp, "chi-prompt") + + do gin = 1, energy_groups + do gout = 1, energy_groups + xs % chi_prompt(gout, gin) = temp_arr(gout) + end do + + ! Normalize chi so its CDF goes to 1 + 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 + + ! Deallocate temporary array for chi-prompt + deallocate(temp_arr) + end if + + ! If chi-delayed provided, set chi-delayed + if (object_exists(xsdata_grp, "chi-delayed")) then + + ! Get the dimensions of the chi-delayed dataset + xsdata = open_dataset(xsdata_grp, "chi-delayed") + call get_ndims(xsdata, ndims) + + ! If chi-delayed is a vector + if (ndims == 1) then + + ! Allocate temporary array for chi-delayed + allocate(temp_arr(energy_groups)) + + ! Get chi-delayed + call read_dataset(temp_arr, xsdata_grp, "chi-delayed") + + do dg = 1, delayed_groups + do gin = 1, energy_groups + do gout = 1, energy_groups + xs % chi_delayed(dg, gout, gin) = temp_arr(gout) + end do + + ! Normalize chi so its CDF goes to 1 + 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 array for chi-delayed + deallocate(temp_arr) + + else if (ndims == 2) then + + ! Allocate temporary array for chi-delayed + allocate(temp_arr(delayed_groups * energy_groups)) + + ! Get chi-delayed + call read_dataset(temp_arr, xsdata_grp, "chi-delayed") + allocate(temp_2d(delayed_groups, energy_groups)) + temp_2d = reshape(temp_arr, (/delayed_groups, energy_groups/)) + + do dg = 1, delayed_groups + do gin = 1, energy_groups + do gout = 1, energy_groups + xs % chi_delayed(dg, gout, gin) = temp_2d(dg, gout) + end do + + ! Normalize chi so its CDF goes to 1 + 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 arrays for chi-delayed + deallocate(temp_arr) + deallocate(temp_2d) + else - call fatal_error("kappa-fission data missing, required due to & - &kappa-fission tallies in tallies.xml file!") + call fatal_error("chi-delayed must be provided as a 1D or 2D & + &array") end if end if - else - xs % nu_fission = ZERO - xs % chi = ZERO + + ! 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") + end if + + ! If delayed-nu-fission provided, set delayed-nu-fission. If + ! delayed-nu-fission is a matrix, set chi-delayed. + if (object_exists(xsdata_grp, "delayed-nu-fission")) then + + ! Get the dimensions of the delayed-nu-fission dataset + xsdata = open_dataset(xsdata_grp, "delayed-nu-fission") + call get_ndims(xsdata, ndims) + + ! If delayed-nu-fission is a vector + if (ndims == 1) then + + ! If beta is zeros, raise error + if (temp_beta(1,1) == ZERO) then + call fatal_error("cannot set delayed-nu-fission with a 1D & + &array if beta not provided") + end if + + ! Allocate temporary array for delayed-nu-fission + allocate(temp_arr(energy_groups)) + + ! Get delayed-nu-fission + call read_dataset(temp_arr, xsdata_grp, "delayed-nu-fission") + + do gin = 1, energy_groups + do dg = 1, delayed_groups + + ! Set delayed-nu-fission using delayed neutron fraction + xs % delayed_nu_fission(dg, gin) = temp_beta(dg, gin) * & + temp_arr(gin) + end do + end do + + ! Deallocate temporary delayed-nu-fission array + deallocate(temp_arr) + + ! If delayed-nu-fission is a (delayed_group, energy_group) + ! matrix, set delayed-nu-fission separately for each delayed + ! group. + else if (ndims == 2) then + + ! Get the shape of delayed-nu-fission + call get_shape(xsdata, dims) + + ! Issue error if 1st dimension not correct + if (dims(1) /= delayed_groups) then + call fatal_error("The delayed-nu-fission matrix was input & + &with a 1st dimension not equal to the number of & + &delayed groups.") + end if + + ! Issue error if 2nd dimension not correct + if (dims(2) /= energy_groups) then + call fatal_error("The delayed-nu-fission matrix was input & + &with a 2nd dimension not equal to the number of & + &energy groups.") + end if + + ! Issue warning if delayed_groups == energy_groups + if (delayed_groups == energy_groups) then + call warning("delayed-nu-fission was input as a dimension & + &2 matrix with the same number of delayed groups and & + &groups. It is important to know that OpenMC assumes & + &the dimensions in the matrix are (delayed_groups, & + &energy_groups). Currently, delayed-nu-fission cannot & + &be set as a group by group matrix.") + end if + + ! Get delayed-nu-fission + allocate(temp_arr(delayed_groups * energy_groups)) + call read_dataset(temp_arr, xsdata_grp, "delayed-nu-fission") + xs % delayed_nu_fission = reshape(temp_arr, (/delayed_groups, & + energy_groups/)) + + ! Deallocate temporary array for delayed-nu-fission matrix + deallocate(temp_arr) + + else + call fatal_error("delayed-nu-fission must be provided as a & + &1D or 2D array") + end if + end if + + ! Deallocate temporary beta array + deallocate(temp_beta) + + ! chi-prompt, chi-delayed, prompt-nu-fission, and delayed-nu-fission + ! have been set; Now we will check for the rest of the XS that are + ! unique to fissionable isotopes + + ! Get fission xs + if (object_exists(xsdata_grp, "fission")) then + call read_dataset(xs % fission, xsdata_grp, "fission") + end if + + ! Get kappa-fission xs + if (object_exists(xsdata_grp, "kappa-fission")) then + call read_dataset(xs % kappa_fission, xsdata_grp, "kappa-fission") + end if + + ! Get decay rate xs + if (object_exists(xsdata_grp, "decay rate")) then + call read_dataset(xs % decay_rate, xsdata_grp, "decay rate") + end if end if + ! All the XS unique to fissionable isotopes have been set; Now set all + ! the generation XS + if (object_exists(xsdata_grp, "absorption")) then - allocate(xs % absorption(groups)) call read_dataset(xs % absorption, xsdata_grp, "absorption") else call fatal_error("Must provide absorption!") end if + ! Get inverse velocity + if (object_exists(xsdata_grp, "inverse-velocity")) then + call read_dataset(xs % inverse_velocity, xsdata_grp, & + "inverse-velocity") + end if + ! Get scattering data if (.not. object_exists(xsdata_grp, "scatter_data")) & call fatal_error("Must provide 'scatter_data'") + scatt_grp = open_group(xsdata_grp, 'scatter_data') + ! First get the outgoing group boundary indices if (object_exists(scatt_grp, "g_min")) then - allocate(gmin(groups)) + allocate(gmin(energy_groups)) call read_dataset(gmin, scatt_grp, "g_min") else call fatal_error("'g_min' for the scatter_data must be provided") end if + if (object_exists(scatt_grp, "g_max")) then - allocate(gmax(groups)) + allocate(gmax(energy_groups)) call read_dataset(gmax, scatt_grp, "g_max") else call fatal_error("'g_max' for the scatter_data must be provided") @@ -538,11 +915,14 @@ module mgxs_header ! Now use this information to find the length of a container array ! to hold the flattened data length = 0 - do gin = 1, groups + + do gin = 1, energy_groups length = length + order_dim * (gmax(gin) - gmin(gin) + 1) end do + ! Allocate flattened array allocate(temp_arr(length)) + if (.not. object_exists(scatt_grp, 'scatter_matrix')) & call fatal_error("'scatter_matrix' must be provided") call read_dataset(temp_arr, scatt_grp, "scatter_matrix") @@ -554,11 +934,12 @@ module mgxs_header order_dim = order + 1 end if - ! Convert temp_arr to a jagged array ((gin) % data(l, gout)) for passing - ! to ScattData - allocate(input_scatt(groups)) + ! Convert temp_arr to a jagged array ((gin) % data(l, gout)) for + ! passing to ScattData + allocate(input_scatt(energy_groups)) + index = 1 - do gin = 1, groups + do gin = 1, energy_groups allocate(input_scatt(gin) % data(order_dim, gmin(gin):gmax(gin))) do gout = gmin(gin), gmax(gin) do l = 1, order_dim @@ -567,21 +948,28 @@ module mgxs_header end do end do end do + deallocate(temp_arr) ! Finally convert the legendre to tabular if needed - allocate(scatt_coeffs(groups)) + allocate(scatt_coeffs(energy_groups)) + if (this % scatter_format == ANGLE_LEGENDRE .and. & legendre_to_tabular) then + this % scatter_format = ANGLE_TABULAR order_dim = legendre_to_tabular_points order = order_dim dmu = TWO / real(order - 1, 8) - do gin = 1, groups + + do gin = 1, energy_groups allocate(scatt_coeffs(gin) % data(order_dim, gmin(gin):gmax(gin))) do gout = gmin(gin), gmax(gin) + norm = ZERO + do imu = 1, order_dim + if (imu == 1) then mu = -ONE else if (imu == order_dim) then @@ -589,11 +977,14 @@ module mgxs_header else mu = -ONE + real(imu - 1, 8) * dmu end if + scatt_coeffs(gin) % data(imu, gout) = & evaluate_legendre(input_scatt(gin) % data(:, gout), mu) + ! Ensure positivity of distribution if (scatt_coeffs(gin) % data(imu, gout) < ZERO) & scatt_coeffs(gin) % data(imu, gout) = ZERO + ! And accrue the integral if (imu > 1) then norm = norm + HALF * dmu * & @@ -601,8 +992,10 @@ module mgxs_header scatt_coeffs(gin) % data(imu, gout)) end if end do ! mu - ! Now that we have the integral, lets ensure that the distribution - ! is normalized such that it preserves the original scattering xs + + ! Now that we have the integral, lets ensure that the + ! distribution is normalized such that it preserves the original + ! scattering xs if (norm > ZERO) then scatt_coeffs(gin) % data(:, gout) = & scatt_coeffs(gin) % data(:, gout) * & @@ -611,33 +1004,41 @@ module mgxs_header end do ! gout end do ! gin else + ! Sticking with current representation - do gin = 1, groups + do gin = 1, energy_groups allocate(scatt_coeffs(gin) % data(order_dim, gmin(gin):gmax(gin))) scatt_coeffs(gin) % data(:, :) = & input_scatt(gin) % data(1:order_dim, :) end do end if + deallocate(input_scatt) ! Now get the multiplication matrix - if (object_exists(scatt_grp, 'multiplicity_matrix')) then + if (object_exists(scatt_grp, 'multiplicity matrix')) then + ! Now use this information to find the length of a container array ! to hold the flattened data length = 0 - do gin = 1, groups + + do gin = 1, energy_groups length = length + (gmax(gin) - gmin(gin) + 1) end do + ! Allocate flattened array allocate(temp_arr(length)) - call read_dataset(temp_arr, scatt_grp, "multiplicity_matrix") + call read_dataset(temp_arr, scatt_grp, "multiplicity matrix") + + ! Convert temp_arr to a jagged array ((gin) % data(gout)) for + ! passing to ScattData + allocate(temp_mult(energy_groups)) - ! Convert temp_arr to a jagged array ((gin) % data(gout)) for passing - ! to ScattData - allocate(temp_mult(groups)) index = 1 - do gin = 1, groups + do gin = 1, energy_groups + allocate(temp_mult(gin) % data(gmin(gin):gmax(gin))) + do gout = gmin(gin), gmax(gin) temp_mult(gin) % data(gout) = temp_arr(index) index = index + 1 @@ -645,9 +1046,11 @@ module mgxs_header end do deallocate(temp_arr) else + ! Default to multiplicities of 1.0 - allocate(temp_mult(groups)) - do gin = 1, groups + allocate(temp_mult(energy_groups)) + + do gin = 1, energy_groups allocate(temp_mult(gin) % data(gmin(gin):gmax(gin))) temp_mult(gin) % data = ONE end do @@ -668,12 +1071,11 @@ module mgxs_header ! Check sigA to ensure it is not 0 since it is ! often divided by in the tally routines ! (This may happen with Helium data) - do gin = 1, groups + do gin = 1, energy_groups if (xs % absorption(gin) == ZERO) xs % absorption(gin) = 1E-10_8 end do ! Get, or infer, total xs data. - allocate(xs % total(groups)) if (object_exists(xsdata_grp, "total")) then call read_dataset(xs % total, xsdata_grp, "total") else @@ -682,46 +1084,42 @@ module mgxs_header ! Check sigT to ensure it is not 0 since it is ! often divided by in the tally routines - do gin = 1, groups + do gin = 1, energy_groups if (xs % total(gin) == ZERO) xs % total(gin) = 1E-10_8 end do - ! Get kinetics data - if (object_exists(xsdata_grp, "inverse_velocities")) then - allocate(xs % inv_vel(groups)) - call read_dataset(xs % inv_vel, xsdata_grp, "inverse_velocities") - end if - ! Close the groups we have opened and deallocate call close_group(xsdata_grp) call close_group(scatt_grp) deallocate(scatt_coeffs, temp_mult) end associate ! xs end do ! Temperatures + end subroutine mgxsiso_from_hdf5 - subroutine mgxsang_from_hdf5(this, xs_id, groups, temperature, method, & - tolerance, get_kfiss, get_fiss, max_order, & - legendre_to_tabular, legendre_to_tabular_points) + subroutine mgxsang_from_hdf5(this, xs_id, energy_groups, delayed_groups, & + temperature, method, tolerance, max_order, legendre_to_tabular, & + legendre_to_tabular_points) class(MgxsAngle), intent(inout) :: this ! Working Object integer(HID_T), intent(in) :: xs_id ! Group in H5 file - integer, intent(in) :: groups ! Number of Energy groups + integer, intent(in) :: energy_groups ! Number of energy groups + integer, intent(in) :: delayed_groups ! Number of energy groups type(VectorReal), intent(in) :: temperature ! list of desired temperatures integer, intent(inout) :: method ! Type of temperature access real(8), intent(in) :: tolerance ! Tolerance on method - logical, intent(in) :: get_kfiss ! Need Kappa-Fission? - logical, intent(in) :: get_fiss ! Should we get fiss data? integer, intent(in) :: max_order ! Maximum requested order logical, intent(in) :: legendre_to_tabular ! Convert Legendres to Tabular? integer, intent(in) :: legendre_to_tabular_points ! Number of points to use ! in that conversion character(MAX_LINE_LEN) :: temp_str - integer(HID_T) :: xsdata_grp, scatt_grp + integer(HID_T) :: xsdata, xsdata_grp, scatt_grp, ndims + integer(HSIZE_T) :: dims(4) integer, allocatable :: int_arr(:) - real(8), allocatable :: temp_arr(:), temp_4d(:, :, :, :) - real(8) :: dmu, mu, norm - integer :: order, order_dim, gin, gout, l, imu + real(8), allocatable :: temp_1d(:), temp_3d(:, :, :) + real(8), allocatable :: temp_4d(:, :, :, :), temp_beta(:, :, :, :) + real(8) :: dmu, mu, norm, chi_sum + integer :: order, order_dim, gin, gout, l, imu, dg type(VectorInt) :: temps_to_read integer :: t, length, ipol, iazi type(Jagged2D), allocatable :: input_scatt(:, :, :), scatt_coeffs(:, :, :) @@ -732,160 +1130,637 @@ module mgxs_header call mgxs_from_hdf5(this, xs_id, temperature, method, tolerance, & temps_to_read, order_dim) + ! Set the number of delayed groups + this % num_delayed_groups = delayed_groups + ! Load the more specific data do t = 1, temps_to_read % size() associate(xs => this % xs(t)) + ! Get temperature as a string temp_str = trim(to_str(temps_to_read % data(t))) // "K" xsdata_grp = open_group(xs_id, trim(temp_str)) - allocate(xs % nu_fission(groups, this % n_azi, this % n_pol)) - allocate(xs % chi(groups, groups, this % n_azi, this % n_pol)) + + ! Load the more specific data + allocate(xs % prompt_nu_fission(energy_groups, this % n_azi, & + this % n_pol)) + allocate(xs % delayed_nu_fission(delayed_groups, energy_groups, & + this % n_azi, this % n_pol)) + allocate(xs % chi_prompt(energy_groups, energy_groups, this % n_azi, & + this % n_pol)) + allocate(xs % chi_delayed(delayed_groups, energy_groups, & + energy_groups, this % n_azi, this % n_pol)) + allocate(xs % total(energy_groups, this % n_azi, this % n_pol)) + allocate(xs % absorption(energy_groups, this % n_azi, this % n_pol)) + allocate(xs % fission(energy_groups, this % n_azi, this % n_pol)) + allocate(xs % kappa_fission(energy_groups, this % n_azi, & + this % n_pol)) + allocate(xs % decay_rate(delayed_groups, this % n_azi, this % n_pol)) + allocate(xs % inverse_velocity(energy_groups, this % n_azi, & + this % n_pol)) + + ! Set all fissionable terms to zero + xs % delayed_nu_fission = ZERO + xs % prompt_nu_fission = ZERO + xs % fission = ZERO + xs % kappa_fission = ZERO + xs % chi_delayed = ZERO + xs % chi_prompt = ZERO + xs % decay_rate = ZERO + xs % inverse_velocity = ZERO + if (this % fissionable) then - if (object_exists(xsdata_grp, "chi")) then - ! Chi was provided, that means we need chi and nu-fission vectors - ! Get chi - allocate(temp_arr(groups * this % n_azi * this % n_pol)) - call read_dataset(temp_arr, xsdata_grp, "chi") - ! Initialize counter for temp_arr - l = 0 - gin = 1 - do ipol = 1, this % n_pol - do iazi = 1, this % n_azi - do gout = 1, groups - l = l + 1 - xs % chi(gout, gin, iazi, ipol) = temp_arr(l) - end do - ! Normalize chi so its CDF goes to 1 - xs % chi(:, gin, iazi, ipol) = & - xs % chi(:, gin, iazi, ipol) / & - sum(xs % chi(:, gin, iazi, ipol)) - end do - end do - ! Now set all the other gin values - do ipol = 1, this % n_pol - do iazi = 1, this % n_azi - do gin = 2, groups - xs % chi(:, gin, iazi, ipol) = & - xs % chi(:, 1, iazi, ipol) + ! Allocate temporary array for beta + allocate(temp_beta(delayed_groups, energy_groups, this % n_azi, & + this % n_pol)) + + ! Set beta + if (object_exists(xsdata_grp, "beta")) then + + ! Get the dimensions of the beta dataset + xsdata = open_dataset(xsdata_grp, "beta") + call get_ndims(xsdata, ndims) + + ! Beta is input as (delayed_groups, n_azi, n_pol) + if (ndims == 3) then + + ! Allocate temporary arrays for beta + allocate(temp_1d(delayed_groups * this % n_azi * this % n_pol)) + allocate(temp_3d(delayed_groups, this % n_azi, this % n_pol)) + + ! Read beta + call read_dataset(temp_1d, xsdata_grp, "beta") + temp_3d = reshape(temp_1d, (/delayed_groups, this % n_azi, & + this % n_pol/)) + + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do dg = 1, delayed_groups + do gin = 1, energy_groups + temp_beta(dg, gin, iazi, ipol) = temp_3d(dg, iazi, ipol) + end do + end do end do end do - end do - deallocate(temp_arr) - ! Get nu_fission (as a vector) - if (object_exists(xsdata_grp, "nu-fission")) then - allocate(temp_arr(groups * this % n_azi * this % n_pol)) - call read_dataset(temp_arr, xsdata_grp, "nu-fission") - xs % nu_fission = reshape(temp_arr, (/groups, this % n_azi, & - this % n_pol/)) - deallocate(temp_arr) + ! Deallocate temporary beta arrays + deallocate(temp_1d) + deallocate(temp_3d) + + ! Beta is input as (delayed_groups, energy_groups, n_azi, n_pol) + else if (ndims == 4) then + + ! Allocate temporary array for beta + allocate(temp_1d(delayed_groups * energy_groups * this % n_azi & + * this % n_pol)) + + ! Read beta + call read_dataset(temp_1d, xsdata_grp, "beta") + + ! Reshape array and set to dedicated beta array + temp_beta = reshape(temp_1d, (/delayed_groups, & + energy_groups, this % n_azi, this % n_pol/)) + + ! Deallocate temporary beta array + deallocate(temp_1d) + else - call fatal_error("If fissionable, must provide nu-fission!") + call fatal_error("beta must be provided as a 3D or 4D array") end if - else - ! chi isnt provided but is within nu_fission, existing as a matrix - ! So, get nu_fission (as a matrix) - if (object_exists(xsdata_grp, "nu-fission")) then - allocate(temp_arr(groups * groups * this % n_azi * this % n_pol)) - call read_dataset(temp_arr, xsdata_grp, "nu-fission") - allocate(temp_4d(groups, groups, this % n_azi, this % n_pol)) - temp_4d(:, :, :, :) = reshape(temp_arr, (/groups, groups, & + temp_beta = ZERO + end if + + ! If chi provided, set chi-prompt and chi-delayed + if (object_exists(xsdata_grp, "chi")) then + + ! Allocate temporary array for chi + allocate(temp_1d(energy_groups * this % n_azi * this % n_pol)) + allocate(temp_3d(energy_groups, this % n_azi, this % n_pol)) + + ! Read chi + call read_dataset(temp_1d, xsdata_grp, "chi") + temp_3d = reshape(temp_1d, (/energy_groups, this % n_azi, & + this % n_pol/)) + + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gin = 1, energy_groups + do gout = 1, energy_groups + xs % chi_prompt(gout, gin, iazi, ipol) = & + temp_3d(gout, iazi, ipol) + end do + + ! Normalize chi-prompt so its CDF goes to 1 + chi_sum = sum(xs % chi_prompt(:, gin, iazi, ipol)) + if (chi_sum == ZERO) then + call fatal_error("Encountered chi 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 + + ! Set chi-delayed to chi-prompt + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do dg = 1, delayed_groups + xs % chi_delayed(dg, :, :, iazi, ipol) = & + xs % chi_prompt(:, :, iazi, ipol) + end do + end do + end do + + ! Deallocate temporary chi arrays + deallocate(temp_1d) + deallocate(temp_3d) + end if + + ! If nu-fission provided, set prompt-nu_-ission and + ! delayed-nu-fission. If nu fission is a matrix, set chi-prompt and + ! chi-delayed. + if (object_exists(xsdata_grp, "nu-fission")) then + + ! Get the dimensions of the nu-fission dataset + xsdata = open_dataset(xsdata_grp, "nu-fission") + call get_ndims(xsdata, ndims) + + ! If nu-fission is a 3D array + if (ndims == 3) then + + ! Get nu-fission + call read_dataset(xs % prompt_nu_fission, xsdata_grp, & + "nu-fission") + + ! Set delayed-nu-fission and correct prompt-nu-fission with + ! beta + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gin = 1, energy_groups + do dg = 1, delayed_groups + + ! Set delayed-nu-fission using delayed neutron fraction + xs % delayed_nu_fission(dg, gin, iazi, ipol) = & + temp_beta(dg, gin, iazi, ipol) * & + xs % prompt_nu_fission(gin, iazi, ipol) + end do + + ! Correct prompt-nu-fission using delayed neutron fraction + if (delayed_groups > 0) then + xs % prompt_nu_fission(gin, iazi, ipol) = & + (1 - sum(temp_beta(:, gin, iazi, ipol))) * & + xs % prompt_nu_fission(gin, iazi, ipol) + end if + end do + end do + end do + + ! If nu-fission is a matrix, set prompt-nu-fission, + ! delayed-nu-fission, chi-prompt, and chi-delayed. + else if (ndims == 4) then + + ! chi is embedded in nu-fission -> extract chi + allocate(temp_1d(energy_groups * energy_groups * & + this % n_azi * this % n_pol)) + call read_dataset(temp_1d, xsdata_grp, "nu-fission") + allocate(temp_4d(energy_groups, energy_groups, this % n_azi, & + this % n_pol)) + temp_4d = reshape(temp_1d, (/energy_groups, energy_groups, & + this % n_azi, this % n_pol /)) + + ! Deallocate temporary 1D array for nu-fission matrix + deallocate(temp_1d) + + ! Set the vector nu-fission from the matrix 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 + + ! Set delayed-nu-fission and correct prompt-nu-fission with + ! beta + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gin = 1, energy_groups + do dg = 1, delayed_groups + + ! Set delayed-nu-fission using delayed neutron fraction + xs % delayed_nu_fission(dg, gin, iazi, ipol) = & + temp_beta(dg, gin, iazi, ipol) * & + xs % prompt_nu_fission(gin, iazi, ipol) + end do + + ! Correct prompt-nu-fission using delayed neutron fraction + if (delayed_groups > 0) then + xs % prompt_nu_fission(gin, iazi, ipol) = & + (1 - sum(temp_beta(:, gin, iazi, ipol))) * & + xs % prompt_nu_fission(gin, iazi, ipol) + end if + 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 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 + + ! Set chi-delayed to chi-prompt + do dg = 1, delayed_groups + xs % chi_delayed(dg, :, :, iazi, ipol) = & + xs % chi_prompt(:, :, iazi, ipol) + end do + end do + end do + else + call fatal_error("nu-fission must be provided as a 3D or & + &4D array") + end if + end if + + ! If chi-prompt provided, set chi-prompt + if (object_exists(xsdata_grp, "chi-prompt")) then + + ! Allocate temporary array for chi-prompt + allocate(temp_1d(energy_groups * this % n_azi * this % n_pol)) + allocate(temp_3d(energy_groups, this % n_azi, this % n_pol)) + + ! Get array with chi-prompt + call read_dataset(temp_1d, xsdata_grp, "chi-prompt") + temp_3d = reshape(temp_1d, (/energy_groups, this % n_azi, & + this % n_pol/)) + + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gin = 1, energy_groups + do gout = 1, energy_groups + xs % chi_prompt(gout, gin, iazi, ipol) = & + temp_3d(gout, iazi, ipol) + end do + + ! Normalize chi so its CDF goes to 1 + 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 + + ! Deallocate temporary arrays for chi-prompt + deallocate(temp_1d) + deallocate(temp_3d) + end if + + ! If chi-delayed provided, set chi-delayed + if (object_exists(xsdata_grp, "chi-delayed")) then + + ! Get the dimensions of the chi-delayed dataset + xsdata = open_dataset(xsdata_grp, "chi-delayed") + call get_ndims(xsdata, ndims) + + ! chi-delayed is input as (energy_groups, n_azi, n_pol) + if (ndims == 3) then + + ! Allocate temporary array for chi-prompt + allocate(temp_1d(energy_groups * this % n_azi * this % n_pol)) + allocate(temp_3d(energy_groups, this % n_azi, this % n_pol)) + + ! Get array with chi-prompt + call read_dataset(temp_1d, xsdata_grp, "chi-delayed") + temp_3d = reshape(temp_1d, (/energy_groups, this % n_azi, & + this % n_pol/)) + + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do dg = 1, delayed_groups + do gin = 1, energy_groups + do gout = 1, energy_groups + xs % chi_delayed(dg, gout, gin, iazi, ipol) = & + temp_3d(gout, iazi, ipol) + end do + + ! Normalize chi so its CDF goes to 1 + 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 arrays for chi-delayed + deallocate(temp_1d) + deallocate(temp_3d) + + ! chi-delayed is input as (delayed_groups, energy_groups, n_azi, + ! n_pol) + else if (ndims == 4) then + + ! Allocate temporary array for chi-delayed + allocate(temp_1d(delayed_groups * energy_groups * this % n_azi & + * this % n_pol)) + allocate(temp_4d(delayed_groups, energy_groups, this % n_azi, & + this % n_pol)) + + ! Get chi-delayed + call read_dataset(temp_1d, xsdata_grp, "chi-delayed") + temp_4d = reshape(temp_1d, (/delayed_groups, energy_groups, & this % n_azi, this % n_pol/)) - deallocate(temp_arr) - else - call fatal_error("If fissionable, must provide nu-fission!") - end if - ! Set the vector nu-fission from the matrix nu-fission - do ipol = 1, this % n_pol - do iazi = 1, this % n_azi - do gin = 1, groups - xs % nu_fission(gin, iazi, ipol) = & - sum(temp_4d(:, gin, iazi, ipol)) + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do dg = 1, delayed_groups + do gin = 1, energy_groups + do gout = 1, energy_groups + xs % chi_delayed(dg, gout, gin, iazi, ipol) = & + temp_4d(dg, gout, iazi, ipol) + end do + + ! Normalize chi so its CDF goes to 1 + 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 - end do - ! Now pull out information needed for chi - xs % chi = 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, groups - xs % chi(:, gin, iazi, ipol) = & - xs % chi(:, gin, iazi, ipol) / & - sum(xs % chi(:, gin, iazi, ipol)) + ! Deallocate temporary arrays for chi-delayed + deallocate(temp_1d) + deallocate(temp_4d) + + else + call fatal_error("chi-delayed must be provided as a 3D or 4D & + &array") + end if + end if + + ! 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)) + + ! 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/)) + + ! Deallocate temporary array for prompt-nu-fission + deallocate(temp_1d) + end if + + ! If delayed-nu-fission provided, set delayed-nu-fission. If + ! delayed-nu-fission is a matrix, set chi-delayed. + if (object_exists(xsdata_grp, "delayed-nu-fission")) then + + ! Get the dimensions of the delayed-nu-fission dataset + xsdata = open_dataset(xsdata_grp, "delayed-nu-fission") + call get_ndims(xsdata, ndims) + + ! delayed-nu-fission is input as (energy_groups, n_azi, n_pol) + if (ndims == 3) then + + ! If beta is zeros, raise error + if (temp_beta(1,1,1,1) == ZERO) then + call fatal_error("cannot set delayed-nu-fission with a 3D & + &array if beta not provided") + end if + + ! Allocate temporary arrays for delayed-nu-fission + allocate(temp_1d(energy_groups * this % n_azi * this % n_pol)) + allocate(temp_3d(energy_groups, this % n_azi, this % n_pol)) + + ! Get delayed-nu-fission + call read_dataset(temp_1d, xsdata_grp, "delayed-nu-fission") + temp_3d = reshape(temp_1d, (/energy_groups, this % n_azi, & + this % n_pol/)) + + do ipol = 1, this % n_pol + do iazi = 1, this % n_azi + do gin = 1, energy_groups + do dg = 1, delayed_groups + + ! Set delayed-nu-fission using delayed neutron fraction + xs % delayed_nu_fission(dg, gin, iazi, ipol) = & + temp_beta(dg, gin, iazi, ipol) * & + temp_3d(gin, iazi, ipol) + end do + end do end do end do - end do - deallocate(temp_4d) + + ! Deallocate temporary delayed-nu-fission arrays + deallocate(temp_1d) + deallocate(temp_3d) + + ! If delayed-nu-fission is a (delayed_group, energy_group, + ! n_azi, n_pol) matrix, set delayed-nu-fission separately for + ! each delayed group. + else if (ndims == 4) then + + ! Get the shape of delayed-nu-fission + call get_shape(xsdata, dims) + + ! Issue error if 1st dimension not correct + if (dims(1) /= delayed_groups) then + call fatal_error("The delayed-nu-fission matrix was input & + &with a 1st dimension not equal to the number of & + &delayed groups.") + end if + + ! Issue error if 2nd dimension not correct + if (dims(2) /= energy_groups) then + call fatal_error("The delayed-nu-fission matrix was input & + &with a 2nd dimension not equal to the number of & + &energy groups.") + end if + + ! Issue warning if delayed_groups == energy_groups + if (delayed_groups == energy_groups) then + call warning("delayed-nu-fission was input as a dimension & + &4 matrix with the same number of delayed groups and & + &groups. It is important to know that OpenMC assumes & + &the dimensions in the matrix are (delayed_groups, & + &energy_groups, n_azi, n_pol). Currently, & + &delayed-nu-fission cannot be set as a group by group & + &matrix.") + end if + + ! Get delayed-nu-fission + allocate(temp_1d(delayed_groups * energy_groups * this % n_azi & + * this % n_pol)) + call read_dataset(temp_1d, xsdata_grp, "delayed-nu-fission") + xs % delayed_nu_fission = reshape(temp_1d, (/delayed_groups, & + energy_groups, this % n_azi, this % n_pol /)) + + ! Deallocate temporary array for delayed-nu-fission matrix + deallocate(temp_1d) + + else + call fatal_error("delayed-nu-fission must be provided as a & + &1D or 2D array") + end if end if - ! If we have a need* for the fission and kappa-fission x/s, get them - ! (*Need is defined as will be using it to tally) - if (get_fiss) then - if (object_exists(xsdata_grp, "fission")) then - allocate(temp_arr(groups * this % n_azi * this % n_pol)) - call read_dataset(temp_arr, xsdata_grp, "fission") - allocate(xs % fission(groups, this % n_azi, this % n_pol)) - xs % fission = reshape(temp_arr, (/groups, this % n_azi, & - this % n_pol/)) - deallocate(temp_arr) - else - call fatal_error("Fission data missing, required due to fission& - & tallies in tallies.xml file!") - end if + ! Deallocate temporary beta array + deallocate(temp_beta) + + ! chi-prompt, chi-delayed, prompt-nu-fission, and delayed-nu-fission + ! have been set; Now we will check for the rest of the XS that are + ! unique to fissionable isotopes + + ! Set fission xs + if (object_exists(xsdata_grp, "fission")) then + + ! Allocate temporary array for fission + allocate(temp_1d(energy_groups * this % n_azi * this % n_pol)) + + ! Get fission array + call read_dataset(temp_1d, xsdata_grp, "fission") + xs % fission(:, :, :) = reshape(temp_1d, (/energy_groups, & + this % n_azi, this % n_pol/)) + + ! Deallocate temporary array for fission + deallocate(temp_1d) end if - if (get_kfiss) then - if (object_exists(xsdata_grp, "kappa-fission")) then - allocate(temp_arr(groups * this % n_azi * this % n_pol)) - call read_dataset(temp_arr, xsdata_grp, "kappa-fission") - allocate(xs % k_fission(groups, this % n_azi, this % n_pol)) - xs % k_fission = reshape(temp_arr, (/groups, this % n_azi, & - this % n_pol/)) - deallocate(temp_arr) - else - call fatal_error("kappa_fission data missing, required due to & - &kappa-fission tallies in tallies.xml file!") - end if + + ! Set kappa-fission xs + if (object_exists(xsdata_grp, "kappa-fission")) then + + ! Allocate temporary array for kappa-fission + allocate(temp_1d(energy_groups * this % n_azi * this % n_pol)) + + ! Get kappa-fission array + call read_dataset(temp_1d, xsdata_grp, "kappa-fission") + xs % kappa_fission(:, :, :) = reshape(temp_1d, (/energy_groups, & + this % n_azi, this % n_pol/)) + + ! Deallocate temporary array for kappa-fission + deallocate(temp_1d) + end if + + ! Set decay rate + if (object_exists(xsdata_grp, "decay rate")) then + + ! Allocate temporary array for decay rate + allocate(temp_1d(this % n_azi * this % n_pol * delayed_groups)) + + ! Get decay rate array + call read_dataset(temp_1d, xsdata_grp, "decay rate") + xs % decay_rate(:, :, :) = reshape(temp_1d, (/delayed_groups, & + this % n_azi, this % n_pol/)) + + ! Deallocate temporary array for decay rate + deallocate(temp_1d) end if - else - xs % nu_fission = ZERO - xs % chi = ZERO end if + ! All the XS unique to fissionable isotopes have been set; Now set all + ! the generation XS + if (object_exists(xsdata_grp, "absorption")) then - allocate(temp_arr(groups * this % n_azi * this % n_pol)) - call read_dataset(temp_arr, xsdata_grp, "absorption") - allocate(xs % absorption(groups, this % n_azi, this % n_pol)) - xs % absorption = reshape(temp_arr, (/groups, this % n_azi, & - this % n_pol/)) - deallocate(temp_arr) + + ! Allocate temporary array for absorption xs + allocate(temp_1d(energy_groups * this % n_azi * this % n_pol)) + + ! Read in absorption xs + call read_dataset(temp_1d, xsdata_grp, "absorption") + + xs % absorption = reshape(temp_1d, (/energy_groups, this % n_azi, & + this % n_pol/)) + + ! Deallocate temporary array for absorption xs + deallocate(temp_1d) else call fatal_error("Must provide absorption!") end if + if (object_exists(xsdata_grp, "inverse-velocity")) then + + ! Allocate temporary array for inverse velocity + allocate(temp_1d(energy_groups * this % n_azi * this % n_pol)) + + ! Read in inverse velocity + call read_dataset(temp_1d, xsdata_grp, "inverse-velocity") + + xs % inverse_velocity = reshape(temp_1d, (/energy_groups, & + this % n_azi, this % n_pol/)) + + ! Deallocate temporary array for inverse velocity + deallocate(temp_1d) + end if + ! Get scattering data if (.not. object_exists(xsdata_grp, "scatter_data")) & call fatal_error("Must provide 'scatter_data'") + scatt_grp = open_group(xsdata_grp, 'scatter_data') + ! First get the outgoing group boundary indices if (object_exists(scatt_grp, "g_min")) then - allocate(int_arr(groups * this % n_azi * this % n_pol)) + + allocate(int_arr(energy_groups * this % n_azi * this % n_pol)) + call read_dataset(int_arr, scatt_grp, "g_min") - allocate(gmin(groups, this % n_azi, this % n_pol)) - gmin = reshape(int_arr, (/groups, this % n_azi, this % n_pol/)) + allocate(gmin(energy_groups, this % n_azi, this % n_pol)) + gmin = reshape(int_arr, (/energy_groups, this % n_azi, & + this % n_pol/)) + deallocate(int_arr) else call fatal_error("'g_min' for the scatter_data must be provided") end if + if (object_exists(scatt_grp, "g_max")) then - allocate(int_arr(groups * this % n_azi * this % n_pol)) + + allocate(int_arr(energy_groups * this % n_azi * this % n_pol)) + call read_dataset(int_arr, scatt_grp, "g_max") - allocate(gmax(groups, this % n_azi, this % n_pol)) - gmax = reshape(int_arr, (/groups, this % n_azi, this % n_pol/)) + allocate(gmax(energy_groups, this % n_azi, this % n_pol)) + gmax = reshape(int_arr, (/energy_groups, this % n_azi, & + this % n_pol/)) + deallocate(int_arr) else call fatal_error("'g_max' for the scatter_data must be provided") @@ -896,17 +1771,19 @@ module mgxs_header length = 0 do ipol = 1, this % n_pol do iazi = 1, this % n_azi - do gin = 1, groups + do gin = 1, energy_groups length = length + order_dim * (gmax(gin, iazi, ipol) - & - gmin(gin, iazi, ipol) + 1) + gmin(gin, iazi, ipol) + 1) end do end do end do + ! Allocate flattened array - allocate(temp_arr(length)) + allocate(temp_1d(length)) + if (.not. object_exists(scatt_grp, 'scatter_matrix')) & call fatal_error("'scatter_matrix' must be provided") - call read_dataset(temp_arr, scatt_grp, "scatter_matrix") + call read_dataset(temp_1d, scatt_grp, "scatter_matrix") ! Compare the number of orders given with the maximum order of the ! problem. Strip off the superfluous orders if needed. @@ -915,41 +1792,47 @@ module mgxs_header order_dim = order + 1 end if - ! Convert temp_arr to a jagged array ((gin) % data(l, gout)) for passing + ! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for passing ! to ScattData - allocate(input_scatt(groups, this % n_azi, this % n_pol)) + allocate(input_scatt(energy_groups, this % n_azi, this % n_pol)) + index = 1 do ipol = 1, this % n_pol do iazi = 1, this % n_azi - do gin = 1, groups + do gin = 1, energy_groups allocate(input_scatt(gin, iazi, ipol) % data(order_dim, & gmin(gin, iazi, ipol):gmax(gin, iazi, ipol))) do gout = gmin(gin, iazi, ipol), gmax(gin, iazi, ipol) do l = 1, order_dim input_scatt(gin, iazi, ipol) % data(l, gout) = & - temp_arr(index) + temp_1d(index) index = index + 1 end do ! gout end do ! order end do ! gin end do ! iazi end do ! ipol - deallocate(temp_arr) + + deallocate(temp_1d) ! Finally convert the legendre to tabular if needed - allocate(scatt_coeffs(groups, this % n_azi, this % n_pol)) + allocate(scatt_coeffs(energy_groups, this % n_azi, this % n_pol)) + if (this % scatter_format == ANGLE_LEGENDRE .and. & legendre_to_tabular) then + this % scatter_format = ANGLE_TABULAR order_dim = legendre_to_tabular_points order = order_dim dmu = TWO / real(order - 1, 8) + do ipol = 1, this % n_pol do iazi = 1, this % n_azi - do gin = 1, groups + do gin = 1, energy_groups allocate(scatt_coeffs(gin, iazi, ipol) % data(& order_dim, & gmin(gin, iazi, ipol):gmax(gin, iazi, ipol))) + do gout = gmin(gin, iazi, ipol), gmax(gin, iazi, ipol) norm = ZERO do imu = 1, order_dim @@ -960,12 +1843,15 @@ module mgxs_header else mu = -ONE + real(imu - 1, 8) * dmu end if + scatt_coeffs(gin, iazi, ipol) % data(imu, gout) = & evaluate_legendre(& input_scatt(gin, iazi, ipol) % data(:, gout), mu) + ! Ensure positivity of distribution if (scatt_coeffs(gin, iazi, ipol) % data(imu, gout) < ZERO) & scatt_coeffs(gin, iazi, ipol) % data(imu, gout) = ZERO + ! And accrue the integral if (imu > 1) then norm = norm + HALF * dmu * & @@ -973,6 +1859,7 @@ module mgxs_header scatt_coeffs(gin, iazi, ipol) % data(imu, gout)) end if end do ! mu + ! Now that we have the integral, lets ensure that the distribution ! is normalized such that it preserves the original scattering xs if (norm > ZERO) then @@ -989,7 +1876,7 @@ module mgxs_header ! Sticking with current representation, carry forward do ipol = 1, this % n_pol do iazi = 1, this % n_azi - do gin = 1, groups + do gin = 1, energy_groups allocate(scatt_coeffs(gin, iazi, ipol) % data(order_dim, & gmin(gin, iazi, ipol):gmax(gin, iazi, ipol))) scatt_coeffs(gin, iazi, ipol) % data(:, :) = & @@ -998,46 +1885,54 @@ module mgxs_header end do end do end if + deallocate(input_scatt) ! Now get the multiplication matrix - if (object_exists(scatt_grp, 'multiplicity_matrix')) then + if (object_exists(scatt_grp, 'multiplicity matrix')) then + ! Now use this information to find the length of a container array ! to hold the flattened data length = 0 + do ipol = 1, this % n_pol do iazi = 1, this % n_azi - do gin = 1, groups + do gin = 1, energy_groups length = length + (gmax(gin, iazi, ipol) - gmin(gin, iazi, ipol) + 1) end do end do end do + ! Allocate flattened array - allocate(temp_arr(length)) - call read_dataset(temp_arr, scatt_grp, "multiplicity_matrix") - ! Convert temp_arr to a jagged array ((gin) % data(gout)) for passing + allocate(temp_1d(length)) + call read_dataset(temp_1d, scatt_grp, "multiplicity matrix") + + ! Convert temp_1d to a jagged array ((gin) % data(gout)) for passing ! to ScattData - allocate(temp_mult(groups, this % n_azi, this % n_pol)) + allocate(temp_mult(energy_groups, this % n_azi, this % n_pol)) + index = 1 do ipol = 1, this % n_pol do iazi = 1, this % n_azi - do gin = 1, groups + do gin = 1, energy_groups allocate(temp_mult(gin, iazi, ipol) % data( & gmin(gin, iazi, ipol):gmax(gin, iazi, ipol))) do gout = gmin(gin, iazi, ipol), gmax(gin, iazi, ipol) - temp_mult(gin, iazi, ipol) % data(gout) = temp_arr(index) + temp_mult(gin, iazi, ipol) % data(gout) = temp_1d(index) index = index + 1 end do end do end do end do - deallocate(temp_arr) + deallocate(temp_1d) else - allocate(temp_mult(groups, this % n_azi, this % n_pol)) + + allocate(temp_mult(energy_groups, this % n_azi, this % n_pol)) + ! Default to multiplicities of 1.0 do ipol = 1, this % n_pol do iazi = 1, this % n_azi - do gin = 1, groups + do gin = 1, energy_groups allocate(temp_mult(gin, iazi, ipol) % data( & gmin(gin, iazi, ipol):gmax(gin, iazi, ipol))) temp_mult(gin, iazi, ipol) % data = ONE @@ -1048,8 +1943,10 @@ module mgxs_header ! Allocate and initialize our ScattData Object. allocate(xs % scatter(this % n_azi, this % n_pol)) + do ipol = 1, this % n_pol do iazi = 1, this % n_azi + ! Allocate and initialize our ScattData Object. if (this % scatter_format == ANGLE_HISTOGRAM) then allocate(ScattDataHistogram :: xs % scatter(iazi, ipol) % obj) @@ -1071,7 +1968,7 @@ module mgxs_header ! (This may happen with Helium data) do ipol = 1, this % n_pol do iazi = 1, this % n_azi - do gin = 1, groups + do gin = 1, energy_groups if (xs % absorption(gin, iazi, ipol) == ZERO) then xs % absorption(gin, iazi, ipol) = 1E-10_8 end if @@ -1079,13 +1976,14 @@ module mgxs_header end do end do - allocate(xs % total(groups, this % n_azi, this % n_pol)) if (object_exists(xsdata_grp, "total")) then - allocate(temp_arr(groups * this % n_azi * this % n_pol)) - call read_dataset(temp_arr, xsdata_grp, "total") - xs % total = reshape(temp_arr, (/groups, this % n_azi, & - this % n_pol/)) - deallocate(temp_arr) + + allocate(temp_1d(energy_groups * this % n_azi * this % n_pol)) + call read_dataset(temp_1d, xsdata_grp, "total") + xs % total = reshape(temp_1d, (/energy_groups, this % n_azi, & + this % n_pol/)) + + deallocate(temp_1d) else do ipol = 1, this % n_pol do iazi = 1, this % n_azi @@ -1099,7 +1997,7 @@ module mgxs_header ! the tally routines do ipol = 1, this % n_pol do iazi = 1, this % n_azi - do gin = 1, groups + do gin = 1, energy_groups if (xs % total(gin, iazi, ipol) == ZERO) then xs % total(gin, iazi, ipol) = 1E-10_8 end if @@ -1107,20 +2005,11 @@ module mgxs_header end do end do - ! Get kinetics data - if (object_exists(xsdata_grp, "inverse_velocities")) then - allocate(xs % inv_vel(groups, this % n_azi, this % n_pol)) - allocate(temp_arr(groups * this % n_azi * this % n_pol)) - call read_dataset(temp_arr, xsdata_grp, "inverse_velocities") - xs % inv_vel = reshape(temp_arr, (/groups, this % n_azi, & - this % n_pol/)) - deallocate(temp_arr) - end if - ! Close the groups we have opened and deallocate call close_group(xsdata_grp) call close_group(scatt_grp) deallocate(scatt_coeffs, temp_mult) + end associate ! xs end do ! Temperatures end subroutine mgxsang_from_hdf5 @@ -1147,6 +2036,8 @@ module mgxs_header else this % name = trim(mat % name) end if + + ! Set whether this material is fissionable this % fissionable = mat % fissionable ! The following info we should initialize, but we dont need it nor @@ -1154,6 +2045,7 @@ module mgxs_header this % awr = -ONE allocate(this % kTs(temps % size())) + do t = 1, temps % size() this % kTs(t) = temps % data(t) end do @@ -1169,12 +2061,14 @@ module mgxs_header ! Determine the scattering type of our data and ensure all scattering orders ! are the same. scatter_format = nuclides(mat % nuclide(1)) % obj % scatter_format + select type(nuc => nuclides(mat % nuclide(1)) % obj) type is (MgxsIso) order = size(nuc % xs(1) % scatter % dist(1) % data, dim=1) type is (MgxsAngle) order = size(nuc % xs(1) % scatter(1, 1) % obj % dist(1) % data, dim=1) end select + ! If we have tabular only data, then make sure all datasets have same size if (scatter_format == ANGLE_HISTOGRAM) then ! Check all scattering data to ensure it is the same size @@ -1190,8 +2084,10 @@ module mgxs_header & same length!") end select end do + ! Ok, got our order, store the dimensionality order_dim = order + else if (scatter_format == ANGLE_TABULAR) then ! Check all scattering data to ensure it is the same size do i = 2, mat % n_nuclides @@ -1206,11 +2102,15 @@ module mgxs_header & same length!") end select end do + ! Ok, got our order, store the dimensionality order_dim = order + else if (scatter_format == ANGLE_LEGENDRE) then + ! Need to determine the maximum scattering order of all data in this material mat_max_order = 0 + do i = 1, mat % n_nuclides select type(nuc => nuclides(mat % nuclide(i)) % obj) type is (MgxsIso) @@ -1230,19 +2130,21 @@ module mgxs_header ! Now need to compare this material maximum scattering order with ! the problem wide max scatt order and use whichever is lower order = min(mat_max_order, max_order) + ! Ok, got our order, store the dimensionality order_dim = order end if end subroutine mgxs_combine - subroutine mgxsiso_combine(this, temps, mat, nuclides, groups, max_order, & - tolerance, method) - class(MgxsIso), intent(inout) :: this ! The Mgxs to initialize - type(VectorReal), intent(in) :: temps ! Temperatures to obtain [eV] - type(Material), pointer, intent(in) :: mat ! base material - type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from - integer, intent(in) :: groups ! Number of E groups + subroutine mgxsiso_combine(this, temps, mat, nuclides, energy_groups, & + delayed_groups, max_order, tolerance, method) + class(MgxsIso), intent(inout) :: this ! The Mgxs to initialize + type(VectorReal), intent(in) :: temps ! Temperatures to obtain [MeV] + type(Material), pointer, intent(in) :: mat ! base material + type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from + integer, intent(in) :: energy_groups ! Number of energy groups + integer, intent(in) :: delayed_groups ! Number of delayed groups integer, intent(in) :: max_order ! Maximum requested order real(8), intent(in) :: tolerance ! Tolerance on method integer, intent(in) :: method ! Type of temperature access @@ -1250,6 +2152,7 @@ module mgxs_header integer :: i ! loop index over nuclides integer :: t ! Index in to temps integer :: gin, gout ! group indices + integer :: dg ! delayed group index real(8) :: atom_density ! atom density of a nuclide real(8) :: norm, nuscatt integer :: order_dim, nuc_order_dim @@ -1268,33 +2171,59 @@ module mgxs_header call mgxs_combine(this, temps, mat, nuclides, max_order, scatter_format, & order_dim) + ! Set the number of delayed groups + this % num_delayed_groups = delayed_groups + ! Create the Xs Data for each temperature TEMP_LOOP: do t = 1, temps % size() + ! Allocate and initialize the data needed for macro_xs(i_mat) object - allocate(this % xs(t) % total(groups)) + allocate(this % xs(t) % total(energy_groups)) this % xs(t) % total(:) = ZERO - allocate(this % xs(t) % absorption(groups)) + + allocate(this % xs(t) % absorption(energy_groups)) this % xs(t) % absorption(:) = ZERO - allocate(this % xs(t) % fission(groups)) + + allocate(this % xs(t) % fission(energy_groups)) this % xs(t) % fission(:) = ZERO - allocate(this % xs(t) % k_fission(groups)) - this % xs(t) % k_fission(:) = ZERO - allocate(this % xs(t) % nu_fission(groups)) - this % xs(t) % nu_fission(:) = ZERO - allocate(this % xs(t) % chi(groups,groups)) - this % xs(t) % chi(:, :) = ZERO - allocate(this % xs(t) % inv_vel(groups)) - this % xs(t) % inv_vel(:) = ZERO - allocate(temp_mult(groups,groups)) + + allocate(this % xs(t) % kappa_fission(energy_groups)) + this % xs(t) % kappa_fission(:) = ZERO + + allocate(this % xs(t) % prompt_nu_fission(energy_groups)) + this % xs(t) % prompt_nu_fission(:) = ZERO + + allocate(this % xs(t) % delayed_nu_fission(delayed_groups, & + energy_groups)) + this % xs(t) % delayed_nu_fission(:, :) = ZERO + + allocate(this % xs(t) % chi_prompt(energy_groups, energy_groups)) + this % xs(t) % chi_prompt(:, :) = ZERO + + allocate(this % xs(t) % chi_delayed(delayed_groups, energy_groups, & + energy_groups)) + this % xs(t) % chi_delayed(:, :, :) = ZERO + + allocate(this % xs(t) % inverse_velocity(energy_groups)) + this % xs(t) % inverse_velocity(:) = ZERO + + allocate(this % xs(t) % decay_rate(delayed_groups)) + this % xs(t) % decay_rate(:) = ZERO + + allocate(temp_mult(energy_groups, energy_groups)) temp_mult(:, :) = ZERO - allocate(mult_num(groups,groups)) + + allocate(mult_num(energy_groups, energy_groups)) mult_num(:, :) = ZERO - allocate(mult_denom(groups,groups)) + + allocate(mult_denom(energy_groups, energy_groups)) mult_denom(:, :) = ZERO - allocate(scatt_coeffs(order_dim,groups,groups)) + + allocate(scatt_coeffs(order_dim, energy_groups, energy_groups)) scatt_coeffs(:, :, :) = ZERO this % scatter_format = scatter_format + if (scatter_format == ANGLE_LEGENDRE) then allocate(ScattDataLegendre :: this % xs(t) % scatter) else if (scatter_format == ANGLE_TABULAR) then @@ -1306,16 +2235,20 @@ module mgxs_header ! Add contribution from each nuclide in material NUC_LOOP: do i = 1, mat % n_nuclides associate(nuc => nuclides(mat % nuclide(i)) % obj) - ! Copy atom density of nuclide in material + + ! Copy atom density of nuclide in material atom_density = mat % atom_density(i) select case (method) case (TEMPERATURE_NEAREST) + ! Determine actual temperatures to read temp_desired = temps % data(i) allocate(nuc_ts(1)) + nuc_ts(1) = minloc(abs(nuc % kTs - temp_desired), dim=1) temp_actual = nuc % kTs(nuc_ts(1)) + if (abs(temp_actual - temp_desired) >= K_BOLTZMANN * tolerance) then call fatal_error("MGXS library does not contain cross sections & &for " // trim(this % name) // " at or near " // & @@ -1323,11 +2256,13 @@ module mgxs_header end if case (TEMPERATURE_INTERPOLATION) + ! If temperature interpolation or multipole is selected, get a ! list of bounding temperatures for each actual temperature ! present in the model temp_desired = temps % data(i) allocate(nuc_ts(2)) + do j = 1, size(nuc % kTs) - 1 if (nuc % kTs(j) <= temp_desired .and. & temp_desired < nuc % kTs(j + 1)) then @@ -1344,7 +2279,9 @@ module mgxs_header select type(nuc) type is (MgxsIso) do j = 1, size(nuc_ts) + nuc_t = nuc_ts(j) + if (size(nuc_ts) == 1) then interp = ONE else if (j == 1) then @@ -1358,25 +2295,38 @@ module mgxs_header ! Add contributions to total, absorption, and fission data (if necessary) this % xs(t) % total = this % xs(t) % total + & atom_density * nuc % xs(nuc_t) % total * interp + this % xs(t) % absorption = this % xs(t) % absorption + & atom_density * nuc % xs(nuc_t) % absorption * interp + + this % xs(t) % decay_rate = this % xs(t) % decay_rate + & + atom_density * nuc % xs(nuc_t) % decay_rate * interp + + this % xs(t) % inverse_velocity = & + this % xs(t) % inverse_velocity + & + atom_density * nuc % xs(nuc_t) % inverse_velocity * interp + if (nuc % fissionable) then - this % xs(t) % chi = this % xs(t) % chi + & - atom_density * nuc % xs(nuc_t) % chi * interp - this % xs(t) % nu_fission = this % xs(t) % nu_fission + & - atom_density * nuc % xs(nuc_t) % nu_fission * interp - if (allocated(nuc % xs(nuc_t) % fission)) then - this % xs(t) % fission = this % xs(t) % fission + & - atom_density * nuc % xs(nuc_t) % fission * interp - end if - if (allocated(nuc % xs(nuc_t) % k_fission)) then - this % xs(t) % k_fission = this % xs(t) % k_fission + & - atom_density * nuc % xs(nuc_t) % k_fission * interp - end if - end if - if (allocated(nuc % xs(nuc_t) % inv_vel)) then - this % xs(t) % inv_vel = this % xs(t) % inv_vel + & - atom_density * nuc % xs(nuc_t) % inv_vel * interp + + this % xs(t) % chi_prompt = this % xs(t) % chi_prompt + & + atom_density * nuc % xs(nuc_t) % chi_prompt * interp + + this % xs(t) % chi_delayed = this % xs(t) % chi_delayed + & + atom_density * nuc % xs(nuc_t) % chi_delayed * interp + + this % xs(t) % prompt_nu_fission = this % xs(t) % & + prompt_nu_fission + atom_density * nuc % xs(nuc_t) % & + prompt_nu_fission * interp + + this % xs(t) % delayed_nu_fission = this % xs(t) % & + delayed_nu_fission + atom_density * nuc % xs(nuc_t) % & + delayed_nu_fission * interp + + this % xs(t) % fission = this % xs(t) % fission + & + atom_density * nuc % xs(nuc_t) % fission * interp + + this % xs(t) % kappa_fission = this % xs(t) % kappa_fission +& + atom_density * nuc % xs(nuc_t) % kappa_fission * interp end if ! We will next gather the multiplicity and scattering matrices. @@ -1396,12 +2346,15 @@ module mgxs_header ! sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'})) ! nuscatt_{i,g,g'} can be reconstructed from scatter % energy and ! scatter % scattxs - do gin = 1, groups + do gin = 1, energy_groups do gout = nuc % xs(nuc_t) % scatter % gmin(gin), nuc % xs(nuc_t) % scatter % gmax(gin) + nuscatt = nuc % xs(nuc_t) % scatter % scattxs(gin) * & nuc % xs(nuc_t) % scatter % energy(gin) % data(gout) + mult_num(gout, gin) = mult_num(gout, gin) + atom_density * & nuscatt * interp + if (nuc % xs(nuc_t) % scatter % mult(gin) % data(gout) > ZERO) then mult_denom(gout, gin) = mult_denom(gout,gin) + atom_density * & nuscatt / nuc % xs(nuc_t) % scatter % mult(gin) % data(gout) * & @@ -1417,9 +2370,11 @@ module mgxs_header ! Get the complete scattering matrix nuc_order_dim = size(nuc % xs(nuc_t) % scatter % dist(1) % data, dim=1) nuc_order_dim = min(nuc_order_dim, order_dim) + call nuc % xs(nuc_t) % scatter % get_matrix(nuc_order_dim, & nuc_matrix) - do gin = 1, groups + + do gin = 1, energy_groups do gout = nuc % xs(nuc_t) % scatter % gmin(gin), & nuc % xs(nuc_t) % scatter % gmax(gin) scatt_coeffs(1:nuc_order_dim, gout, gin) = & @@ -1435,8 +2390,8 @@ module mgxs_header end select ! Obtain temp_mult - do gin = 1, groups - do gout = 1, groups + do gin = 1, energy_groups + do gout = 1, energy_groups if (mult_denom(gout, gin) > ZERO) then temp_mult(gout, gin) = mult_num(gout, gin) / mult_denom(gout, gin) else @@ -1455,12 +2410,23 @@ module mgxs_header ! Now normalize chi if (mat % fissionable) then - do gin = 1, groups - norm = sum(this % xs(t) % chi(:, gin)) + do gin = 1, energy_groups + norm = sum(this % xs(t) % chi_prompt(:, gin)) if (norm > ZERO) then - this % xs(t) % chi(:, gin) = this % xs(t) % chi(:, gin) / norm + this % xs(t) % chi_prompt(:, gin) = & + this % xs(t) % chi_prompt(:, gin) / norm end if end do + + do dg = 1, delayed_groups + do gin = 1, energy_groups + norm = sum(this % xs(t) % chi_delayed(dg, :, gin)) + if (norm > ZERO) then + this % xs(t) % chi_delayed(dg, :, gin) = & + this % xs(t) % chi_delayed(dg, :, gin) / norm + end if + end do + end do end if ! Deallocate temporaries @@ -1472,13 +2438,14 @@ module mgxs_header end subroutine mgxsiso_combine - subroutine mgxsang_combine(this, temps, mat, nuclides, groups, max_order, & - tolerance, method) + subroutine mgxsang_combine(this, temps, mat, nuclides, energy_groups, & + delayed_groups, max_order, tolerance, method) class(MgxsAngle), intent(inout) :: this ! The Mgxs to initialize type(VectorReal), intent(in) :: temps ! Temperatures to obtain type(Material), pointer, intent(in) :: mat ! base material type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from - integer, intent(in) :: groups ! Number of E groups + integer, intent(in) :: energy_groups ! Number of energy groups + integer, intent(in) :: delayed_groups ! Number of delayed groups integer, intent(in) :: max_order ! Maximum requested order real(8), intent(in) :: tolerance ! Tolerance on method integer, intent(in) :: method ! Type of temperature access @@ -1486,6 +2453,7 @@ module mgxs_header integer :: i ! loop index over nuclides integer :: t ! temperature loop index integer :: gin, gout ! group indices + integer :: dg ! delayed group index real(8) :: atom_density ! atom density of a nuclide integer :: ipol, iazi, n_pol, n_azi real(8) :: norm, nuscatt @@ -1503,20 +2471,27 @@ module mgxs_header ! Set the meta-data call mgxs_combine(this, temps, mat, nuclides, max_order, scatter_format, & - order_dim) + order_dim) + + ! Set the number of delayed groups + this % num_delayed_groups = delayed_groups ! Get the number of each polar and azi angles and make sure all the ! NuclideAngle types have the same number of these angles n_pol = -1 n_azi = -1 + do i = 1, mat % n_nuclides select type(nuc => nuclides(mat % nuclide(i)) % obj) type is (MgxsAngle) + if (n_pol == -1) then n_pol = nuc % n_pol n_azi = nuc % n_azi + allocate(this % polar(n_pol)) this % polar(:) = nuc % polar(:) + allocate(this % azimuthal(n_azi)) this % azimuthal(:) = nuc % azimuthal(:) else @@ -1529,31 +2504,55 @@ module mgxs_header ! Create the Xs Data for each temperature TEMP_LOOP: do t = 1, temps % size() + ! Allocate and initialize the data needed for macro_xs(i_mat) object - allocate(this % xs(t) % total(groups, n_azi, n_pol)) + allocate(this % xs(t) % total(energy_groups, n_azi, n_pol)) this % xs(t) % total = ZERO - allocate(this % xs(t) % absorption(groups, n_azi, n_pol)) + + allocate(this % xs(t) % absorption(energy_groups, n_azi, n_pol)) this % xs(t) % absorption = ZERO - allocate(this % xs(t) % fission(groups, n_azi, n_pol)) + + allocate(this % xs(t) % fission(energy_groups, n_azi, n_pol)) this % xs(t) % fission = ZERO - allocate(this % xs(t) % k_fission(groups, n_azi, n_pol)) - this % xs(t) % k_fission = ZERO - allocate(this % xs(t) % nu_fission(groups, n_azi, n_pol)) - this % xs(t) % nu_fission = ZERO - allocate(this % xs(t) % chi(groups, groups, n_azi, n_pol)) - this % xs(t) % chi = ZERO - allocate(this % xs(t) % inv_vel(groups, n_azi, n_pol)) - this % xs(t) % inv_vel = ZERO - allocate(temp_mult(groups, groups, n_azi, n_pol)) + + allocate(this % xs(t) % decay_rate(delayed_groups, n_azi, n_pol)) + this % xs(t) % decay_rate = ZERO + + allocate(this % xs(t) % inverse_velocity(energy_groups, n_azi, n_pol)) + this % xs(t) % inverse_velocity = ZERO + + allocate(this % xs(t) % kappa_fission(energy_groups, n_azi, n_pol)) + this % xs(t) % kappa_fission = ZERO + + allocate(this % xs(t) % prompt_nu_fission(energy_groups, n_azi, n_pol)) + this % xs(t) % prompt_nu_fission = ZERO + + allocate(this % xs(t) % delayed_nu_fission(delayed_groups, & + energy_groups, n_azi, n_pol)) + this % xs(t) % delayed_nu_fission = ZERO + + allocate(this % xs(t) % chi_prompt(energy_groups, energy_groups, & + n_azi, n_pol)) + this % xs(t) % chi_prompt = ZERO + + allocate(this % xs(t) % chi_delayed(delayed_groups, energy_groups, & + energy_groups, n_azi, n_pol)) + this % xs(t) % chi_delayed = ZERO + + allocate(temp_mult(energy_groups, energy_groups, n_azi, n_pol)) temp_mult = ZERO - allocate(mult_num(groups, groups, n_azi, n_pol)) + + allocate(mult_num(energy_groups, energy_groups, n_azi, n_pol)) mult_num = ZERO - allocate(mult_denom(groups, groups, n_azi, n_pol)) + + allocate(mult_denom(energy_groups, energy_groups, n_azi, n_pol)) mult_denom = ZERO - allocate(scatt_coeffs(order_dim, groups, groups, n_azi, n_pol)) + + allocate(scatt_coeffs(order_dim, energy_groups, energy_groups, n_azi, n_pol)) scatt_coeffs = ZERO allocate(this % xs(t) % scatter(n_azi, n_pol)) + do ipol = 1, n_pol do iazi = 1, n_azi if (scatter_format == ANGLE_LEGENDRE) then @@ -1575,11 +2574,14 @@ module mgxs_header select case (method) case (TEMPERATURE_NEAREST) + ! Determine actual temperatures to read temp_desired = temps % data(i) allocate(nuc_ts(1)) + nuc_ts(1) = minloc(abs(nuc % kTs - temp_desired), dim=1) temp_actual = nuc % kTs(nuc_ts(1)) + if (abs(temp_actual - temp_desired) >= K_BOLTZMANN * tolerance) then call fatal_error("MGXS library does not contain cross sections & &for " // trim(this % name) // " at or near " // & @@ -1587,11 +2589,13 @@ module mgxs_header end if case (TEMPERATURE_INTERPOLATION) + ! If temperature interpolation or multipole is selected, get a ! list of bounding temperatures for each actual temperature ! present in the model temp_desired = temps % data(i) allocate(nuc_ts(2)) + do j = 1, size(nuc % kTs) - 1 if (nuc % kTs(j) <= temp_desired .and. & temp_desired < nuc % kTs(j + 1)) then @@ -1611,7 +2615,9 @@ module mgxs_header select type(nuc) type is (MgxsAngle) do j = 1, size(nuc_ts) + nuc_t = nuc_ts(j) + if (size(nuc_ts) == 1) then interp = ONE else if (j == 1) then @@ -1620,30 +2626,45 @@ module mgxs_header else interp = ONE - interp end if + ! Perform our operations which depend upon the type ! Add contributions to total, absorption, and fission data ! (if necessary) this % xs(t) % total = this % xs(t) % total + & atom_density * nuc % xs(nuc_t) % total * interp + this % xs(t) % absorption = this % xs(t) % absorption + & atom_density * nuc % xs(nuc_t) % absorption * interp - if (allocated(nuc % xs(nuc_t) % inv_vel)) then - this % xs(t) % inv_vel = this % xs(t) % inv_vel + & - atom_density * nuc % xs(nuc_t) % inv_vel * interp - end if + + this % xs(t) % decay_rate = this % xs(t) % decay_rate + & + atom_density * nuc % xs(nuc_t) % decay_rate * interp + + this % xs(t) % inverse_velocity = & + this % xs(t) % inverse_velocity + & + atom_density * nuc % xs(nuc_t) % inverse_velocity * interp + if (nuc % fissionable) then - this % xs(t) % chi = this % xs(t) % chi + & - atom_density * nuc % xs(nuc_t) % chi * interp - this % xs(t) % nu_fission = this % xs(t) % nu_fission + & - atom_density * nuc % xs(nuc_t) % nu_fission * interp - if (allocated(nuc % xs(nuc_t) % fission)) then - this % xs(t) % fission = this % xs(t) % fission + & - atom_density * nuc % xs(nuc_t) % fission * interp - end if - if (allocated(nuc % xs(nuc_t) % k_fission)) then - this % xs(t) % k_fission = this % xs(t) % k_fission + & - atom_density * nuc % xs(nuc_t) % k_fission * interp - end if + + this % xs(t) % chi_prompt = this % xs(t) % chi_prompt + & + atom_density * nuc % xs(nuc_t) % chi_prompt * interp + + this % xs(t) % chi_delayed = this % xs(t) % chi_delayed + & + atom_density * nuc % xs(nuc_t) % chi_delayed * interp + + this % xs(t) % prompt_nu_fission = & + this % xs(t) % prompt_nu_fission + atom_density * & + nuc % xs(nuc_t) % prompt_nu_fission * interp + + this % xs(t) % delayed_nu_fission = & + this % xs(t) % delayed_nu_fission + atom_density * & + nuc % xs(nuc_t) % delayed_nu_fission * interp + + this % xs(t) % fission = this % xs(t) % fission + & + atom_density * nuc % xs(nuc_t) % fission * interp + + this % xs(t) % kappa_fission = this % xs(t) % kappa_fission & + + atom_density * nuc % xs(nuc_t) % kappa_fission * interp + end if ! We will next gather the multiplicity and scattering matrices. @@ -1665,14 +2686,17 @@ module mgxs_header ! scatter % scattxs do ipol = 1, n_pol do iazi = 1, n_azi - do gin = 1, groups + do gin = 1, energy_groups do gout = nuc % xs(nuc_t) % scatter(iazi, ipol) %obj % gmin(gin), & nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % gmax(gin) + nuscatt = nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % scattxs(gin) * & nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % energy(gin) % data(gout) + mult_num(gout, gin, iazi, ipol) = & mult_num(gout, gin, iazi, ipol) + atom_density * & nuscatt * interp + if (nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % mult(gin) % data(gout) > ZERO) then mult_denom(gout, gin, iazi, ipol) = & mult_denom(gout, gin, iazi, ipol) + atom_density * & @@ -1692,9 +2716,11 @@ module mgxs_header size(nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % & dist(1) % data, dim=1) nuc_order_dim = min(nuc_order_dim, order_dim) + call nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % & get_matrix(nuc_order_dim, nuc_matrix) - do gin = 1, groups + + do gin = 1, energy_groups do gout = & nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % gmin(gin), & nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % gmax(gin) @@ -1716,8 +2742,8 @@ module mgxs_header do ipol = 1, n_pol do iazi = 1, n_azi ! Obtain temp_mult - do gin = 1, groups - do gout = 1, groups + do gin = 1, energy_groups + do gout = 1, energy_groups if (mult_denom(gout, gin, iazi, ipol) > ZERO) then temp_mult(gout, gin, iazi, ipol) = & mult_num(gout, gin, iazi, ipol) / & @@ -1745,21 +2771,38 @@ module mgxs_header if (mat % fissionable) then do ipol = 1, n_pol do iazi = 1, n_azi - do gin = 1, groups - norm = sum(this % xs(t) % chi(:, gin, iazi, ipol)) + do gin = 1, energy_groups + norm = sum(this % xs(t) % chi_prompt(:, gin, iazi, ipol)) if (norm > ZERO) then - this % xs(t) % chi(:, gin, iazi, ipol) = & - this % xs(t) % chi(:, gin, iazi, ipol) / norm + this % xs(t) % chi_prompt(:, gin, iazi, ipol) = & + this % xs(t) % chi_prompt(:, gin, iazi, ipol) / norm end if end do end do end do + + do dg = 1, delayed_groups + do ipol = 1, n_pol + do iazi = 1, n_azi + do gin = 1, energy_groups + norm = sum(this % xs(t) % chi_delayed(dg, :, gin, iazi, ipol)) + if (norm > ZERO) then + this % xs(t) % chi_delayed(dg, :, gin, iazi, ipol) = & + this % xs(t) % chi_delayed(dg, :, gin, iazi, ipol)& + / norm + end if + end do + end do + end do + end do + end if ! Deallocate temporaries deallocate(scatt_coeffs, temp_mult, mult_num, mult_denom) end associate ! nuc end do NUC_LOOP + end do TEMP_LOOP end subroutine mgxsang_combine @@ -1768,44 +2811,80 @@ module mgxs_header ! MGXS*_GET_XS returns the requested data cross section data !=============================================================================== - pure function mgxsiso_get_xs(this, xstype, gin, gout, uvw, mu) result(xs) + pure function mgxsiso_get_xs(this, xstype, gin, gout, uvw, mu, dg) result(xs) class(MgxsIso), intent(in) :: this ! The Xs to get data from character(*) , intent(in) :: xstype ! Type of xs requested integer, intent(in) :: gin ! Incoming Energy group integer, optional, intent(in) :: gout ! Outgoing Energy group real(8), optional, intent(in) :: uvw(3) ! Requested Angle real(8), optional, intent(in) :: mu ! Change in angle + integer, optional, intent(in) :: dg ! Delayed group real(8) :: xs ! Requested x/s integer :: t ! temperature index t = this % index_temp select case(xstype) + case('total') xs = this % xs(t) % total(gin) + case('absorption') xs = this % xs(t) % absorption(gin) + case('fission') - if (allocated(this % xs(t) % fission)) then - xs = this % xs(t) % fission(gin) + xs = this % xs(t) % fission(gin) + + case('kappa-fission') + xs = this % xs(t) % kappa_fission(gin) + + case('inverse-velocity') + xs = this % xs(t) % inverse_velocity(gin) + + case('decay rate') + if (present(dg)) then + xs = this % xs(t) % decay_rate(dg) else - xs = ZERO + xs = this % xs(t) % decay_rate(1) end if - case('kappa_fission') - if (allocated(this % xs(t) % k_fission)) then - xs = this % xs(t) % k_fission(gin) + + case('prompt-nu-fission') + xs = this % xs(t) % prompt_nu_fission(gin) + + case('delayed-nu-fission') + if (present(dg)) then + xs = this % xs(t) % delayed_nu_fission(dg, gin) else - xs = ZERO + xs = sum(this % xs(t) % delayed_nu_fission(:, gin)) end if - case('nu_fission') - xs = this % xs(t) % nu_fission(gin) - case('chi') + + case('nu-fission') + xs = this % xs(t) % prompt_nu_fission(gin) + & + sum(this % xs(t) % delayed_nu_fission(:, gin)) + + case('chi-prompt') if (present(gout)) then - xs = this % xs(t) % chi(gout,gin) + xs = this % xs(t) % chi_prompt(gout,gin) else ! Not sure youd want a 1 or a 0, but here you go! - xs = sum(this % xs(t) % chi(:, gin)) + xs = sum(this % xs(t) % chi_prompt(:, gin)) end if + + case('chi-delayed') + if (present(gout)) then + if (present(dg)) then + xs = this % xs(t) % chi_delayed(dg, gout, gin) + else + xs = this % xs(t) % chi_delayed(1, gout, gin) + end if + else + if (present(dg)) then + xs = sum(this % xs(t) % chi_delayed(dg, :, gin)) + else + xs = sum(this % xs(t) % chi_delayed(dg, :, gin)) + end if + end if + case('scatter') if (present(gout)) then if (gout < this % xs(t) % scatter % gmin(gin) .or. & @@ -1818,6 +2897,7 @@ module mgxs_header else xs = this % xs(t) % scatter % scattxs(gin) end if + case('scatter/mult') if (present(gout)) then if (gout < this % xs(t) % scatter % gmin(gin) .or. & @@ -1833,6 +2913,7 @@ module mgxs_header (dot_product(this % xs(t) % scatter % mult(gin) % data, & this % xs(t) % scatter % energy(gin) % data)) end if + case('scatter*f_mu/mult','scatter*f_mu') if (present(gout)) then if (gout < this % xs(t) % scatter % gmin(gin) .or. & @@ -1853,20 +2934,21 @@ module mgxs_header ! user of this code wants the complete 1-outgoing group distribution ! which Im not sure what they would do with that. end if - case('inv_vel') - xs = this % xs(t) % inv_vel(gin) + case default xs = ZERO end select + end function mgxsiso_get_xs - pure function mgxsang_get_xs(this, xstype, gin, gout, uvw, mu) result(xs) + pure function mgxsang_get_xs(this, xstype, gin, gout, uvw, mu, dg) result(xs) class(MgxsAngle), intent(in) :: this ! The Mgxs to initialize character(*) , intent(in) :: xstype ! Type of xs requested integer, intent(in) :: gin ! Incoming Energy group integer, optional, intent(in) :: gout ! Outgoing Energy group real(8), optional, intent(in) :: uvw(3) ! Requested Angle real(8), optional, intent(in) :: mu ! Change in angle + integer, optional, intent(in) :: dg ! Delayed group real(8) :: xs ! Requested x/s integer :: iazi, ipol, t @@ -1874,33 +2956,70 @@ module mgxs_header t = this % index_temp if (present(uvw)) then + call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol) + select case(xstype) + case('total') xs = this % xs(t) % total(gin, iazi, ipol) + case('absorption') xs = this % xs(t) % absorption(gin, iazi, ipol) + case('fission') - if (allocated(this % xs(t) % fission)) then - xs = this % xs(t) % fission(gin, iazi, ipol) + xs = this % xs(t) % fission(gin, iazi, ipol) + + case('kappa-fission') + xs = this % xs(t) % kappa_fission(gin, iazi, ipol) + + case('prompt-nu-fission') + xs = this % xs(t) % prompt_nu_fission(gin, iazi, ipol) + + case('delayed-nu-fission') + if (present(dg)) then + xs = this % xs(t) % delayed_nu_fission(dg, gin, iazi, ipol) else - xs = ZERO + xs = sum(this % xs(t) % delayed_nu_fission(:, gin, iazi, ipol)) end if - case('kappa_fission') - if (allocated(this % xs(t) % k_fission)) then - xs = this % xs(t) % k_fission(gin, iazi, ipol) - else - xs = ZERO - end if - case('nu_fission') - xs = this % xs(t) % nu_fission(gin, iazi, ipol) - case('chi') + + case('nu-fission') + xs = this % xs(t) % prompt_nu_fission(gin, iazi, ipol) + & + sum(this % xs(t) % delayed_nu_fission(:, gin, iazi, ipol)) + + case('chi-prompt') if (present(gout)) then - xs = this % xs(t) % chi(gout, gin, iazi, ipol) + xs = this % xs(t) % chi_prompt(gout, gin, iazi, ipol) else ! Not sure you would want a 1 or a 0, but here you go! - xs = sum(this % xs(t) % chi(:, gin, iazi, ipol)) + xs = sum(this % xs(t) % chi_prompt(:, gin, iazi, ipol)) end if + + case('chi-delayed') + if (present(gout)) then + if (present(dg)) then + xs = this % xs(t) % chi_delayed(dg, gout, gin, iazi, ipol) + else + xs = this % xs(t) % chi_delayed(1, gout, gin, iazi, ipol) + end if + else + if (present(dg)) then + xs = sum(this % xs(t) % chi_delayed(dg, :, gin, iazi, ipol)) + else + xs = sum(this % xs(t) % chi_delayed(1, :, gin, iazi, ipol)) + end if + end if + + case('decay rate') + if (present(dg)) then + xs = this % xs(t) % decay_rate(iazi, ipol, dg) + else + xs = this % xs(t) % decay_rate(iazi, ipol, 1) + end if + + case('inverse-velocity') + xs = this % xs(t) % inverse_velocity(gin, iazi, ipol) + case('scatter') if (present(gout)) then if (gout < this % xs(t) % scatter(iazi, ipol) % obj % gmin(gin) .or. & @@ -1913,6 +3032,7 @@ module mgxs_header else xs = this % xs(t) % scatter(iazi, ipol) % obj % scattxs(gin) end if + case('scatter/mult') if (present(gout)) then if (gout < this % xs(t) % scatter(iazi, ipol) % obj % gmin(gin) .or. & @@ -1928,6 +3048,7 @@ module mgxs_header (dot_product(this % xs(t) % scatter(iazi, ipol) % obj % mult(gin) % data, & this % xs(t) % scatter(iazi, ipol) % obj % energy(gin) % data)) end if + case('scatter*f_mu/mult','scatter*f_mu') if (present(gout)) then if (gout < this % xs(t) % scatter(iazi, ipol) % obj % gmin(gin) .or. & @@ -1949,60 +3070,147 @@ module mgxs_header ! user of this code wants the complete 1-outgoing group distribution ! which Im not sure what they would do with that. end if - case('inv_vel') - xs = this % xs(t) % inv_vel(gin, iazi, ipol) + case default xs = ZERO + end select + else xs = ZERO end if + end function mgxsang_get_xs !=============================================================================== ! MGXS*_SAMPLE_FISSION_ENERGY samples the outgoing energy from a fission event !=============================================================================== - function mgxsiso_sample_fission_energy(this, gin, uvw) result(gout) - class(MgxsIso), intent(in) :: this ! Data to work with - integer, intent(in) :: gin ! Incoming energy group - real(8), intent(in) :: uvw(3) ! Particle Direction - integer :: gout ! Sampled outgoing group - real(8) :: xi ! Our random number - real(8) :: prob ! Running probability + subroutine mgxsiso_sample_fission_energy(this, gin, uvw, dg, gout) - xi = prn() - gout = 1 - prob = this % xs(this % index_temp) % chi(gout,gin) + class(MgxsIso), intent(in) :: this ! Data to work with + integer, intent(in) :: gin ! Incoming energy group + real(8), intent(in) :: uvw(3) ! Particle Direction + integer, intent(out) :: dg ! Delayed group + integer, intent(out) :: gout ! Sampled outgoing group + real(8) :: xi_pd ! Our random number for prompt/delayed + real(8) :: xi_gout ! Our random number for gout + real(8) :: prob_gout ! Running probability for gout - do while (prob < xi) - gout = gout + 1 - prob = prob + this % xs(this % index_temp) % chi(gout,gin) - end do + ! Get nu and nu_prompt + real(8) :: prob_prompt - end function mgxsiso_sample_fission_energy + prob_prompt = this % get_xs('prompt-nu-fission', gin) / & + this % get_xs('nu-fission', gin) - function mgxsang_sample_fission_energy(this, gin, uvw) result(gout) + ! Sample random numbers + xi_pd = prn() + xi_gout = prn() + + ! Neutron is born prompt + if (xi_pd <= prob_prompt) then + + ! set the delayed group for the particle born from fission to 0 + dg = 0 + + gout = 1 + prob_gout = this % get_xs('chi-prompt', gin, gout) + + do while (prob_gout < xi_gout) + gout = gout + 1 + prob_gout = prob_gout + this % get_xs('chi-prompt', gin, gout) + end do + + ! Neutron is born delayed + else + + ! Get the delayed group + dg = 0 + + do while (xi_pd >= prob_prompt) + dg = dg + 1 + prob_prompt = prob_prompt + & + this % get_xs('delayed-nu-fission', gin, dg=dg) & + / this % get_xs('nu-fission', gin) + end do + + ! Adjust dg in case of round off error + dg = min(dg, this % num_delayed_groups) + + ! Get the outgoing group + gout = 1 + prob_gout = this % get_xs('chi-delayed', gin, gout, dg=dg) + + do while (prob_gout < xi_gout) + gout = gout + 1 + prob_gout = prob_gout + this % get_xs('chi-delayed', gin, gout, dg=dg) + end do + end if + + end subroutine mgxsiso_sample_fission_energy + + subroutine mgxsang_sample_fission_energy(this, gin, uvw, dg, gout) class(MgxsAngle), intent(in) :: this ! Data to work with integer, intent(in) :: gin ! Incoming energy group - real(8), intent(in) :: uvw(3) ! Particle Direction - integer :: gout ! Sampled outgoing group - real(8) :: xi ! Our random number - real(8) :: prob ! Running probability - integer :: iazi, ipol ! Angle indices + real(8), intent(in) :: uvw(3) ! Direction vector + integer, intent(out) :: dg ! Delayed group + integer, intent(out) :: gout ! Sampled outgoing group + real(8) :: xi_pd ! Our random number for prompt/delayed + real(8) :: xi_gout ! Our random number for gout + real(8) :: prob_gout ! Running probability for gout + real(8) :: prob_prompt - call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol) + ! Get nu and nu_prompt + prob_prompt = this % get_xs('prompt-nu-fission', gin, uvw=uvw) / & + this % get_xs('nu-fission', gin, uvw=uvw) - xi = prn() - gout = 1 - prob = this % xs(this % index_temp) % chi(gout, gin, iazi, ipol) + ! Sample random numbers + xi_pd = prn() + xi_gout = prn() - do while (prob < xi) - gout = gout + 1 - prob = prob + this % xs(this % index_temp) % chi(gout, gin, iazi, ipol) - end do + ! Neutron is born prompt + if (xi_pd <= prob_prompt) then - end function mgxsang_sample_fission_energy + ! set the delayed group for the particle born from fission to 0 + dg = 0 + + gout = 1 + prob_gout = this % get_xs('chi-prompt', gin, gout, uvw=uvw) + + do while (prob_gout < xi_gout) + gout = gout + 1 + prob_gout = prob_gout + & + this % get_xs('chi-prompt', gin, gout, uvw=uvw) + end do + + ! Neutron is born delayed + else + + ! Get the delayed group + dg = 0 + + do while (xi_pd >= prob_prompt) + dg = dg + 1 + prob_prompt = prob_prompt + & + this % get_xs('delayed-nu-fission', gin, uvw=uvw, dg=dg) / & + this % get_xs('nu-fission', gin, uvw=uvw) + end do + + ! Adjust dg in case of round off error + dg = min(dg, this % num_delayed_groups) + + ! Get the outgoing group + gout = 1 + prob_gout = this % get_xs('chi-delayed', gin, gout, uvw=uvw, dg=dg) + + do while (prob_gout < xi_gout) + gout = gout + 1 + prob_gout = prob_gout + & + this % get_xs('chi-delayed', gin, gout, uvw=uvw, dg=dg) + end do + end if + + end subroutine mgxsang_sample_fission_energy !=============================================================================== ! MGXS*_SAMPLE_SCATTER Selects outgoing energy and angle after a scatter event @@ -2051,7 +3259,9 @@ module mgxs_header xs % elastic = this % xs(this % index_temp) % scatter % scattxs(gin) xs % absorption = this % xs(this % index_temp) % absorption(gin) xs % fission = this % xs(this % index_temp) % fission(gin) - xs % nu_fission = this % xs(this % index_temp) % nu_fission(gin) + xs % nu_fission = & + this % xs(this % index_temp) % prompt_nu_fission(gin) + & + sum(this % xs(this % index_temp) % delayed_nu_fission(:, gin)) end subroutine mgxsiso_calculate_xs @@ -2073,7 +3283,9 @@ module mgxs_header xs % fission = this % xs(this % index_temp) % & fission(gin, iazi, ipol) xs % nu_fission = this % xs(this % index_temp) % & - nu_fission(gin, iazi, ipol) + prompt_nu_fission(gin, iazi, ipol) + & + sum(this % xs(this % index_temp) % & + delayed_nu_fission(:, gin, iazi, ipol)) end subroutine mgxsang_calculate_xs diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index c60fc6e2e1..776c7905aa 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -171,7 +171,10 @@ contains type(Bank), intent(inout) :: bank_array(:) integer(8), intent(inout) :: size_bank + integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born integer :: i ! loop index + integer :: dg ! delayed group + integer :: gout ! group out integer :: nu ! actual number of neutrons produced integer :: ijk(3) ! indices in ufs mesh real(8) :: nu_t ! total nu @@ -190,8 +193,10 @@ contains ! the expected number of fission sites produced if (ufs) then + ! Determine indices on ufs mesh for current location call get_mesh_indices(ufs_mesh, p % coord(1) % xyz, ijk, in_mesh) + if (.not. in_mesh) then call write_particle_restart(p) call fatal_error("Source site outside UFS mesh!") @@ -209,6 +214,7 @@ contains ! Determine expected number of neutrons produced nu_t = p % wgt / keff * weight * & material_xs % nu_fission / material_xs % total + ! Sample number of neutrons produced if (prn() > nu_t - int(nu_t)) then nu = int(nu_t) @@ -234,6 +240,10 @@ contains ! Bank source neutrons if (nu == 0 .or. size_bank == size(bank_array)) return + ! Initialize counter of delayed neutrons encountered for each delayed group + ! to zero. + nu_d(:) = 0 + p % fission = .true. ! Fission neutrons will be banked do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4) ! Bank source neutrons by copying particle data @@ -254,8 +264,19 @@ contains ! Sample secondary energy distribution for fission reaction and set energy ! in fission bank - bank_array(i) % E = & - real(xs % sample_fission_energy(p % g, bank_array(i) % uvw), 8) + call xs % sample_fission_energy(p % g, bank_array(i) % uvw, dg, gout) + + bank_array(i) % E = real(gout, 8) + bank_array(i) % delayed_group = dg + + ! Set delayed group on particle too + p % delayed_group = dg + + ! Increment the number of neutrons born delayed + if (p % delayed_group > 0) then + nu_d(p % delayed_group) = nu_d(p % delayed_group) + 1 + end if + end do ! increment number of bank sites @@ -264,6 +285,7 @@ contains ! Store total weight banked for analog fission tallies p % n_bank = nu p % wgt_bank = nu/weight + p % n_delayed_bank(:) = nu_d(:) end subroutine create_fission_sites diff --git a/src/source.F90 b/src/source.F90 index 9aeccde156..cc85bf6eca 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -199,10 +199,10 @@ contains if (.not. run_CE) then if (site % E <= energy_bins(1)) then site % E = real(1, 8) - else if (site % E > energy_bins(energy_groups + 1)) then - site % E = real(energy_groups, 8) + else if (site % E > energy_bins(num_energy_groups + 1)) then + site % E = real(num_energy_groups, 8) else - site % E = real(binary_search(energy_bins, energy_groups + 1, & + site % E = real(binary_search(energy_bins, num_energy_groups + 1, & site % E), 8) end if end if diff --git a/src/tally.F90 b/src/tally.F90 index 0e0c6dd7e8..69e4e8ecdc 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -164,6 +164,7 @@ contains case (SCORE_INVERSE_VELOCITY) + ! make sure the correct energy is used if (t % estimator == ESTIMATOR_TRACKLENGTH) then E = p % E @@ -382,7 +383,7 @@ contains ! neutrons were emitted with different energies, multiple ! outgoing energy bins may have been scored to. The following ! logic treats this special case and results to multiple bins - call score_fission_eout_ce(p, t, score_index, score_bin) + call score_fission_eout(p, t, score_index, score_bin) cycle SCORE_LOOP end if end if @@ -432,7 +433,7 @@ contains ! neutrons were emitted with different energies, multiple ! outgoing energy bins may have been scored to. The following ! logic treats this special case and results to multiple bins - call score_fission_eout_ce(p, t, score_index, score_bin) + call score_fission_eout(p, t, score_index, score_bin) cycle SCORE_LOOP end if end if @@ -504,7 +505,7 @@ contains ! neutrons were emitted with different energies, multiple ! outgoing energy bins may have been scored to. The following ! logic treats this special case and results to multiple bins - call score_fission_eout_ce(p, t, score_index, score_bin) + call score_fission_eout(p, t, score_index, score_bin) cycle SCORE_LOOP end if end if @@ -512,7 +513,8 @@ contains ! No fission events occur if survival biasing is on -- need to ! calculate fraction of absorptions that would have resulted in ! delayed-nu-fission - if (micro_xs(p % event_nuclide) % absorption > ZERO) then + if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & + nuclides(p % event_nuclide) % fissionable) then ! Check if the delayed group filter is present if (dg_filter > 0) then @@ -572,7 +574,9 @@ contains d = filt % groups(d_bin) ! Compute the score and tally to bin - score = keff * p % wgt_bank / p % n_bank * p % n_delayed_bank(d) + score = keff * p % wgt_bank / p % n_bank * & + p % n_delayed_bank(d) + call score_fission_delayed_dg(t, d_bin, score, score_index) end do cycle SCORE_LOOP @@ -1104,6 +1108,11 @@ contains integer :: q ! loop index for scoring bins integer :: score_bin ! scoring bin, e.g. SCORE_FLUX integer :: score_index ! scoring bin index + integer :: d ! delayed neutron index + integer :: g ! delayed neutron index + integer :: k ! loop index for bank sites + integer :: d_bin ! delayed group bin index + integer :: dg_filter ! index of delayed group filter real(8) :: score ! analog tally score real(8) :: p_uvw(3) ! Particle's current uvw integer :: p_g ! Particle group to use for getting info @@ -1115,7 +1124,9 @@ contains ! this only depends on if we if (t % estimator == ESTIMATOR_ANALOG .or. & t % estimator == ESTIMATOR_COLLISION) then + if (survival_biasing) then + ! Then we either are alive and had a scatter (and so g changed), ! or are dead and g did not change if (p % alive) then @@ -1126,16 +1137,19 @@ contains p_g = p % g end if else if (p % event == EVENT_SCATTER) then + ! Then the energy group has been changed by the scattering routine ! meaning gin is now in p % last_g p_uvw = p % last_uvw p_g = p % last_g else + ! No scatter, no change in g. p_uvw = p % coord(p % n_coord) % uvw p_g = p % g end if else + ! No actual collision so g has not changed. p_uvw = p % coord(p % n_coord) % uvw p_g = p % g @@ -1144,6 +1158,7 @@ contains ! To significantly reduce de-referencing, point matxs to the ! macroscopic Mgxs for the material of interest matxs => macro_xs(p % material) % obj + ! Do same for nucxs, point it to the microscopic nuclide data of interest if (i_nuclide > 0) then nucxs => nuclides_MG(i_nuclide) % obj @@ -1172,6 +1187,7 @@ contains ! All events score to a flux bin. We actually use a collision ! estimator in place of an analog one since there is no way to count ! 'events' exactly for the flux + if (survival_biasing) then ! We need to account for the fact that some weight was already ! absorbed @@ -1199,6 +1215,7 @@ contains else score = p % last_wgt end if + if (i_nuclide > 0) then score = score * atom_density * & nucxs % get_xs('total', p_g, UVW=p_uvw) / & @@ -1228,20 +1245,21 @@ contains else score = p % last_wgt end if + if (i_nuclide > 0) then - score = score * nucxs % get_xs('inv_vel', p_g, UVW=p_uvw) / & - matxs % get_xs('total', p_g, UVW=p_uvw) * flux + score = score * nucxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) & + / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux else - score = matxs % get_xs('inv_vel', p_g, UVW=p_uvw) * flux + score = score * matxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) & + / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux end if else - ! For inverse velocity, we need no cross section + if (i_nuclide > 0) then - score = score * nucxs % get_xs('inv_vel', p_g, UVW=p_uvw) * & - atom_density * flux + score = flux * nucxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) else - score = flux * matxs % get_xs('inv_vel', p_g, UVW=p_uvw) + score = flux * matxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) end if end if @@ -1358,6 +1376,14 @@ contains case (SCORE_FISSION) + + ! make sure the correct energy is used + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + p_g = p % g + else + p_g = p % last_g + end if + if (t % estimator == ESTIMATOR_ANALOG) then if (survival_biasing) then ! No fission events occur if survival biasing is on -- need to @@ -1386,13 +1412,20 @@ contains score = nucxs % get_xs('fission', p_g, UVW=p_uvw) * & atom_density * flux else - score = flux * material_xs % fission - + score = matxs % get_xs('fission', p_g, UVW=p_uvw) * flux end if end if case (SCORE_NU_FISSION) + + ! make sure the correct energy is used + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + p_g = p % g + else + p_g = p % last_g + end if + if (t % estimator == ESTIMATOR_ANALOG) then if (survival_biasing .or. p % fission) then if (t % find_filter(FILTER_ENERGYOUT) > 0) then @@ -1401,8 +1434,7 @@ contains ! neutrons were emitted with different energies, multiple ! outgoing energy bins may have been scored to. The following ! logic treats this special case and results to multiple bins - call score_fission_eout_mg(p, t, score_index, i_nuclide, & - atom_density) + call score_fission_eout(p, t, score_index, score_bin) cycle SCORE_LOOP end if end if @@ -1413,11 +1445,11 @@ contains score = p % absorb_wgt * flux if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('nu_fission', p_g, UVW=p_uvw) / & + nucxs % get_xs('nu-fission', p_g, UVW=p_uvw) / & matxs % get_xs('absorption', p_g, UVW=p_uvw) else score = score * & - matxs % get_xs('nu_fission', p_g, UVW=p_uvw) / & + matxs % get_xs('nu-fission', p_g, UVW=p_uvw) / & matxs % get_xs('absorption', p_g, UVW=p_uvw) end if else @@ -1438,15 +1470,421 @@ contains else if (i_nuclide > 0) then - score = nucxs % get_xs('nu_fission', p_g, UVW=p_uvw) * & + score = nucxs % get_xs('nu-fission', p_g, UVW=p_uvw) * & atom_density * flux else - score = material_xs % nu_fission * flux + score = matxs % get_xs('nu-fission', p_g, UVW=p_uvw) * flux + end if + end if + + + case (SCORE_PROMPT_NU_FISSION) + + ! make sure the correct energy is used + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + p_g = p % g + else + p_g = p % last_g + end if + + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing .or. p % fission) then + if (t % find_filter(FILTER_ENERGYOUT) > 0) then + ! Normally, we only need to make contributions to one scoring + ! bin. However, in the case of fission, since multiple fission + ! neutrons were emitted with different energies, multiple + ! outgoing energy bins may have been scored to. The following + ! logic treats this special case and results to multiple bins + call score_fission_eout(p, t, score_index, score_bin) + cycle SCORE_LOOP + end if + end if + if (survival_biasing) then + ! No fission events occur if survival biasing is on -- need to + ! calculate fraction of absorptions that would have resulted in + ! nu-fission + score = p % absorb_wgt * flux + if (i_nuclide > 0) then + score = score * atom_density * & + nucxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) / & + matxs % get_xs('absorption', p_g, UVW=p_uvw) + else + score = score * & + matxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) / & + matxs % get_xs('absorption', p_g, UVW=p_uvw) + end if + else + ! Skip any non-fission events + if (.not. p % fission) cycle SCORE_LOOP + ! If there is no outgoing energy filter, than we only need to + ! score to one bin. For the score to be 'analog', we need to + ! score the number of particles that were banked in the fission + ! bank. Since this was weighted by 1/keff, we multiply by keff + ! to get the proper score. + score = keff * p % wgt_bank * (ONE - sum(p % n_delayed_bank) & + / real(p % n_bank, 8)) + if (i_nuclide > 0) then + score = score * atom_density * & + nucxs % get_xs('fission', p_g, UVW=p_uvw) / & + matxs % get_xs('fission', p_g, UVW=p_uvw) + end if + end if + + else + if (i_nuclide > 0) then + score = nucxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) * & + atom_density * flux + else + score = matxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) * flux + end if + end if + + + case (SCORE_DELAYED_NU_FISSION) + + ! make sure the correct energy is used + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + p_g = p % g + else + p_g = p % last_g + end if + + ! Set the delayedgroup filter index and the number of delayed group bins + dg_filter = t % find_filter(FILTER_DELAYEDGROUP) + + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing .or. p % fission) then + if (t % find_filter(FILTER_ENERGYOUT) > 0) then + ! Normally, we only need to make contributions to one scoring + ! bin. However, in the case of fission, since multiple fission + ! neutrons were emitted with different energies, multiple + ! outgoing energy bins may have been scored to. The following + ! logic treats this special case and results to multiple bins + call score_fission_eout(p, t, score_index, score_bin) + cycle SCORE_LOOP + end if + end if + if (survival_biasing) then + ! No fission events occur if survival biasing is on -- need to + ! calculate fraction of absorptions that would have resulted in + ! nu-fission + if (matxs % get_xs('absorption', p_g, UVW=p_uvw) > ZERO) then + + if (dg_filter > 0) then + select type(filt => t % filters(dg_filter) % obj) + type is (DelayedGroupFilter) + + ! Loop over all delayed group bins and tally to them + ! individually + do d_bin = 1, filt % n_bins + + ! Get the delayed group for this bin + d = filt % groups(d_bin) + + score = p % absorb_wgt * flux + if (i_nuclide > 0) then + score = score * nucxs % get_xs('delayed-nu-fission', & + p_g, UVW=p_uvw, dg=d) / & + matxs % get_xs('absorption', p_g, UVW=p_uvw) + else + score = score * matxs % get_xs('delayed-nu-fission', & + p_g, UVW=p_uvw, dg=d) / & + matxs % get_xs('absorption', p_g, UVW=p_uvw) + end if + + call score_fission_delayed_dg(t, d_bin, score, score_index) + end do + cycle SCORE_LOOP + end select + else + score = p % absorb_wgt * flux + if (i_nuclide > 0) then + score = score * nucxs % get_xs('delayed-nu-fission', p_g, & + UVW=p_uvw) / matxs % get_xs('absorption', p_g, UVW=p_uvw) + else + score = score * matxs % get_xs('delayed-nu-fission', p_g, & + UVW=p_uvw) / matxs % get_xs('absorption', p_g, UVW=p_uvw) + end if + end if + end if + else + ! Skip any non-fission events + if (.not. p % fission) cycle SCORE_LOOP + ! If there is no outgoing energy filter, than we only need to + ! score to one bin. For the score to be 'analog', we need to + ! score the number of particles that were banked in the fission + ! bank. Since this was weighted by 1/keff, we multiply by keff + ! to get the proper score. + + ! Check if the delayed group filter is present + if (dg_filter > 0) then + select type(filt => t % filters(dg_filter) % obj) + type is (DelayedGroupFilter) + + ! Loop over all delayed group bins and tally to them + ! individually + do d_bin = 1, filt % n_bins + + ! Get the delayed group for this bin + d = filt % groups(d_bin) + + score = keff * p % wgt_bank / p % n_bank * & + p % n_delayed_bank(d) + + if (i_nuclide > 0) then + score = score * atom_density * & + nucxs % get_xs('fission', p_g, UVW=p_uvw) / & + matxs % get_xs('fission', p_g, UVW=p_uvw) + end if + + call score_fission_delayed_dg(t, d_bin, score, score_index) + end do + cycle SCORE_LOOP + end select + else + score = keff * p % wgt_bank / p % n_bank * sum(p % n_delayed_bank) + if (i_nuclide > 0) then + score = score * atom_density * & + nucxs % get_xs('fission', p_g, UVW=p_uvw) / & + matxs % get_xs('fission', p_g, UVW=p_uvw) + end if + end if + end if + else + + ! Check if the delayed group filter is present + if (dg_filter > 0) then + select type(filt => t % filters(dg_filter) % obj) + type is (DelayedGroupFilter) + + ! Loop over all delayed group bins and tally to them + ! individually + do d_bin = 1, filt % n_bins + + ! Get the delayed group for this bin + d = filt % groups(d_bin) + + if (i_nuclide > 0) then + score = nucxs % get_xs('delayed-nu-fission', p_g, & + UVW=p_uvw, dg=d) * atom_density * flux + else + score = matxs % get_xs('delayed-nu-fission', p_g, & + UVW=p_uvw, dg=d) * flux + end if + + call score_fission_delayed_dg(t, d_bin, score, score_index) + end do + cycle SCORE_LOOP + end select + else + if (i_nuclide > 0) then + score = nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw) & + * atom_density * flux + else + score = matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw) & + * flux + end if + end if + end if + + case (SCORE_DECAY_RATE) + + ! make sure the correct energy is used + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + p_g = p % g + else + p_g = p % last_g + end if + + ! Set the delayedgroup filter index and the number of delayed group bins + dg_filter = t % find_filter(FILTER_DELAYEDGROUP) + + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing) then + ! No fission events occur if survival biasing is on -- need to + ! calculate fraction of absorptions that would have resulted in + ! nu-fission + if (matxs % get_xs('absorption', p_g, UVW=p_uvw) > ZERO) then + + if (dg_filter > 0) then + select type(filt => t % filters(dg_filter) % obj) + type is (DelayedGroupFilter) + + ! Loop over all delayed group bins and tally to them + ! individually + do d_bin = 1, filt % n_bins + + ! Get the delayed group for this bin + d = filt % groups(d_bin) + + score = p % absorb_wgt * flux + if (i_nuclide > 0) then + score = score * nucxs % get_xs('decay rate', p_g, & + UVW=p_uvw, dg=d) * & + nucxs % get_xs('delayed-nu-fission', p_g, & + UVW=p_uvw, dg=d) / matxs % get_xs('absorption', & + p_g, UVW=p_uvw) + else + score = score * matxs % get_xs('decay rate', p_g, & + UVW=p_uvw, dg=d) * & + matxs % get_xs('delayed-nu-fission', p_g, & + UVW=p_uvw, dg=d) / matxs % get_xs('absorption', & + p_g, UVW=p_uvw) + end if + + call score_fission_delayed_dg(t, d_bin, score, score_index) + end do + cycle SCORE_LOOP + end select + else + + score = ZERO + + ! If the delayed group filter is not present, compute the score + ! by accumulating the absorbed weight times the decay rate times + ! the fraction of the delayed-nu-fission xs to the absorption xs + ! for all delayed groups. + do d = 1, num_delayed_groups + if (i_nuclide > 0) then + score = score + p % absorb_wgt * & + nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & + nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, & + dg=d) / matxs % get_xs('absorption', p_g, UVW=p_uvw) + else + score = score + p % absorb_wgt * & + matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & + matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, & + dg=d) / matxs % get_xs('absorption', p_g, UVW=p_uvw) + end if + end do + end if + end if + else + ! Skip any non-fission events + if (.not. p % fission) cycle SCORE_LOOP + ! If there is no outgoing energy filter, than we only need to + ! score to one bin. For the score to be 'analog', we need to + ! score the number of particles that were banked in the fission + ! bank. Since this was weighted by 1/keff, we multiply by keff + ! to get the proper score. + + score = ZERO + + ! loop over number of particles banked + do k = 1, p % n_bank + + ! get the delayed group + g = fission_bank(n_bank - p % n_bank + k) % delayed_group + + ! Case for tallying delayed emissions + if (g /= 0) then + + ! determine score based on bank site weight and keff. + if (i_nuclide > 0) then + score = score + keff * atom_density * & + fission_bank(n_bank - p % n_bank + k) % wgt * & + nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=g) * & + nucxs % get_xs('fission', p_g, UVW=p_uvw) / & + matxs % get_xs('fission', p_g, UVW=p_uvw) + else + score = score + keff * & + fission_bank(n_bank - p % n_bank + k) % wgt * & + matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=g) + end if + + ! if the delayed group filter is present, tally to corresponding + ! delayed group bin if it exists + if (dg_filter > 0) then + + ! declare the delayed group filter type + select type(filt => t % filters(dg_filter) % obj) + type is (DelayedGroupFilter) + + ! loop over delayed group bins until the corresponding bin + ! is found + do d_bin = 1, filt % n_bins + d = filt % groups(d_bin) + + ! check whether the delayed group of the particle is equal + ! to the delayed group of this bin + if (d == g) then + call score_fission_delayed_dg(t, d_bin, score, & + score_index) + end if + end do + end select + + ! Reset the score to zero + score = ZERO + end if + end if + end do + + ! If the delayed group filter is present, cycle because the + ! score_fission_delayed_dg(...) has already tallied the score + if (dg_filter > 0) then + cycle SCORE_LOOP + end if + end if + else + + ! Check if the delayed group filter is present + if (dg_filter > 0) then + select type(filt => t % filters(dg_filter) % obj) + type is (DelayedGroupFilter) + + ! Loop over all delayed group bins and tally to them + ! individually + do d_bin = 1, filt % n_bins + + ! Get the delayed group for this bin + d = filt % groups(d_bin) + + if (i_nuclide > 0) then + score = nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & + nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, & + dg=d) * atom_density * flux + else + score = matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & + matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, & + dg=d) * flux + end if + + call score_fission_delayed_dg(t, d_bin, score, score_index) + end do + cycle SCORE_LOOP + end select + else + score = ZERO + + ! If the delayed group filter is not present, compute the score + ! by accumulating the absorbed weight times the decay rate times + ! the fraction of the delayed-nu-fission xs to the absorption xs + ! for all delayed groups. + do d = 1, num_delayed_groups + if (i_nuclide > 0) then + score = score + atom_density * flux * & + nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & + nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, dg=d) + else + score = score + flux * & + matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & + matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, dg=d) + end if + end do end if end if case (SCORE_KAPPA_FISSION) + + ! make sure the correct energy is used + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + p_g = p % g + else + p_g = p % last_g + end if + if (t % estimator == ESTIMATOR_ANALOG) then if (survival_biasing) then ! No fission events occur if survival biasing is on -- need to @@ -1463,19 +1901,19 @@ contains end if if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('kappa_fission', p_g, UVW=p_uvw) / & + nucxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / & matxs % get_xs('absorption', p_g, UVW=p_uvw) else score = score * & - matxs % get_xs('kappa_fission', p_g, UVW=p_uvw) / & + matxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / & matxs % get_xs('absorption', p_g, UVW=p_uvw) end if else if (i_nuclide > 0) then - score = flux * atom_density * & - nucxs % get_xs('kappa_fission', p_g, UVW=p_uvw) + score = nucxs % get_xs('kappa-fission', p_g, UVW=p_uvw) * & + atom_density * flux else - score = flux * matxs % get_xs('kappa_fission', p_g, UVW=p_uvw) + score = matxs % get_xs('kappa-fission', p_g, UVW=p_uvw) * flux end if end if @@ -1937,7 +2375,7 @@ contains ! neutrons produced with different energies. !=============================================================================== - subroutine score_fission_eout_ce(p, t, i_score, score_bin) + subroutine score_fission_eout(p, t, i_score, score_bin) type(Particle), intent(in) :: p type(TallyObject), intent(inout) :: t @@ -1956,6 +2394,7 @@ contains real(8) :: filter_weight ! combined weight of all filters real(8) :: score ! actual score real(8) :: E_out ! energy of fission bank site + integer :: g_out ! energy group of fission bank site ! save original outgoing energy bin and score index i = t % find_filter(FILTER_ENERGYOUT) @@ -1982,14 +2421,31 @@ contains ! determine score based on bank site weight and keff score = keff * fission_bank(n_bank - p % n_bank + k) % wgt - ! determine outgoing energy from fission bank - E_out = fission_bank(n_bank - p % n_bank + k) % E + if (.not. run_CE .and. eo_filt % matches_transport_groups) then - ! check if outgoing energy is within specified range on filter - if (E_out < eo_filt % bins(1) .or. E_out > eo_filt % bins(n)) cycle + ! determine outgoing energy from fission bank + g_out = int(fission_bank(n_bank - p % n_bank + k) % E) - ! change outgoing energy bin - Matching_bins(i) = binary_search(eo_filt % bins, n, E_out) + ! change outgoing energy bin + matching_bins(i) = g_out + + else + + ! determine outgoing energy from fission bank + if (run_CE) then + E_out = fission_bank(n_bank - p % n_bank + k) % E + else + E_out = energy_bin_avg(int(fission_bank(n_bank - p % n_bank + k) & + % E)) + end if + + ! check if outgoing energy is within specified range on filter + if (E_out < eo_filt % bins(1) .or. E_out > eo_filt % bins(n)) cycle + + ! change outgoing energy bin + matching_bins(i) = binary_search(eo_filt % bins, n, E_out) + + end if ! Case for tallying prompt neutrons if (score_bin == SCORE_NU_FISSION .or. & @@ -2052,91 +2508,7 @@ contains ! reset outgoing energy bin and score index matching_bins(i) = bin_energyout - end subroutine score_fission_eout_ce - - subroutine score_fission_eout_mg(p, t, i_score, i_nuclide, atom_density) - type(Particle), intent(in) :: p - type(TallyObject), intent(inout) :: t - integer, intent(in) :: i_score ! index for score - integer, intent(in) :: i_nuclide ! index for nuclide - real(8), intent(in) :: atom_density - - integer :: i ! index of outgoing energy filter - integer :: n ! number of energies on filter - integer :: k ! loop index for bank sites - integer :: bin_energyout ! original outgoing energy bin - integer :: i_filter ! index for matching filter bin combination - real(8) :: filter_weight ! combined weight of all filters - real(8) :: score ! actual score - integer :: gout ! energy group of fission bank site - integer :: gin ! energy group of incident particle - real(8) :: E_out - - ! save original outgoing energy bin and score index - i = t % find_filter(FILTER_ENERGYOUT) - bin_energyout = matching_bins(i) - - ! Declare the filter type - select type(filt => t % filters(i) % obj) - type is (EnergyoutFilter) - - ! Get number of energies on filter - n = size(filt % bins) - - ! Since the creation of fission sites is weighted such that it is - ! expected to create n_particles sites, we need to multiply the - ! score by keff to get the true nu-fission rate. Otherwise, the sum - ! of all nu-fission rates would be ~1.0. - - ! loop over number of particles banked - do k = 1, p % n_bank - ! determine score based on bank site weight and keff - score = keff * fission_bank(n_bank - p % n_bank + k) % wgt - if (i_nuclide > 0) then - if (survival_biasing) then - gin = p % g - else - gin = p % last_g - end if - score = score * atom_density * & - nuclides_MG(i_nuclide) % obj % get_xs('fission', gin, & - UVW=p % last_uvw) / & - macro_xs(p % material) % obj % get_xs('fission', gin, & - UVW=p % last_uvw) - end if - - if (filt % matches_transport_groups) then - ! determine outgoing energy from fission bank - gout = int(fission_bank(n_bank - p % n_bank + k) % E) - - ! change outgoing energy bin - matching_bins(i) = gout - else - ! determine outgoing energy from fission bank - E_out = energy_bin_avg(int(fission_bank(n_bank - p % n_bank + k) % E)) - - ! check if outgoing energy is within specified range on filter - if (E_out < filt % bins(1) .or. E_out > filt % bins(n)) cycle - - ! change outgoing energy bin - matching_bins(i) = binary_search(filt % bins, n, E_out) - end if - - ! determine scoring index and weight for this filter combination - i_filter = sum((matching_bins(1:size(t%filters)) - 1) * t % stride) + 1 - filter_weight = product(filter_weights(:size(t % filters))) - - ! Add score to tally -!$omp atomic - t % results(i_score, i_filter) % value = & - t % results(i_score, i_filter) % value + score * filter_weight - end do - - ! reset outgoing energy bin and score index - matching_bins(i) = bin_energyout - end select - - end subroutine score_fission_eout_mg + end subroutine score_fission_eout !=============================================================================== ! SCORE_FISSION_DELAYED_DG helper function used to increment the tally when a diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 4dd9bf129d..37265cc5ff 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -890,7 +890,7 @@ contains ! Tallies are ordered in increasing groups, group indices ! however are the opposite, so switch - next_bin = energy_groups - next_bin + 1 + next_bin = num_energy_groups - next_bin + 1 else ! Make sure the correct energy is used. @@ -959,7 +959,7 @@ contains ! Tallies are ordered in increasing groups, group indices ! however are the opposite, so switch - next_bin = energy_groups - next_bin + 1 + next_bin = num_energy_groups - next_bin + 1 else ! Check if energy of the particle is within energy bins. diff --git a/tests/1d_mgxs.h5 b/tests/1d_mgxs.h5 index d77c19c119..e3b85c2416 100644 Binary files a/tests/1d_mgxs.h5 and b/tests/1d_mgxs.h5 differ diff --git a/tests/test_mg_basic/results_true.dat b/tests/test_mg_basic/results_true.dat index 01b7925469..2192c3746d 100644 --- a/tests/test_mg_basic/results_true.dat +++ b/tests/test_mg_basic/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.160653E+00 6.227102E-02 +1.086852E+00 2.677280E-02 diff --git a/tests/test_mg_max_order/results_true.dat b/tests/test_mg_max_order/results_true.dat index 0e470e2cdb..2189587d4f 100644 --- a/tests/test_mg_max_order/results_true.dat +++ b/tests/test_mg_max_order/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.098323E+00 1.979833E-02 +1.077247E+00 2.412834E-02 diff --git a/tests/test_mg_nuclide/results_true.dat b/tests/test_mg_nuclide/results_true.dat index 01b7925469..2192c3746d 100644 --- a/tests/test_mg_nuclide/results_true.dat +++ b/tests/test_mg_nuclide/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.160653E+00 6.227102E-02 +1.086852E+00 2.677280E-02 diff --git a/tests/test_mg_tallies/results_true.dat b/tests/test_mg_tallies/results_true.dat index b540d9a492..a765c0ba02 100644 --- a/tests/test_mg_tallies/results_true.dat +++ b/tests/test_mg_tallies/results_true.dat @@ -1 +1 @@ -6ab4a178b642f886eea9fddbdf67a8d6d3e0d1968f06ee52d7fd5ecdc1b4a796398e9b3a8c15ac72fb1a011d13a5fcc59a868781d38fb11787e060bb2e602a7b \ No newline at end of file +3d2ce1b8bdd558fe9f8560e8bb91455e1fedd80a47349344399e22f5b0fac1391f1721e7f12a42bb0c458fd6c87ebf9fa38fe9539fcf69292b21cebf8e5f8989 \ No newline at end of file diff --git a/tests/test_mgxs_library_ce_to_mg/results_true.dat b/tests/test_mgxs_library_ce_to_mg/results_true.dat index 0b73724c9c..27db7f3f66 100644 --- a/tests/test_mgxs_library_ce_to_mg/results_true.dat +++ b/tests/test_mgxs_library_ce_to_mg/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.140724E+00 2.952379E-02 +1.133362E+00 1.746521E-02 diff --git a/tests/test_mgxs_library_condense/test_mgxs_library_condense.py b/tests/test_mgxs_library_condense/test_mgxs_library_condense.py index a0d10e6265..0adedac1c6 100644 --- a/tests/test_mgxs_library_condense/test_mgxs_library_condense.py +++ b/tests/test_mgxs_library_condense/test_mgxs_library_condense.py @@ -22,9 +22,6 @@ class MGXSTestHarness(PyAPITestHarness): # Initialize a two-group structure energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 0.625, 20.e6]) - # Initialize a six-delayed-group structure - delayed_groups = list(range(1,7)) - # Initialize MGXS Library for a few cross section types self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry) self.mgxs_lib.by_nuclide = False @@ -33,7 +30,7 @@ class MGXSTestHarness(PyAPITestHarness): self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \ openmc.mgxs.MDGXS_TYPES self.mgxs_lib.energy_groups = energy_groups - self.mgxs_lib.delayed_groups = delayed_groups + self.mgxs_lib.num_delayed_groups = 6 self.mgxs_lib.legendre_order = 3 self.mgxs_lib.domain_type = 'material' self.mgxs_lib.build_library() diff --git a/tests/test_mgxs_library_distribcell/test_mgxs_library_distribcell.py b/tests/test_mgxs_library_distribcell/test_mgxs_library_distribcell.py index a37c1b83fa..cb64d9d76a 100644 --- a/tests/test_mgxs_library_distribcell/test_mgxs_library_distribcell.py +++ b/tests/test_mgxs_library_distribcell/test_mgxs_library_distribcell.py @@ -22,9 +22,6 @@ class MGXSTestHarness(PyAPITestHarness): # Initialize a one-group structure energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 20.e6]) - # Initialize a six-delayed-group structure - delayed_groups = list(range(1,7)) - # Initialize MGXS Library for a few cross section types # for one material-filled cell in the geometry self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry) @@ -34,7 +31,7 @@ class MGXSTestHarness(PyAPITestHarness): self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \ openmc.mgxs.MDGXS_TYPES self.mgxs_lib.energy_groups = energy_groups - self.mgxs_lib.delayed_groups = delayed_groups + self.mgxs_lib.num_delayed_groups = 6 self.mgxs_lib.legendre_order = 3 self.mgxs_lib.domain_type = 'distribcell' cells = self.mgxs_lib.openmc_geometry.get_all_material_cells() diff --git a/tests/test_mgxs_library_hdf5/test_mgxs_library_hdf5.py b/tests/test_mgxs_library_hdf5/test_mgxs_library_hdf5.py index db8dc33201..9197f549c1 100644 --- a/tests/test_mgxs_library_hdf5/test_mgxs_library_hdf5.py +++ b/tests/test_mgxs_library_hdf5/test_mgxs_library_hdf5.py @@ -29,9 +29,6 @@ class MGXSTestHarness(PyAPITestHarness): # Initialize a two-group structure energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 0.625, 20.e6]) - # Initialize a six-delayed-group structure - delayed_groups = list(range(1,7)) - # Initialize MGXS Library for a few cross section types self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry) self.mgxs_lib.by_nuclide = False @@ -40,7 +37,7 @@ class MGXSTestHarness(PyAPITestHarness): self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \ openmc.mgxs.MDGXS_TYPES self.mgxs_lib.energy_groups = energy_groups - self.mgxs_lib.delayed_groups = delayed_groups + self.mgxs_lib.num_delayed_groups = 6 self.mgxs_lib.legendre_order = 3 self.mgxs_lib.domain_type = 'material' self.mgxs_lib.build_library() diff --git a/tests/test_mgxs_library_mesh/test_mgxs_library_mesh.py b/tests/test_mgxs_library_mesh/test_mgxs_library_mesh.py index 5dd1848bae..f3f62d7f1e 100644 --- a/tests/test_mgxs_library_mesh/test_mgxs_library_mesh.py +++ b/tests/test_mgxs_library_mesh/test_mgxs_library_mesh.py @@ -18,9 +18,6 @@ class MGXSTestHarness(PyAPITestHarness): # Initialize a one-group structure energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 20.e6]) - # Initialize a six-delayed-group structure - delayed_groups = list(range(1,7)) - # Initialize MGXS Library for a few cross section types # for one material-filled cell in the geometry self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry) @@ -30,7 +27,7 @@ class MGXSTestHarness(PyAPITestHarness): self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \ openmc.mgxs.MDGXS_TYPES self.mgxs_lib.energy_groups = energy_groups - self.mgxs_lib.delayed_groups = delayed_groups + self.mgxs_lib.num_delayed_groups = 6 self.mgxs_lib.legendre_order = 3 self.mgxs_lib.domain_type = 'mesh' diff --git a/tests/test_mgxs_library_no_nuclides/test_mgxs_library_no_nuclides.py b/tests/test_mgxs_library_no_nuclides/test_mgxs_library_no_nuclides.py index aa4374e608..a5f6d2ccb0 100644 --- a/tests/test_mgxs_library_no_nuclides/test_mgxs_library_no_nuclides.py +++ b/tests/test_mgxs_library_no_nuclides/test_mgxs_library_no_nuclides.py @@ -22,9 +22,6 @@ class MGXSTestHarness(PyAPITestHarness): # Initialize a two-group structure energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 0.625, 20.e6]) - # Initialize a six-delayed-group structure - delayed_groups = list(range(1,7)) - # Initialize MGXS Library for a few cross section types self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry) self.mgxs_lib.by_nuclide = False @@ -33,7 +30,7 @@ class MGXSTestHarness(PyAPITestHarness): self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \ openmc.mgxs.MDGXS_TYPES self.mgxs_lib.energy_groups = energy_groups - self.mgxs_lib.delayed_groups = delayed_groups + self.mgxs_lib.num_delayed_groups = 6 self.mgxs_lib.legendre_order = 3 self.mgxs_lib.domain_type = 'material' self.mgxs_lib.build_library()